First Principle Surface Analysis of YF3 and Isostructural HoF3

The trifluorides of the two high field strength elements yttrium and holmium are studied by periodic density functional theory. As a lanthanide, holmium also belongs to the group of rare earth elements (REE). Due to their equivalent geochemical behavior, both elements form a geochemical twin pair and consequently, yttrium is generally associated with the REE as REE+Y. Interestingly, it has been found that DFT/DFT+U describe bulk HoF3 best, when the 4f-electrons are excluded from the valence region. An extensive surface stability analysis of YF3 (PBE) and HoF3 (PBE+U d /3 eV/4f-in-core) using two-dimensional surface models (slabs) is performed. All seven low-lying Miller indices surfaces are considered with all possible stoichiometric or substoichiometric terminations with a maximal fluorine-deficit of two. This leads to a scope of 24 terminations per compound. The resulting Wulff plots consists of seven surfaces with 5–26% abundance for YF3 and six surfaces with 6–34% for HoF3. The stoichiometric (010) surface is dominating in both compounds. However, subtle differences have been found between these two geochemical twins.


Introduction
Yttrium and holmium form a geochemical twin pair. The term emphasizes their identical geochemical behavior caused by the equal ratio of charge to radius in their only stable oxidation state +III. According to their small ionic radii of 1.075 Å (Y) and 1.072 Å (Ho) in nine-fold coordination [1] and their high oxidation state, both belong to the interesting group of high field strength elements (HFSE). As a lanthanide, holmium also belongs to the rare earth elements (REE). Due to their twin character, yttrium is also often associated with that group [2][3][4].
As fluorides, both metals can be used for different specific applications. The wideband-gap material YF 3 has very good properties for laser applications [5][6][7][8]. Doped with trivalent REE cations, YF 3 is also applicable as an optical filter in 157-nm photolithography [9]. Another emerging field of application is solid-state fluoride batteries, resulting from the very high conductivity of fluoride anions [10][11][12][13][14]. HoF 3 is interesting for magnetic high-field applications as, e.g., a contrast agent, due to the very high magnetic moment of holmium [15,16]. Moreover, YF 3 and HoF 3 are important precursors for the synthesis of the respective pure metallic compounds [17,18]. In nature, YF 3 is found within the mineral waimirite-(Y), which contains high concentrations of other REE [8]. Fluoride plays a significant role in accumulating HFSE and REE within hydrothermal fluids, as these cations do not form such stable complexes with chloride [2,[19][20][21]. Interestingly, those fluoride-rich hydrothermal fluids produce ores with a non-chondritic excess of yttrium over holmium. It is suggested that one underlying reason for the twin separation is their different affinity to fluorine, which was found in dissolving experiments of YF 3 and HoF 3 in diluted hydrofluoric acid [20]. To lay one foundation for future quantum chemical studies on the different fluorine-affinity of yttrium and holmium, we started with an investigation of the respective trifluorides and their surfaces.

YF 3 (PBE)
HoF 3 (PBE+U d /3 eV/4f-in-Core) To the best knowledge of the authors, no first principle surface stability analysis of any compound within the whole structure type is available in the literature. The only surface calculation of any compound of β-YF 3 structure was published in 2013 by Ye et al. [39] on two selected surfaces of DyF 3 (001) and (101). They only calculated the surfaces matching their experimentally obtained nano-plates. However, these and other experiments on this class of compounds clearly demonstrate that the obtained surface structures are very dependent on the experimental conditions, especially on the utilized nature and geometry of the substrate, as well as on the solvent and fluoride concentration [14,16,[40][41][42][43][44][45]. The present work analyzes the inherent quantum chemical stability of all of the seven low Miller indices (hkl) surfaces, namely (001), (010), (100), (011), (101), (110) and (111). A previous study on another metal trifluoride, AlF 3 revealed stoichiometric or substoichiometric surfaces with a small fluorine-deficit as the most stable terminations [46]. Additionally, a substoichiometric fluorine content has also been found for YbF 3 thin films made from ion assisted deposition [47]. Consequently, this study includes all possible stoichiometric terminations and those with a small-to-moderate fluorine-deficit of 1-2 fluorine atoms per surface unit cell. This results in a scope of 24 terminations. The obtained surface energy results are combined with the geometry of the surface cut by a Wulff analysis to examine the expected surface abundance [48,49].

Computational Details
All calculations were performed in the Vienna Ab Initio Simulation Package (VASP, version 5.4.4) [50] on the supercomputer cluster HLRN in Berlin and Göttingen, Germany. using periodic density functional theory (DFT) with a generalised gradient approximation (GGA). As an exchange-correlation functional, the one of Perdew-Burke-Ernzerhof (PBE) is applied [51]. The inner shell electrons were described by the projector augmented wave (PAW) method [52,53]. The outer shell electrons were expanded in plane waves.
For converged YF 3 total bulk energies, the VASP potential files F_h ("hard", 7 electrons) and Y_sv (11 electrons) were applied together with a 9 × 9 × 9 Monkhorst-Pack grid. In accordance with the F_h potential file, a kinetic energy cut-off of 772.6 eV was used. For HoF 3 , both available Ho potential files Ho_3 (9 electrons, 4f-in-core) and Ho (21 electrons, 4fin-valence) were evaluated with respective grid sizes of 7 × 7 × 7 and 3 × 3 × 3. On holmium, the Hubbard-type correction in the simple Dudarev formalism was applied [54]. In a test series of 1-10 eV in 1 eV steps with U d (with Ho_3) and U f (with Ho), PBE+U d with 3 eV agreed best with the crystal structure and the presumed electronic structure (Table 1 and Figures S1 and S2). As an electronic structure reference, bulk HoF 3 was also calculated with the Heyd-Scuseria-Ernzerhof hybrid functional (HSE06) [55].
For electron smearing, tests on several bulk and slabs structures of both trifluorides were performed, comparing Gaussian smearing with the tetrahedron method with Blöchl correction [56]. No energy difference within the applied self-consistent field (SCF) convergence criteria could be found.
We therefore used Gaussian smearing on our insulating trifluorides.
Apart from the trifluorides, molecular fluorine, as well as metallic yttrium and holmium, were also considered. The first was calculated in a cubic box of 25 Å length. For the latter two, Gaussian smearing could not be applied. A convergence test with 1st and 2nd-order Methfessel-Paxton smearing with widths of 0.05-0.35 eV yielded 2ndorder Methfessel-Paxton smearing with a width of 0.10 eV (Y) or 0.15 eV (Ho) as the best combination to minimize the difference between total energy and free energy.
Each bulk structure started from the respective, experimental crystal structure (YF 3 [23], HoF 3 [22], Y [57], Ho [58]) and was fully relaxed in atomic positions, lattice constants and volume. The accurate precision setting was applied. As convergence criteria, 0.01 meV per unit cell was used for SCF total energies and 0.1 meV per unit cell for the difference in total energy between two ionic steps. Final total energies, density of states (DOS) and Bader charges were performed with an SCF criteria of 0.001 meV. All DOS plots and Bader charges, as well as all HoF 3 data, were calculated with allowed spin polarization. To aid SCF convergence, an additional support grid (.ADDGRID.) and/or a reduced minimal mixing parameter for Kerker's initial approximation [59] (AMIN) of <0.01 were applied on most slabs.
Symmetric slabs were built from the relaxed bulk structure with the Python package pymatgen [60,61]. The vacuum height perpendicular to the surface was tested for one stoichiometric termination of YF 3 (001). The converged value of 25 Å was applied for all slabs. For slab calculations, only one k-point was used perpendicular to the surface. For the other two directions, we applied the same k-point grid size as in bulk. The complete slabs were relaxed in atomic positions.
DOS plots and band structures were generated with pymatgen. Wulff plots were constructed with the WulffPack Python package [62]. Atomic structures were visualized in VESTA [63].

Choice of Electronic Structure Method
The effect of dispersion was tested by applying Grimme's dispersion correction with Becke-Johnson damping (D3(BJ)) [64]. From PBE to PBE+D3(BJ), the lattice constants changed only by 1.9-4.5 pm or 0.3-1.0% during the full optimization of atomic positions, lattice constants and volume of YF 3 . Due to this small deviation, we neglected dispersion correction for our highly ionic systems.
For HoF 3 , a test series was performed to decide whether to treat the 4f-electrons inside the core or at the valence level. Hubbard-type Coulomb parameters of 1-10 eV were scanned for the 4f-in-core with U d acting on Ho-d orbitals, as well as for 4f-in-valence with U f acting on Ho-f orbitals. It should be noted that the Ho-5d orbitals mainly constituted the broad conduction band in both approaches. Yet, they also hybridized in the valence band mainly constructed by F-2p ( Figure S3). The PBE+U benchmark plots for unit cell parameters and band gaps are given in the SI with further discussion (Figures S1 and S2). All HoF 3 (PBE+U d /4f-in-core) band structures resembled the YF 3 (PBE) one and produced comparable F-2p to Ho-5d or Y-4d charge transfer band gaps of 7-8 eV ( Figure S3). By adding exact exchange via HSE06/4f-in-core, these bands were further separated to 11 eV. Whereas, HSE06/4f-in-valence predicted an Ho-4f to Ho-4f transition of 8 eV. In contrast, PBE/4f-in-valence was not able to separate the partially filled 4f 10 into un-/occupied bands. Instead, it placed the Fermi-level (E F ) inside the 4f band, predicting a pseudo-metal. When introducing the additional Coulomb potential of 1-10 eV onto the 4f in PBE+U f , this 4f-4f gap was tuneable from 1 eV to a maximum of 6 eV. At U f ≥ 5 eV, the nature of the band gap changed to a charge transfer of F-2p to Ho-4f. Unfortunately, no measured band gap exists in the literature for HoF 3 . Therefore, it was not possible to pin-point the true band gap, nor to evaluate the correct nature of that transition. Nevertheless, based on a purely empirical model derived from other lanthanide compounds, HoF 3 is expected to have a band gap of ca. 9 eV [65]. This empirically estimated band gap, as well as the calculated HSE06 reference, were best reproduced without including the 4f-electrons explicitly.
Another quantity upon which to judge the applied electronic structure method was the Bader charges obtained by applying the atoms in molecules (AIM) population analysis [66][67][68][69][70]. For both bulk materials of YF 3 and HoF 3 , all tested methods predicted a metal charge of 2.2-2.4 e and fluorine charge of −(0.7-0.8) e. The Bader charges of all applied methods with 4f-in-core or valence agreed well with each other and thus suggested that including 4f explicitly was not necessary for HoF 3 .
Furthermore, all methods used with 4f-in-valence predicted a high-spin bulk unit cell with all four holmium aligned resulting in an electronic magnetic moment of 16 µ B . This ferromagnetic result was obtained even when starting from anti-ferromagnetic spin arrangements. According to the experimentally known magnetic structures, the physically correct spin arrangement is anti-ferromagnetic below 0.53 K or paramagnetic above [71].
To summarize, not including the 4f-electrons explicitly provided the best electronic structure results. The differences between simple PBE and PBE+U d were minor. When considering the unit cell parameters given in Table 1, PBE+U d /3 eV/4f-in-core performed best with deviations of as little as 0.1-0.8%.

Surface Energy
The surface formation energy (E surf ) is generally calculated from the total energy of the 2D-periodic slab (E n ), the energy of the 3D-periodic bulk unit cell (E bulk ) and the surface area of the slab (A): n is the slab thickness measured in unit cells. We label this bulk-derived surface energy E bd surf . Equation (1) is used for all YF 3 surface energies. In this work, we also considered surfaces with a substoichiometric amount of fluorine. For these, the fluorine potential µ F for each missing fluorine was added to the numerator of Equation (1). µ F was obtained from E bulk , the bulk energy per atom of the pure metal of yttrium or holmium (µ M ), as well as the number of metal (n M = 4) and fluorine (n F = 12) atoms within the bulk MF 3 unit cell: Yet, as pointed out by Boettger, this bulk-derived surface energy (E bd surf ) can lead to diverging E surf with respect to n [72]. This can be avoided by using slab-derived (sd) energies only: E bulk is then replaced by the difference of E n to the total energy of the next smaller slab (E n−1 ). For HoF 3 , we indeed observed linearly diverging E bd surf when applying Equation (1), despite system sizes of up to 7 UC or Ho 28 F 84 . Depending on the (hkl), this stoichiometry corresponds to 12, 24 or 26 HoF 3 -layers. Likely, this is a result of the allowed spinpolarization with Hubbard-type correction and atomic relaxation of the whole slab. It can be seen in Table S3, that this linear divergence only appears after relaxation in E bd surf,opt . The unrelaxed surface energies E bd surf,SP show no divergence. In YF 3 , no divergent E bd surf,opt are observed. Here, no Hubbard-type correction is applied and the atomic relaxation is performed without spin polarization. A comparison of slab convergence by both equations is given in Tables S2 and S3. Due to the divergence issue, all HoF 3 surface energies given within this paper are slab-derived using Equation (3), which nicely converge. As each E sd surf is derived from two slabs differing by one unit cell in size, at least three slabs are needed to determine convergence. Whereas for E bd surf , these are just two. Due to the observed convergence of E bd surf in YF 3 , only two slab thicknesses are modeled for many terminations. Therefore, the convergence of the respective E sd surf cannot be evaluated. As a consequence, we used the converged E bd surf for YF 3 to compare with the converged HoF 3 E sd surf . All YF 3 bulk-derived surface energies converged within 0.03 J m −2 at slab thickness of about 5-5.5 UC or 10-22 YF 3 -layers (Table S2). The HoF 3 slab-derived surface energies of 14 terminations, including all of the most stable ones per Miller indices, converged to 0.01 J m −2 or less within a slab thickness of about 6-6.5 UC or 12-26 HoF 3 -layers (Table S3). Some of the higher energy terminations converged only to 0.02-0.04 J m −2 at that thickness, while four high energy terminations did not converge even to 0.1 J m −2 . Fortunately, it is clear from their surface energies that even within the present uncertainty, those high energy terminations do not compete with the lowest energy ones. The slab thickness convergence for HoF 3 is visualized by error bars in Figure S5.

Bulk Properties
For YF 3 , the PBE relaxed lattice constants, given in Table 1, agree very well with both experimental values, which are underestimated by as little as 0.5-1.5% [22,23]. The resulting unit cell volume is underestimated by 2.6-2.7%, which is still in good agreement for a GGA functional. The best performing HoF 3 method against the only available experimental unit cell data and the calculated HSE06 band gaps was found to be PBE+U d /3 eV/4f-in-core. The resulting unit cell parameters deviate by as little as 0. Before we come to the surfaces, we evaluate possible energetic differences between the two geochemical twins as bulk materials. We calculate the electronic contribution to the formation enthalpies (∆H f ) according to: The electronic energies are taken from the bulk metals in hcp (P63/mmc) structure and the bulk trifluorides, as well as molecular fluorine. For YF 3 , we obtained an electronic contribution of −1591.1 kJ mol −1 versus −1587.3 kJ mol −1 for HoF 3 . Thus, judged by the electronic energies only, both trifluorides are equally strong bound with a very small favor of −3.7 kJ mol −1 or 0.2% for YF 3 over HoF 3 .

Surface Energies
The surface energies of all calculated terminations are given in Table 2. The given metal surface coordination number (CN surf ) is determined with a bond length cut-off of R F-M ≤ 2.60 Å. Table 2 also includes the nominal net surface charge (q surf ) caused by substoichiometric fluoride. Finally, the last column includes the surface abundance for each respective most stable termination predicted by Wulff construction (% surf ). The two terminations, (110)-1 and -2 greatly illustrate the importance of atomic relaxation of the surface prior analysis. Before relaxation, nothing but the very surface layer differs within each (hkl) cut. As both terminations are stoichiometric, they are also identical in composition. However, for both trifluorides, the unrelaxed (110)-1 surface is by 0.6 J m −2 more stable than the one of (110)-2 (see E unrel. surf in Table 2). When allowed to relax in atomic positions, the {5,9,8} surface coordinations of (110)-2 rearrange into {6,8,8} (Table S1 and Figure S4). Hence, the surface energy reduces by as much as 1.41 J m −2 for YF 3 or 1.18 J m −2 for HoF 3 . In contrast, termination (110)-1 already starts at a higher surface coordination of {6,9,8}, before it also rearranges into {6,8,8}. According to the lesser degree of rearrangement, its surface energy only reduces by 0.79 J m −2 for YF 3 or 0.60 J m −2 for HoF 3 . Both rearranged terminations are structurally equivalent.
The argumentation in CN surf cannot only be applied to explain the high E unrel. surf of some terminations, but is also partially applicable to the relaxed E surf . Within all YF 3 (hkl) subsets, except those of (101) and (111), the respective smallest CN surf value correlates with E surf . Thus, the smaller the smallest coordination polyhedron, the less stable the surface and the higher its surface energy. For example, a 4-fold coordination, as present in many surfaces with a fluorine-deficit of two, is only found for the highest E surf within the (hkl) subset. Yet, this correlation holds only for the very minimal value within a set of CN surf . No correlation can be found for the remaining, higher CN surf values of the same surface. Therefore, these cannot explain the energetic order of two terminations showing the same smallest CN surf value (as seen e.g., in YF 3 (100)-1 and -3 in Table 2). For HoF 3 , this correlation of surface coordination and stability has two more exceptions. Here, the stability of the least and second-least stable (100) and (001) terminations flip compared to YF 3 , without any change in CN surf . Prior to surface relaxation, all coordination polyhedrons of YF 3 and HoF 3 are identical as they share the same crystal bulk structure. After relaxation, this is still true for twenty terminations (Table S1). Only four rearranged terminations differ slightly in surface coordination between YF 3 and HoF 3 . All of these four terminations belong to the less stable surfaces within the respective (hkl). All most or second-most stable terminations are identical in surface coordination between YF 3 and HoF 3 . The most stable surface structure termination for each of the seven Miller indices is shown in Figure 2. As shown in Figure 3, the obtained E surf are similar in magnitude and, within most Miller indices, the order of terminations is equal between YF 3 and HoF 3 . Within convergence, this is also true for the two stoichiometric terminations of (110) and (101), which are very similar in surface energy. For (100) and (001), the least and second-least stable terminations switch their order between YF 3 and HoF 3 . Here, HoF 3 prefers the surface with a nominal surface net charge of +2 over the stoichiometric one. The only difference in termination order between the two compounds, which also affects the most stable surface, is found in (111). For YF 3 , the most stable surface is (111)-2, which shows a surface coordination of CN surf = {6, 5, 8, 8}. Whereas, HoF 3 prefers (111)-3 with CN surf = {6, 7, 6, 9} ( Figure 2). However, both of these terminations are equal in constitution with a fluorinedeficit of 1 per surface. Even though the order within one (hkl) is largely the same between YF 3 and HoF 3 , the order between the different (hkl) does change. For YF 3 , the overall two most stable surfaces are (010)-1 and (001)-2, which are both stoichiometric and give a surface energy of 0.58 J m −2 . This is closely followed by the stoichiometric surface (011)-2. Medium stable surfaces are found for (101)-3 and (111)-2, which are both substoichiometric surfaces missing a single fluorine. The two least stable surfaces (110)-1/-2 and (100)-2 prefer a stoichiometric termination again.
(010)-1 is also the overall most stable surface for HoF 3 , but the remaining surfaces differ in order. (100)-2, which is the most unstable (hkl) in YF 3 , is the second-most stable one in HoF 3   and HoF 3 (right, PBE+U d /3 eV/4f-in-core) from relaxed surfaces. The percentage shows the relative abundance of each surface, which is also given in Table 2.
From the corresponding energies of the respective most stable surfaces shown in Figure 2, the Wulff plots are constructed. A Wulff plot visualizes the thermodynamically most stable crystal shape at quantum chemical conditions of 0 K and vacuum. To test the dependence of surface ratios on the slab thickness convergence, an estimation on the maximum possible error is given in the SI (Table S4).
The largest surface area of over one quarter in YF 3 or one third in HoF 3 is formed by (010). That HoF 3 prefers (010) even stronger is a consequence of its surface energy being 0.11 J m −2 more stable than any other. This contrasts YF 3 , for which (001) has the same surface energy as (010), as well as a very closely (0.03 J m −2 ) following (011). Nonetheless, the geometric interdependence of surfaces cause a much smaller abundance of only 10% for (001) versus more than double for (011) in YF 3 . The third most abundant surface in YF 3 is (101) with 20%, which is one of the two obtained substoichiometric surfaces. As these same three surfaces (001), (011) and (101) are only medium stable in HoF 3 , they also only constitute 6%, 13% and 14% of the overall surface. The second substoichiometric surface present in both Wulff plots is (111), which forms an area of 10% in YF 3 and 7% in HoF 3 . Thus, almost one third of the YF 3 crystal is made from terminations with a nominal positive net charge of +1. Whereas for HoF 3 , these are only about a fifth. The two least stable surfaces of YF 3 are (100) and (110), which constitute about 7% and 5%. In HoF 3 , the latter is to such an extend energetically unstable, that it is completely excluded from the Wulff plot. (100), on the other hand, turns to be the second most stable and second most abundant surface in HoF 3 . It constructs one quarter of the total surface.
The overall scale of most stable surface energies per Miller indices is comparable between YF 3 and HoF 3 . The respective ranges are 0.58-1.03 J m −2 for YF 3 and 0.47-1.00 J m −2 for HoF 3 . However, the resulting average by the Wulff plot is 16% higher for YF 3 with ∅E surf = 0.70 J m −2 than for HoF 3 with ∅E surf = 0.59 J m −2 . This means, that forming surfaces from the bulk crystal involves a higher thermodynamic barrier in YF 3 than in HoF 3 . This is interesting, as there is no significant difference within the formation enthalpies of the bulk. The difference in average surface energy also hints, that the thermodynamic barrier of crystal nucleation is also higher in YF 3 than in HoF 3 . Though, to accurately predict the nucleation, the nature of the respective precursors and the media needs to be considered.

Surface Band Gaps
Investigating the electronic properties, we found that the band gap, total DOS and projected DOS of atoms central within the slab converged already at the smallest slab thickness. All band gaps are plotted in Figure S6. For the most stable termination of each (hkl), the DOS near the Fermi level are also shown in Figure S7. The surface DOS narrow the band gap within all terminations. For stoichiometric slabs, the direct band gap is reduced from the bulk value of ca. 8 eV to 4-7 eV. Thus, all stoichiometric surfaces remain fully within the insulating regime. In contrast, for substoichiometric terminations, the direct band gap collapses to 0-1 eV, predicting a pseudo-metallic or narrow-band-gap surface. For HoF 3 (101) and (111), only one spin direction shows a nearly metallic character. The other stays insulating (5-6 eV). It should be noted that this pseudo-metallic or narrow-band-gap electronic structure at the substoichiometric surfaces might be strongly effected by the chosen neutral 2D-periodic model. However, as this paper is focusing on the relative stability of surfaces, we are not investigating the nature of the band gaps, nor the observed surface magnetism or spin-asymmetric band gaps of some HoF 3 surfaces further.

Conclusions
The aim of this study was to obtain the relative surface stabilities, in order to find the most abundant surfaces of the two REE trifluorides, YF 3 and HoF 3 according to their inherent quantum chemical stability. While YF 3 can be treated on the DFT level, the 4felectrons of HoF 3 required an extensive electronic structure benchmark evaluating DFT, DFT+U and hybrid DFT against the crystal unit cell, band gaps and Bader charges. On the DFT or DFT+U level, our results show that including the 4f-electrons explicitly within the plane wave expansion worsens the geometrical and band gap results, while the Bader charges stay unaffected. Considering also the experimentally not observed high-spin preference of the 4f-electrons, as well as the strongly increased computational demand, we treated HoF 3 by a 4f-in-core DFT+U approach, in which the Hubbard-type correction is applied on the Ho-5d orbitals, which mix into the valence band mainly constructed by F-2p.
From the relaxed bulk, surface models were created for any of the seven low-lying Miller indices. Our analysis included all possible stoichiometric terminations, as well as those showing a small to moderate fluorine-deficit. The surfaces were quantified by Bader charges, band gaps and DOS. From the resulting scope of 24 surfaces, we constructed the first Wulff plots for the whole class of β-YF 3 -structured compounds.
We found that, within each Miller indices, both trifluorides prefer the same termination with the exception of (111), in which different surface coordinations are favored.
Comparing the different Miller indices, both compounds clearly show stoichiometric (010) as the most stable surface. The preference of the other surfaces, though, varies between the two. The greatest difference is found for (100), which is the second-most stable surface for HoF 3 , but the second-least stable one for YF 3 . On average, the surface energy predicted by the Wulff plot is higher for YF 3 than for HoF 3 . This suggests a higher thermodynamical barrier for the formation of YF 3 surfaces from the bulk.
In total, one third of the predicted equilibrium crystal shape of YF 3 is made from the substoichiometric terminations of (101) and (111) missing a single fluorine per surface. In HoF 3 , these only constitute a fifth. In the search for the underlying reason between the different fluorine affinity of the two compounds, this different availability of substoichiometric surfaces is an interesting finding. However, to evaluate possible effects, further studies are needed that actually model binding interactions with these surfaces. These should also apply more elaborate binding analysis tools than simple population analysis.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ma15176048/s1, Figure S1: Calculated HoF 3 direct band gaps with HSE06 (blue), PBE+U (green) and pure PBE (red) applied on 4f-in-core (full markers) or 4f-invalence (crosses); HSE06/4f-in-valence is not relaxed but done on-top of the crystal structure; the area between the two HSE06 values is highlighted in blue; Figure S2: Calculated HoF 3 relaxed unit cell parameters with HSE06 (blue), PBE+U (green) and pure PBE (red) applied on 4f-in-core (full markers) or 4f-in-valence (crosses) compared to the experimental values (horizontal line); Figure S3: Bulk band structure, total DOS (tDOS: gray) and DOS projected onto the metal d band (blue) or fluorine 2p band (green): (a) YF 3 (PBE) and (b) HoF 3 (PBE+U d /3 eV/4f-in-core); Table S1: Comparison of unrelaxed versus relaxed (or rearranged) slabs in metal coordination number at the surface (CN surf ), as well as in metal centers of the non-surface layers (CN non-surf ) as determined with a bond distance cut-off of 2.6 Å; Figure S4: Effect of surface rearrangement on the stoichiometric surface terminations of (110)-1 (left), (110)-2 (middle) and (101)-2 (right). Atomic positions are shown before (gray) and after relaxation (M: blue, F: green). For the latter, all polyhedra are shown but the one from the initially lowest surface coordination number (CN unrel surf ). Given are the surface energies in J m −2 of the unrelaxed surfaces (E unrel surf ) for YF 3 (first) and HoF 3 (second); Table S2: YF 3 (PBE) bulk-derived (E bd surf ) and slab-derived (E sd surf ) surface energies without (SP) and with atomic position relaxation (OPT); all energies in J m −2 ; the E bd surf,opt values are used within the main paper; Table S3: HoF 3 (PBE+U d /3 eV/4f-in-core) bulk-derived (E bd surf ) and slab-derived (E sd surf ) surface energies without (SP) and with atomic position relaxation (OPT); all energies in J m −2 ; all magnetic moments in µ B ; the E sd surf,opt values are used within the main paper; Figure S5: Relaxed slab-derived surface energies of HoF 3 (PBE+U d /3 eV/4f-in-core). The uncertainty due to slab thickness convergence is given by error bars on each termination; Table S4: Effect of maximal error accumulation due to the convergence in slab thickness of maximal ±0.03 J m −2 for YF 3 and ±0.01 J m −2 for HoF 3 onto Wulff construction; i denotes the initial value of average surface energy (∅E surf ) or surface abundance (% surf ) given by the Wulff plots in the main paper Figure 4; Figure S6: YF 3 (left, PBE) and HoF 3 (right, PBE+U d /3 eV/4fin-core) band gaps of surfaces compared with the respective bulk value (gray). Minimal band gaps, direct or indirect are given by solid bars. In the case, the minimal band gap was found to be indirect, also the direct band gap is given by a transparent bar. For HoF 3 (101) and (111)