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+Ud/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.

1. 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 wide-band-gap material YF3 has very good properties for laser applications [5,6,7,8]. Doped with trivalent REE cations, YF3 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]. HoF3 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, YF3 and HoF3 are important precursors for the synthesis of the respective pure metallic compounds [17,18]. In nature, YF3 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 YF3 and HoF3 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.

In accordance with their twin behavior, solid YF3 and HoF3 occur in the same crystal structure type of -YF3 (space group Pnma, fully occupied Wyckoff positions (Y) 4c, (F) 4c, (F) 8d (Figure 1)) [18,22,23]. This is their only stable phase up to 1343 K (HoF3) or 1350 K (YF3), which is well beyond the temperature regime of hydrothermal fluids of typically 323–873 K [7,18,24,25]. The ionic radii of the middle and late lanthanides Sm(III)–Lu(III) differ by only ≤6 pm compared to Y(III) [1]. Consequently, all their trifluorides crystallize as well in -YF3 [18,22,23,26,27,28,29,30,31]. The same low-temperature phase is also found for the two actinides Bk(III) and Cf(III) [32,33,34]. Due to the often observed analogy of actinides to lanthanides, the same crystal structure is assumed for the, so far, experimentally unknown heavier actinide trifluorides of Es(III)–Lr(III) [35]. Interestingly, the known orthorhombic low-temperature phase of plutonium trihydride is also reasoned to be an exotic example of a -YF3 structure [36]. Outside the f-block, the -YF3-structure is experimentally known for bismuth trifluoride [23,37,38] and predicted as an accessible meta-stable phase for the trichlorides of Y(III) and Bi(III) [38].

Figure 1

The orthorhombic unit cell of -YF3 with Pnma symmetry. Lattice constants are given in Table 1. The distorted tricapped trigonal prisms formed by nine fluorides (green) around each yttrium (cyan) are visualized by transparent, cyan planes.

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 -YF3 structure was published in 2013 by Ye et al. [39] on two selected surfaces of DyF3 (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 () surfaces, namely (001), (010), (100), (011), (101), (110) and (111). A previous study on another metal trifluoride, AlF 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 YbF3 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].

2. Methodology

2.1. 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 YF3 total bulk energies, the VASP potential files F_h (“hard”, 7 electrons) and Y_sv (11 electrons) were applied together with a Monkhorst–Pack grid. In accordance with the F_h potential file, a kinetic energy cut-off of 772.6 eV was used. For HoF3, both available Ho potential files Ho_3 (9 electrons, 4f-in-core) and Ho (21 electrons, 4f-in-valence) were evaluated with respective grid sizes of and . 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 (with Ho_3) and U (with Ho), PBE+U 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 HoF3 was also calculated with the Heyd–Scuseria–Ernzerhof hybrid functional (HSE06) [55].

Table 1

Comparison of the relaxed unit cells to experiment including standard deviation in parentheses [22,23]. Given is the absolute difference (), as well as the deviation from experiment in percentage ().

	YF3 (PBE)				HoF3 (PBE+U_d/3 eV/4f-in-Core)
	a ( Å )	b (Å)	c (Å)	V (Å ^3 )	a (Å)	b (Å)	c (Å)	V (Å ^3 )
calc.	6.3215	6.8059	4.3300	186.29	6.4164	6.8796	4.3440	191.76
exp. [22]	6.353(3)	6.850(3)	4.393(3)	191.2	6.404(3)	6.875(3)	4.379(3)	192.8
Δ	0.032	0.044	0.063	4.9	0.012	0.005	0.035	1.0
Δ_%	0.5%	0.6%	1.4%	2.6%	0.2%	0.1%	0.8%	0.5%
exp. [23]	6.3537(7)	6.8545(7)	4.3953(5)	191.42
Δ	0.0322	0.0486	0.0653	5.13	—
Δ_%	0.5%	0.7%	1.5%	2.7%

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 2nd-order 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 (YF3 [23], HoF3 [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 HoF3 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 YF3 (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].

2.2. 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 YF3. Due to this small deviation, we neglected dispersion correction for our highly ionic systems.

For HoF3, 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 acting on Ho-d orbitals, as well as for 4f-in-valence with U 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 HoF3 (PBE+U/4f-in-core) band structures resembled the YF3 (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 into un-/occupied bands. Instead, it placed the Fermi-level () inside the 4f band, predicting a pseudo-metal. When introducing the additional Coulomb potential of 1–10 eV onto the 4f in PBE+U, this 4f–4f gap was tuneable from 1 eV to a maximum of 6 eV. At U 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 HoF3. 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, HoF3 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 YF3 and HoF3, 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 HoF3.

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 µ. 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 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 were minor. When considering the unit cell parameters given in Table 1, PBE+U/3 eV/4f-in-core performed best with deviations of as little as 0.1–0.8%.

2.3. Surface Energy

The surface formation energy () is generally calculated from the total energy of the 2D-periodic slab (), the energy of the 3D-periodic bulk unit cell () and the surface area of the slab (A):

n is the slab thickness measured in unit cells. We label this bulk-derived surface energy . Equation (1) is used for all YF3 surface energies. In this work, we also considered surfaces with a substoichiometric amount of fluorine. For these, the fluorine potential µ for each missing fluorine was added to the numerator of Equation (1). µ was obtained from , the bulk energy per atom of the pure metal of yttrium or holmium (µ), as well as the number of metal () and fluorine () atoms within the bulk MF3 unit cell:

Yet, as pointed out by Boettger, this bulk-derived surface energy () can lead to diverging with respect to n [72]. This can be avoided by using slab-derived (sd) energies only:

is then replaced by the difference of to the total energy of the next smaller slab (). For HoF3, we indeed observed linearly diverging when applying Equation (1), despite system sizes of up to 7 UC or Ho28F84. Depending on the (), this stoichiometry corresponds to 12, 24 or 26 HoF3-layers. Likely, this is a result of the allowed spin-polarization 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 . The unrelaxed surface energies show no divergence. In YF3, no divergent 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 HoF3 surface energies given within this paper are slab-derived using Equation (3), which nicely converge. As each is derived from two slabs differing by one unit cell in size, at least three slabs are needed to determine convergence. Whereas for , these are just two. Due to the observed convergence of in YF3, only two slab thicknesses are modeled for many terminations. Therefore, the convergence of the respective cannot be evaluated. As a consequence, we used the converged for YF3 to compare with the converged HoF3 . All YF3 bulk-derived surface energies converged within 0.03 J m at slab thickness of about 5–5.5 UC or 10–22 YF3-layers (Table S2). The HoF3 slab-derived surface energies of 14 terminations, including all of the most stable ones per Miller indices, converged to 0.01 J m or less within a slab thickness of about 6–6.5 UC or 12–26 HoF3-layers (Table S3). Some of the higher energy terminations converged only to 0.02–0.04 J m at that thickness, while four high energy terminations did not converge even to 0.1 J m. 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 HoF3 is visualized by error bars in Figure S5.

3. Results and Discussion

3.1. Bulk Properties

For YF3, 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 HoF3 method against the only available experimental unit cell data and the calculated HSE06 band gaps was found to be PBE+U/3 eV/4f-in-core. The resulting unit cell parameters deviate by as little as 0.1–0.8%.

The respective F–M bond length on the PBE (YF3) and PBE+U level (HoF3) are Å and Å. These agree perfectly with the measured interatomic distances of 2.3–2.6 Å [22,23].

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 () 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 YF3, we obtained an electronic contribution of kJ mol versus kJ mol for HoF3. Thus, judged by the electronic energies only, both trifluorides are equally strong bound with a very small favor of kJ mol or 0.2% for YF3 over HoF3.

3.2. Surface Energies

The surface energies of all calculated terminations are given in Table 2. The given metal surface coordination number (CN) is determined with a bond length cut-off of 2.60 Å. Table 2 also includes the nominal net surface charge () caused by substoichiometric fluoride. Finally, the last column includes the surface abundance for each respective most stable termination predicted by Wulff construction (%).

Table 2

The YF3 (PBE) and HoF3 (PBE+U/3 eV/4f-in-core) surfaces with respective terminations (term.), slab thickness in layers of formula units without terminal F-deficit (), nominal surface net charge () in e, surface energies of relaxed () and unrelaxed slabs () in J m, as well as the relaxed surface metal coordination number (CN). The lowest surface energies per () cut are highlighted in bold. For these, also the abundance obtained by the Wulff plot (%) is given.

			L_MF_3		CN_surf		E_surf (Esurfunrel.)		%_surf
( hkl )	term.	q_surf	YF3	HoF3	YF3	HoF3	YF3           HoF3		YF3	HoF3
(100)	1	0	20	24	5, 9		1.61 (2.87)	0.93 (1.48)
	2	0	22	26	6, 9		1.03 (2.02)	0.58 (0.96)	7%	25%
	3	+1	20	24	5,8		1.24 (1.61)	0.62 (0.68)
	4	+2	22	26	4, 7		1.79 (2.14)	0.87 (0.90)
(010)	1	0	10	12	8,8		0.58 (0.84)	0.47 (0.49)	26%	34%
	2	+2	10	12	6, 6		1.80 (2.05)	1.52 (1.52)
(001)	1	0	20	24	5, 8, 8, 9		1.23 (2.45)	1.37 (2.25)
	2	0	22	26	6, 7, 8, 9		0.58 (1.39)	0.67 (1.16)	10%	6%
	3	+2	22	26	4, 5, 8, 9		1.27 (1.70)	1.23 (1.29)
(110)	1	0	20	24	6, 8, 8		1.01 (1.80)	0.99 (1.59)	5%	0%
	2	0	22	26	6, 8, 8		1.00 (2.41)	1.00 (2.18)
	3	+2	22	26	4, 6, 9	4, 6, 8	1.42 (1.73)	2.09 (1.36)
(101)	1	0	20	24	6, 7, 8, 8		0.82 (1.48)	0.89 (1.33)
	2	0	20	24	6, 6, 8, 8		0.82 (3.34)	0.88 (3.17)
	3	+1	20	24	6, 7, 8, 8		0.76 (1.16)	0.69 (0.89)	20%	14%
	4	+1	22	26	5, 6, 7, 9	5, 6, 8, 8	1.07 (2.10)	1.03 (1.70)
	5	+2	20	24	4, 5, 8, 8	5, 6, 8, 8	0.98 (1.39)	0.99 (0.99)
(011)	1	0	10	12	6, 6, 8, 8		0.78 (1.30)	0.81 (1.14)
	2	0	10	12	7, 7, 9, 9		0.61 (1.32)	0.68 (1.15)	22%	13%
	3	+2	10	12	4, 4, 8, 8		1.25 (1.68)	1.35 (1.38)
(111)	1	0	20	24	6, 7, 7, 8	7, 7, 8, 8	1.02 (3.46)	0.87 (3.29)
	2	+1	20	24	5, 6, 8, 8		0.83 (1.30)	0.82 (1.04)	10%
	3	+1	22	26	6, 6, 7, 9		1.05 (1.70)	0.75 (1.11)		7%
	4	+2	20	24	5, 5, 7, 7		0.93 (1.22)	0.95 (1.13)

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 () 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 more stable than the one of (110)-2 (see 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 for YF3 or 1.18 J m for HoF3. 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 for YF3 or 0.60 J m for HoF3. Both rearranged terminations are structurally equivalent.

The argumentation in CN cannot only be applied to explain the high of some terminations, but is also partially applicable to the relaxed . Within all YF3 () subsets, except those of (101) and (111), the respective smallest CN value correlates with . 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 within the () subset. Yet, this correlation holds only for the very minimal value within a set of CN. No correlation can be found for the remaining, higher CN values of the same surface. Therefore, these cannot explain the energetic order of two terminations showing the same smallest CN value (as seen e.g., in YF3 (100)-1 and -3 in Table 2). For HoF3, 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 YF3, without any change in CN. Prior to surface relaxation, all coordination polyhedrons of YF3 and HoF3 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 YF3 and HoF3. All of these four terminations belong to the less stable surfaces within the respective (). All most or second-most stable terminations are identical in surface coordination between YF3 and HoF3. The most stable surface structure termination for each of the seven Miller indices is shown in Figure 2.

Figure 2

Relaxed, most stable surface structures: the coordination number of the surface metals (CN) and the surface energies in J m () are given. The first entry corresponds to YF3, the second to HoF3. The mean of both values corresponds to the given order from top left to bottom right. Each () slab is rotated in a way to show the surface coordination best. For (111), two surfaces are given, as (111)-2 is preferred by YF3 and (111)-3 by HoF3.

As shown in Figure 3, the obtained are similar in magnitude and, within most Miller indices, the order of terminations is equal between YF3 and HoF3. 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 YF3 and HoF3. Here, HoF3 prefers the surface with a nominal surface net charge of 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 YF3, the most stable surface is (111)-2, which shows a surface coordination of CN. Whereas, HoF3 prefers (111)-3 with CN (Figure 2). However, both of these terminations are equal in constitution with a fluorine-deficit of 1 per surface.

Figure 3

Relaxed surface energies of YF3 (left, PBE) and HoF3 (right, PBE+U/3 eV/4f-in-core) of all 24 terminations. The surfaces are color-coded to their nominal surface net charge () in e of 0 (blue), (green) or (orange). This magnified plot does not show the HoF3 (110)-3 value of 2.09 J m.

Even though the order within one () is largely the same between YF3 and HoF3, the order between the different () does change. For YF3, 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. 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 HoF3, but the remaining surfaces differ in order. (100)-2, which is the most unstable () in YF3, is the second-most stable one in HoF3. The moderately stable surfaces (001)-2, (011)-2 and (101)-3 have equivalent surface energies within the slab thickness convergence of 0.01 J m. (111)-3 is the second least stable surface. The least stable surface is stoichiometric (110)-1/-2, which is not even part of the Wulff plot (Figure 4).

Figure 4

Wulff Plots of YF3 (left, PBE) and HoF3 (right, PBE+U/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 YF3 or one third in HoF3 is formed by (010). That HoF3 prefers (010) even stronger is a consequence of its surface energy being 0.11 J m more stable than any other. This contrasts YF3, for which (001) has the same surface energy as (010), as well as a very closely (0.03 J m) 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 YF3. The third most abundant surface in YF3 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 HoF3, 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 YF3 and 7% in HoF3. Thus, almost one third of the YF3 crystal is made from terminations with a nominal positive net charge of +1. Whereas for HoF3, these are only about a fifth. The two least stable surfaces of YF3 are (100) and (110), which constitute about 7% and 5%. In HoF3, 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 HoF3. It constructs one quarter of the total surface.

The overall scale of most stable surface energies per Miller indices is comparable between YF3 and HoF3. The respective ranges are 0.58–1.03 J m for YF3 and 0.47–1.00 J m for HoF3. However, the resulting average by the Wulff plot is 16% higher for YF3 with J m than for HoF3 with J m. This means, that forming surfaces from the bulk crystal involves a higher thermodynamic barrier in YF3 than in HoF3. 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 YF3 than in HoF3. Though, to accurately predict the nucleation, the nature of the respective precursors and the media needs to be considered.

3.2.1. Bader Charges

Figure 5 shows the partial charges obtained by Bader analysis for all 24 thickness-converged slabs of YF3 or HoF3. As these slabs are built from up to 104 atoms, a large number of very similar charges are obtained. To ease comparison, the Bader charges given in Figure 5 are rounded to 0.1 e. The Bader charges of the central slab atoms reproduce the bulk values with 2.4 e for yttrium, 2.3 e for holmium, as well as e for fluorine in both compounds. Only HoF3 (011)-3 shows a marginally increased central slab value of 2.4 e for holmium. In general, the Miller indices do not seem to affect the Bader charges within the analyzed accuracy. The highly ionic partial charge on fluorine does not significantly change for any slab with an overall range of −(0.7–1.0) e. Moreover, all stoichiometric surfaces have practically the same metal partial charges with 2.3–2.4 e. On the contrary, substoichiometric surfaces with a fluorine-deficit of one do all have at least one metal center charged less at the surface. These might be as low as 1.7 e, as found for the third termination of (100) in both trifluorides. Subsequently, surfaces with a fluorine-deficit of two contain even less ionic metal centers at the surface. The least charged metal center is observed again in (100) with only 1.2–1.3 e by the forth termination. The Bader charges suggest that all substoichiometric surfaces missing two fluorine have at at least a single surface metal center in an oxidation state of +II.

Figure 5

Bader charges rounded to 0.1 e for all slabs of YF3 (left, PBE) and HoF3 (right, PBE+U/3 eV/4f-in-core): the Bader charges of the bulk are highlighted by the dotted line (Y: 2.4 e, Ho: 2.3 e, F: e). Terminations are differently colored by their formal surface net charge () in e for 0 (blue cycle), (green diamond) and (orange cross).

3.2.2. 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 (), 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 HoF3 (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 HoF3 surfaces further.

4. 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, YF3 and HoF3 according to their inherent quantum chemical stability. While YF3 can be treated on the DFT level, the 4f-electrons of HoF3 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 HoF3 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 -YF3-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 HoF3, but the second-least stable one for YF3. On average, the surface energy predicted by the Wulff plot is higher for YF3 than for HoF3. This suggests a higher thermodynamical barrier for the formation of YF3 surfaces from the bulk.

In total, one third of the predicted equilibrium crystal shape of YF3 is made from the substoichiometric terminations of (101) and (111) missing a single fluorine per surface. In HoF3, 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.