Effective chemical potential for non-equilibrium systems and its application to molecular beam epitaxy of Bi₂Se₃

Abstract

First-principles studies often rely on the assumption of equilibrium, which can be a poor approximation, e.g., for growth. Here, an effective chemical potential ([small mu, Greek, macron]) method for non-equilibrium systems is developed. A salient feature of the theory is that it maintains the equilibrium limits as the correct limit. In application to molecular beam epitaxy, rate equations are solved for the concentrations of small clusters, which serve as feedstock for growth. We find that [small mu, Greek, macron] is determined by the most probable, rather than by the lowest-energy, cluster. In the case of Bi₂Se₃, [small mu, Greek, macron] is found to be highly supersaturated, leading to a high nucleus concentration in agreement with experiment.

---

The concept of chemical potential is one of the most fundamental quantities associated with thermodynamic equilibrium. It generalizes the effect of the environment on the energetics of a single type of particle. Through the use of atomic chemical potentials, we can compare the energies of systems consisting of different numbers of atoms. In defect physics, this plays a central role in the calculation of the formation energies of defects and allows for a statistical determination of defect concentrations which determine the electronic properties of semiconductors. Further, surface reconstruction and even the shape of nanocrystals are intimately tied with the availability of different chemical species, with the surface energy generally depending on the atomic chemical potentials of the various species. While kinetic Monte Carlo (kMC) or molecular dynamics (MD) can be used to investigate system evolution, when comparing the energies of systems containing different numbers of particles, calculation within the grand canonical ensemble requires knowledge of the chemical potentials of the atomic species. The most ubiquitous approximation used to determine the range of chemical potentials is that of so-called near equilibrium growth (NEG),^1,3 in which the stoichiometric sum of the atomic chemical potentials is assumed to be equal to that of the bulk crystal. This assumption, however, is an outstanding issue in the theoretical community as growth itself is a highly non-equilibrium process.

The chemical potential, μ, is a thermodynamic quantity, defined by , where G is the Gibbs free energy and n is the number of the n^th species.⁴ For equilibrium or near-equilibrium growth (NEG), there is a well-defined relationship between the μ′s of the species in the system.⁴ Consider, for example, a A_mB_n binary, we have

mμ_A + nμ_B = μ_AmBn. (1)

As a result of eqn (1), only one μ is an independent variable, i.e., one can use μ_B for this purpose and μ_A will be given by μA = (μA_mB_n − nμB)/m. Should μ for a species i exceeds its upper bound μ_i0, defined as the chemical potential of the corresponding elemental solid or gas, this species will precipitate,^1,5 thereby preventing any further increase of μ. In other words,

μ_i ≤ μ_i0 (i = A, B). (2)

Using the μ′s, the state of the system, whether it is a defect, a surface, or a crystal structure (denoted here as X) is determined by its formation energy with respect to bulk,

ΔH_f[X] = E_tot[X] − E_tot[bulk] − ∑n_iμ_i, (3)

where E_tot is the total energy.¹

This NEG theory has been applied to numerous materials behaviors with ample successes, ranging from surface reconstruction^5,6 and defect physics,^1,7 to nanostructures.^8,9 For systems in which a subsystem of concern can be isolated from the rest due, for example, to high reaction barriers, one can also use the NEG theory, provided that one can establish an approximate equilibrium within the subsystem.^10–12 However, when such a division is not obvious and the system is far away from equilibrium, the NEG theory can run into problems. For instance, recent interest in the physics of topological insulators (TIs) has spurred considerable activity to grow high-quality TIs such as Bi₂Se₃ with controlled properties.^13–18 Molecular beam epitaxy (MBE) is widely accepted as an effective method to grow high-quality films.^15–18 However, the defect density is still far from being satisfactory, due in part to the small size of the islands. Strictly speaking, MBE is not an equilibrium process;¹⁹ its outcome may rely heavily on the dynamics of the source atoms/molecules. On the other hand, MBE is relatively simple for its slow growth rate and the simplicity of the sources. This raises the question: can one incorporate kinetics into the chemical potential-based approach for non-equilibrium process?

Note that considerable theoretical efforts have been made in the past to study growth. Analytical rate equation approach, based on empirical kinetic parameters, has been employed to study island formation and growth.^20–26 Statistical approaches such as kinetic Monte-Carlo simulation and molecular dynamics have been used to study film growth.^27–32 Molecular dynamics based on first-principles calculations could provide a vast amount of information, however, it is generally limited to smaller systems and shorter times due to the high computational cost. kMC offers a great advantage for efficiency and capability by simulating the time evolution of the system via statistically evaluating the probability of selected events. kMC relies on knowing the rate and mechanism of all the relevant transitions from a given initial state, with loss of reality. Moreover, it appears that these non-equilibrium approaches do not use the concept of chemical potential, which had served as a basis for NEG, except perhaps for the thermodynamic theory of nucleation.² Recently, a renewed interest in the chemical potential has emerged, as it can help choosing lowest-energy paths for growth,³³ without knowing relevant transitions in advance during a significant time scale.

In this paper, we derive a generic effective chemical potential ([small mu, Greek, macron]) theory, which is suited for non-equilibrium physics. Taking the MBE growth of Bi₂Se₃ as an example, the use of [small mu, Greek, macron] allows us to simplify the overall complex process. We divide the non-equilibrium growth (nonEG) into three stages: pre-nucleation, nucleation, and island growth. For pre-nucleation, we determine [small mu, Greek, macron] by explicitly solving rate equations. First-principles calculation shows that [small mu, Greek, macron] is solely determined by the most probable, rather than by the lowest-energy, cluster(s) on the surface. For nucleation, it is found that the nucleation barrier vanishes when the critical size of the clusters is only a few molecule large due to the highly supersaturated nature of [small mu, Greek, macron]. This results in a high concentration of nuclei – a conclusion that contradicts the NEG model, but is in qualitative agreement with experiment.

One can define chemical potentials for an individual cluster α, irrespective of equilibrium,^4,34

where c_α is the areal concentration, G_α is the Gibbs free energy, E_0α is the energy of the cluster on the surface relative to the isolated constituent atoms, N_S is the number of sites per area that can hold the cluster (which, for simplicity, has been assumed to be much larger than the number of clusters), and k and T have the standard definitions. Here, considering a binary AB system (α = A_sB_t) as in NEG, we postulate that [small mu, Greek, macron] that enters eqn (3) is given bywhere the sums over s and t run over all possible AB clusters. As the concentrations of the clusters depend on the kinetics of the system, this weighted average incorporates both the energetics and kinetics in a simple way. When the system contains only a single element, eqn (5) is reduced towhere the subscript 0 denotes single-element quantity, as in eqn (2). For a binary system, as long as the formation of A_sB_t from A and B clusters is exothermic (i.e., with ΔE < 0), eqn (6) cannot hold except when species A (or B) starts to precipitate as pure clusters. Hence, it sets an upper bound for [small mu, Greek, macron]_A. The same is true for species B. Therefore,

[small mu, Greek, macron]_i ≤ [small mu, Greek, macron]_i0 (i = A, B). (7)

It is necessary to point out that eqn (5) and (7) are analogous to eqn (1) and (2), respectively.

To calculate [small mu, Greek, macron], we solve the following rate equations,²

k^des_α = ν exp(−φ^des_α/kT), (9)

k⁻ⁱ_α = ν exp(−φ⁻ⁱ_α/kT), (10)

k⁺ⁱ_α = 2π(D_α + D_i), (11)

D_α/i = a²ν exp(−φ^diff_α/i/kT), (12)

where F_α is the deposition rate, k^des_α is the desorption rate, k⁻ⁱ_α and k⁺ⁱ_α are the dissociation and association rates by emitting and absorbing a cluster i, respectively, D_α/i is the diffusion coefficient of cluster α and i³⁵ with a being the in-plane lattice constant, and ν ∼ 10¹³ s⁻¹ is the vibrational frequency. The use of eqn (11) above assumes that the process is diffusion limited [see discussion in ESI and Fig. S1† for the general case]. For simplicity, in the following we consider only a flat surface, so the effect of the Ehrlich–Schwoebel barrier² can be ignored. We also assume a steady-state, i.e., letting ċ_α≡0 in eqn (8).

Note that in the limit when all the barriers in eqn (8)–(12) can be ignored, the system is dominated by its lowest-energy form, namely, the concentrations for all clusters, as well as for all finite-sized islands diminish, leaving only the A_mB_n bulk. In this case, eqn (5) and (7) become eqn (1) and (2) respectively, i.e., the system approaches its NEG limit. Also note that chemical potential is not an easily measurable quantity. People may, instead, use the off-stoichiometry of the system, which can be directly measured, as an indirect indicator of μ_i. This is true not only for the NEG model³⁶ but also for the nonEG model here.

First-principles calculations were performed using the VASP code³⁷ to determine the relevant cluster size, structure, and energy barrier. The projector augmented wave (PAW) method³⁸ and the local density approximation (LDA) to the exchange–correlation functional³⁹ were used. Plane waves with a kinetic energy cutoff of 210 eV were used as the basis set. Integration of the Brillouin zone was done with sufficient k-point sampling⁴⁰ such that the numerical error is less than 0.01 eV. All atoms were fully relaxed until forces on the atoms were less than 0.01 eV Å⁻¹. The improved tangent finding method,⁴¹ within the framework of the nudged elastic band (NEB) method, was used to determine the energy barriers. Experimental growth temperature of 500 K (ref. 15–17) was assumed in the rate calculations.

Denoting t₀ as the time when nucleation takes place, the aforementioned 3 growth stages can be restated as follows: stage 1: t ≤ t₀ (pre-nucleation), stage 2: t ∼ t₀ (nucleation), stage 3: t ≫ t₀ (island growth). Fig. 1 illustrates schematically the evolution of the size of the clusters, and those of the nuclei and islands, and the corresponding μ′s.

Fig. 1 A schematic illustration of cluster/island distribution and the corresponding chemical potentials in three different stages of Bi₂Se₃ growth: (a) pre-nucleation, (b) nucleation, and (c) island growth. t₀ is the characteristic nucleation starting time.

In stage 1, there are only clusters, no nuclei, on the surface of Bi₂Se₃, whose thickness in the calculation is taken as five atomic layers, or called one quintuple layer (QL). Because the binding between QLs is van der Waals (vdW),^13–15 adding QLs to the thickness of the substrate is unlikely to affect our results. An in-plane supercell of 6 × 6 is used to conduct the cluster-adsorption calculation to avoid interactions among periodic images. Fig. 2 depicts the clusters we considered in our ab initio calculations. They are organized according to the chart in Fig. 2(a), where dashed line denotes the cutoff size of the clusters considered in this study (see details in Fig. S2, ESI†). Our calculations show that stable clusters on the surface are the same as those in vacuum. This can be expected as Bi tends to form 3 bonds, whereas Se tends to form 2 bonds, in line with the octet rule.⁴²Fig. 3 shows the calculated diffusion barrier (φ^diff_α), desorption barrier (φ^des_α), and dissociation barrier (φ⁻ⁱ_α) for the clusters: φ^diff_α is generally small, often a couple tenth of an eV, except for atomic Se and BiSe₂. φ^des_α is higher, especially for smaller clusters. φ⁻ⁱ_α is also generally higher and depends sensitively on the reaction pathway. Both φ^des_α and φ⁻ⁱ_α exhibit a large variation over a couple eV.

Fig. 2 (a) Small clusters for Bi₂Se₃ growth. They are organized in terms of size: horizontal = increasing number of Se atoms; vertical = increasing number of Bi atoms. Dashed zigzag line is discussed in the main text. (b)–(d) Calculated most-stable cluster structures on a Bi₂Se₃ (0001) substrate: (b) pure Se, (c) pure Bi, and (d) mixed Bi_sSe_t clusters. Red atoms are Se whereas green atoms are Bi.

Fig. 3 Calculated barriers for (a) diffusion (φ^diff_α), (b) desorption (φ^des_α), and (c) dissociation (φ⁻ⁱ_α) of the clusters in Fig. 2. Cluster names for both (a) and (b) are given in (b). For cluster dissociation in (c), a cluster labeled in the horizontal axis can have multi-reaction pathways, which are labeled inside or by the vertical bars.

Fig. 4(a and b) show c_α and the corresponding μ_α, which reveal that the highest concentrations of clusters are that of atomic Se and BiSe₂, respectively. At first glance, it may appear counterintuitive that some clusters with lower μ_α, e.g., Bi₂Se₂ and Bi₂Se₃, also have low c_α. This is because of the low-desorption barrier of Bi₂Se₃ in Fig. 3, which depletes not only Bi₂Se₃ but also Bi₂Se₂ (via an association with atomic Se, see Fig. S1†). In MBE, the Se partial pressure is much higher than that of Bi. As a result, the population of Se is determined almost entirely by pure Se clusters viaeqn (7),

[small mu, Greek, macron]_Se = [small mu, Greek, macron]_Se0 = μ_Se(atom). (13)

Fig. 4 (a) Logarithmic concentration of small clusters in stage 1. (b) The corresponding chemical potentials by eqn (4) and [small mu, Greek, macron]'s by eqn (13) and (14), respectively. (c) Logarithmic concentration of small clusters in stage 3. Solid squares are based on the nonEG model with c_island = 10⁹ cm⁻², whereas open squares are the prediction by the NEG model.

And for Bi, one has

These [small mu, Greek, macron] values, [small mu, Greek, macron]_Bi = −4.44 eV and [small mu, Greek, macron]_Se = −3.82 eV, are given in Fig. 4(b) as horizontal dashed lines. Our study reveals that it is the most-probable, not the lowest-energy, clusters that determine [small mu, Greek, macron]. This conclusion is in line with the basic principles of statistical physics, irrespective of the computational details, whereby it reinforces the notion that using c_α as the weighting factor for [small mu, Greek, macron] in eqn (5) is physically correct. It has been experimentally realized that the observed clusters were not the lowest energy ones in silicon cluster anions under non-equilibrium conditions.

Note that there is a large (3.08 eV) difference between the current model (2[small mu, Greek, macron]_Bi + 3[small mu, Greek, macron]_Se = −20.34 eV) and the NEG model (2μ_Bi + 3μ_Se = −23.42 eV). The fact that the former is significantly higher than the latter gives rise to supersaturation during the growth. In other words, eqn (7) sets a new set of boundaries for chemical potentials as the prime reason that account for the differences between nonEG and NEG. Importantly, this finding allows one to reexamine defects in Bi₂Se₃, as it suggests that native defects may be more easily formed than predicted by NEG theory,⁴⁴ highlighting experimental challenges in growing high-quality epitaxial films.⁴⁵

Before finishing up the discussion of stage 1, we would like to point out that our consideration of the availability of clusters is important, but only a first, step in improving upon NEG theory. Within the framework of the current development, higher order effects (such as the kinetics associated with the incorporation of available species into the growth front) can also be incorporated, which however will be deferred to future work. For stages 2 and 3, i.e., nucleation and island growth, due to continued deposition, we do not expect the concentrations of the small clusters to change considerably (as will be shown later).

In stage 2, owing to the formation of nuclei, bifurcation of μ will take place – one for the nuclei (μ_nuc) and one for the small clusters ([small mu, Greek, macron]), as depicted in Fig. 1(b). Bifurcation reflects the fact that once nuclei are formed, it is energetically favorable for them to grow further by consuming available smaller clusters on the surface, rather than consuming themselves. This should be contrasted to NEG where bifurcation of μ is strictly forbidden. Once [small mu, Greek, macron] for clusters are determined, one may approach the nucleation problem in different ways. For example, to study the temperature T dependence of nucleus density, one could apply eqn (4)–(12) to potential nuclei [which should be among the larger clusters in Fig. 2 and beyond (=embryo nuclei)] to determine μ_nuc and c_α,nuc, respectively. The non-negligible association barriers [see, e.g., Fig. S1(a)† for clusters] suggest that increasing T could help increasing the nucleus size and hence decreasing its density. However, the window for the increase can be limited, as although dissociation barriers are usually higher than association barriers [e.g., by comparing Fig. S1(a)† with Fig. 3(c)], the difference may not be that dramatic. As a result, further increasing T could also lead to the dissociation of already-formed nuclei, thereby increasing the nucleus density. Although rather qualitative, the conclusion here is already in line with experiments.^16,43

If, however, our interest is only on the critical size of the nuclei, not on its T dependence, explicitly invoking kinetic theory here can be an overkill. Instead, we only need calculate the nucleation barrier, which is defined as the maximum of the Gibbs free-energy,²

ΔG = ΔE_tot − (n_Se[small mu, Greek, macron]_Se − n_Bi[small mu, Greek, macron]_Bi). (15)

After reaching the critical size, adding additional atoms/clusters to the nuclei is energetically favored. Using [small mu, Greek, macron] in stage 1, we obtain Fig. 5, showing that ΔG is remarkably small and even negative for (Bi₂Se₃)_n clusters with n > 3. Usually, ΔG for cluster of these sizes is significantly positive due to its high surface energy. Due to high [small mu, Greek, macron] and the relatively low energy of the vdW surface, however, this is no longer the case here. Fig. S3† shows that except for n = 1 and 3, the desorption barriers are reasonably large so the nuclei in Fig. 5 do grow into islands. Consequently, the island concentration is expected to be high. In other words, the small nuclei size is consistent with measured high island concentration, c_isl ∼ 10⁹ to 10¹⁰ cm⁻².^16,17

Fig. 5 Free-energy barrier ΔG for nucleation on various clusters. Shaded area is where only stoichiometric molecular clusters (Bi₂Se₃)_n are considered in the calculation.

In stage 3, the basic physics should be the same as in stage 2. Due to the presence of the islands, however, one should recalculate c_α and [small mu, Greek, macron] by adding a term −k⁺ⁱ_islc_islc_i to eqn (8), where c_isl is the concentration of the islands and k⁺ⁱ_isl = 2π(D_i). Fig. 4(c) shows the results for a typical c_isl = 10⁹ cm⁻². More results on the c_isl-dependence can be found in Fig. S4, ESI.† As it turns out, neither c_α nor [small mu, Greek, macron] is changed significantly. As a comparison, Fig. 4(c) also shows c_α under the NEG condition. They are many orders of magnitude smaller.

In summary, an effective chemical potential approach for non-equilibrium growth is developed. We find that [small mu, Greek, macron] for an atomic species is determined primarily by the most probable cluster in which the atom resides during growth. Because “most probable” is a balance between energetics and kinetics, our findings thus set a new criterion for most relevant events during growth. Application to MBE growth of Bi₂Se₃ suggests that [small mu, Greek, macron] is highly supersaturated, resulting in an exceedingly small critical size of nuclei. While a high concentration of islands is in agreement with experiment, our results also reveal that to grow better-quality films requires finding ways to stabilize the most probable clusters, thereby lowering [small mu, Greek, macron]. Our formulation is general. It may be used to study structure, surface morphology, and shape of a nanoparticle, as well as defect and impurity, in non-equilibrium-grown solids. While kinetic theory has been around for a long time, most first-principles studies today are still based on the NEG model. Our chemical potential-based development here offers a natural vehicle to transform such bulk studies into the more experimentally relevant non-equilibrium regime.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank X. Liu and J. Wang for stimulating discussions. Chemical potential theory was supported by the US Department of Energy (DOE) under Grant No. DESC0002623. Bi₂Se₃ growth study was supported by NSF Award No. EFMA-1542798. Work in China was supported by the Ministry of Science and Technology of China (Grant No. 2011CB921901 and 2011CB606405), and the National Natural Science Foundation of China (Grant No. 11334006). Supercomputer time was provided by NERSC under the Grant No. DE-AC02-05CH11231 and the Center for Computational Innovations (CCI) at RPI.

References

1. S. B. Zhang and J. E. Northrup, Phys. Rev. Lett., 1991, 67, 2339.
2. I. V. Markov, Crystal Growth for Beginners: Fundamentals of Nucleation, Crystal Growth and Epitaxy, World Scientific, New Jersey, 2nd edn, 2003.
3. C. G. Van de Walle and J. Neugebauer, J. Appl. Phys., 2004, 95, 3851.
4. P. Attard, in Thermodynamics and Statistical Mechanics, Academic Press, London, 2002, p. 1.
5. G.-X. Qian, R. M. Martin and D. J. Chadi, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 38, 7649.
6. C. G. Van de Walle and J. Neugebauer, Phys. Rev. Lett., 2002, 88, 066103.
7. C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti and C. G. Van de Walle, Rev. Mod. Phys., 2014, 86, 253.
8. N. Moll, A. Kley, E. Pehlke and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 8844.
9. N. Wang, Y. Sun, Y. Zhang, D. West, W. Duan and S. Zhang, Phys. Rev. B, 2016, 93, 115306.
10. S. B. Zhang and S.-H. Wei, Phys. Rev. Lett., 2001, 86, 1789.
11. Y. Yan, S. B. Zhang and S. T. Pantelides, Phys. Rev. Lett., 2001, 86, 5723.
12. Z. Zheng, Y.-Y. Sun, W. Xie, J. Zhao and S. Zhang, J. Phys. Chem. C, 2016, 120, 27085.
13. W. Lu, Y. Ding, Y. Chen, Z. L. Wang and J. Fang, J. Am. Chem. Soc., 2005, 127, 10112.
14. D. Kong, W. Dang, J. J. Cha, H. Li, S. Meister, H. Peng, Z. Liu and Y. Cui, Nano Lett., 2010, 10, 2245.
15. X. Chen, X.-C. Ma, K. He, J.-F. Jia and Q.-K. Xue, Adv. Mater., 2011, 23, 1162.
16. L. He, X. Kou and K. L. Wang, Phys. Status Solidi RRL, 2013, 7, 50.
17. M.-H. Xie, X. Guo, Z.-J. Xu and W.-K. Ho, Chin. Phys. B, 2013, 22, 068101.
18. J. Dai, D. West, X. Wang, Y. Wang, D. Kwok, S. W. Cheong, S. B. Zhang and W. Wu, Phys. Rev. Lett., 2016, 117, 106401.
19. M. A. Herman and H. Sitter, Molecular Beam Epitaxy: Fundamentals and Current Status, Springer Berlin Heidelberg, Berlin, Heidelberg, 1989.
20. G. Zinsmeister, Thin Solid Films, 1968, 2, 497.
21. J. A. Venables, G. D. T. Spiller and M. Hanbucken, Rep. Prog. Phys., 1984, 47, 399.
22. J. A. Venables, Phys. Rev. B: Condens. Matter Mater. Phys., 1987, 36, 4153.
23. J. A. Nieminen and K. Kaski, Phys. Rev. A, 1989, 40, 2096.
24. M. C. Bartelt and J. W. Evans, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 12675.
25. G. S. Bales and D. C. Chrzan, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 6057.
26. J. G. Amar and F. Family, Phys. Rev. Lett., 1995, 74, 2066.
27. F. Family and P. Meakin, Phys. Rev. Lett., 1988, 61, 428.
28. R. Hrach and M. Sobotka, Int. J. Electron., 1990, 69, 49.
29. K. Reuter and M. Scheffler, Phys. Rev. Lett., 2003, 90, 046103.
30. K. Reuter, D. Frenkel and M. Scheffler, Phys. Rev. Lett., 2004, 93, 116105.
31. T. J. Woehl, C. Park, J. E. Evans, I. Arslan, W. D. Ristenpart and N. D. Browning, Nano Lett., 2014, 14, 373.
32. J. Zeng, P. Cui and Z. Zhang, Phys. Rev. Lett., 2017, 118, 046101.
33. V. I. Artyukhov, Y. Liu and B. I. Yakobson, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 15136.
34. R. Ouyang, J.-X. Liu and W.-X. Li, J. Am. Chem. Soc., 2013, 135, 1760.
35. P. Atkins and J. D. Paula, Physical Chemistry, Oxford University Press, Great Britain, 8th edn, 2006.
36. A. Samanta, W. E and S. B. Zhang, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 195107.
37. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15.
38. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
39. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133.
40. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188.
41. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978.
42. L. Pauling, The Nature of the Chemical Bond and the Structure of Molecules and Crystals: An Introduction to Modern Structural Chemistry, Cornell University Press, 1960.
43. C.-Z. Chang, K. E. He, L.-L. Wang, X.-C. Ma, M.-H. Liu, Z.-C. Zhang, X. I. Chen, Y.-Y. Wang and Q.-K. Xue, SPIN, 2011, 01, 21.
44. D. West, Y. Y. Sun, H. Wang, J. Bang and S. B. Zhang, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86(12), 121201.
45. J. Dai, D. West, X. Wang, Y. Wang, D. Kwok, S.-W. Cheong, S. B. Zhang and W. Wu, Phys. Rev. Lett., 2016, 117(10), 106401.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8na00136g

---