Electron g-factor in nanostructures: continuum media and atomistic approach

We report studies of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{k}}$$\end{document}k-dependent Landé g-factor, performed by both continuous media approximation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{k}}{\varvec{\cdot }}{\varvec{p}}$$\end{document}k·p method, and atomistic tight-binding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {sp}^3\hbox {d}^5\hbox {s}^*$$\end{document}sp3d5s∗ approach. We propose an effective, mesoscopic model for InAs that we are able to successfully compare with atomistic calculations, for both very small and very large nanostructures, with a number of atoms reaching over 60 million. Finally, for nanostructure dimensions corresponding to near-zero g-factor we report electron spin states anti-crossing as a function of system size, despite no shape-anisotropy nor strain effects included, and merely due to breaking of atomistic symmetry of cation/anion planes constituting the system.


Bulk material
The first part of this section contains a general introduction devoted to the linear response theory describing the Bloch state g-factors. Then, we use this approach to calculate g(k) in InAs within the eight-band k·p and sp 3 d 5 s * tight-binding models. In the last part, we introduce a simple mesoscopic model g eff (k) for InAs.
Introduction: magnetic-field dependence of the Bloch states. The electron Landé g-factor in a bulk can be calculated within the linear response theory 3,19 or using the Berry phase formula 20,21 . In the first approach, magnetic field enters as a perturbation (in the first order) to the Hamiltonian at B = 0 . For the magnetic field oriented along the z direction, this reads 3,19 where g 0 ≈ 2 , the states |ϕ ↑ / ↓� = |ϕ� ⊗ |↑ / ↓� (the orbital part is represented by |ϕ� and |↑ / ↓� is a spinor), and L z is the axial component of the angular momentum operator.
According to the Bloch's theorem, the states in a bulk crystal can be written as where α is the band index (with the corresponding energy E α (k) ). The L z matrix elements can be calculated by 5 where P i,αβ (k) ≡ �� αk |p i |� βk � and the position operator matrix elements are calculated from 22

Conduction band g-factors.
In presence of the spin-orbit interaction, the spin operator Ŝ z does not commute with the Hamiltonian and states described by some nominal spins, contains admixtures of the opposite orientation 23 . Hence, a state can not be (generally) expressed as a single product of the orbital and spin part. However, for the (mainly s-type) conduction band in InAs, the effect of the SO coupling is relatively weak in the vicinity of k = 0 . One can build superpositions for which �S z � ≈ ± /2 . To this end, the Ŝ z operator is diagonalized 19 . In the basis of the CB Bloch states this gives where a ± , b ± are coefficients found from the diagonalization and |� αk �, |� α ′ k � is a pair of CB states which are slightly splitted ( E α (k) ≈ E α ′ (k) ) due to the Dresselhaus spin-orbit coupling. The effective bulk CB g-factor calculated between such configurations is given by The Bloch states can be expressed by where c (α) m (k) are coefficients, |m� are basis states at a chosen k 0 point (in the k·p approach) or atomic orbitals (in TB models). The index m carries information about the orbital and spin part {ϕ, ↑ / ↓} . Then the momentum matrix elements can be written as where P i,nm (k) can be calculated using the Hellmann-Feynman theorem 19,22,24 where H nm (k) ≡ �n|H(k)|m� are the bulk Hamiltonian (at B = 0 ) matrix elements.
(1)  [25][26][27] where 6c , 8v , and 7v are related to the conduction-and valence band blocks (this notation corresponds to the irreducible representations of the T d point group), E g is the energy gap, {Â,B} = (ÂB +BÂ)/2 is the symmetric product, 0 is a parameter accounting for the spin-orbit interaction in the valence band, γ ′ 1−3 are the reduced (by subtracting the contributions coming from the 6c band block) Luttinger parameters, P 0 is a parameter proportional to the interband momentum matrix element; σ i are the Pauli matrices, while the explicit definitions of matrices J i , T i , and T ij are provided in Refs. 26,27 . We use the values of material parameters given in the Appendix of Ref. 14 , except for the P 0 = E p 2 /(2m 0 ) which is taken with E p = 21.2 . The remote band contributions to the electron effective mass are neglected. We also neglect the Dresselhaus SO terms and small k-linear terms related to the inversion asymmetry (the C k parameter). Finally, the Hamiltonian is diagonalized (at each considered k , separately), which gives the c (α) m (k) coefficients. To calculate the g cb (k) , one needs to find the momentum and position matrix elements constituting the L z,αβ (k) (see Eq. 2). The momentum block matrices are calculated using Eq. (6) where P y (k) and P z (k) matrices can be found from cyclic permutations. Then, the P i,αβ (k) are calculated from Eq. (5), while the R i,αβ (k) result from Eq. (3).
For k = 0 the problem is reduced to the well known case of the band-edge g-factor 2,3 , where the states in electron Zeeman doublet are purely S-type and (at B = 0 ) belong to the two-dimensional Ŵ 6c representation, giving The electron g-factor is then given by which reproduces the well-known Roth-Lax-Zwerdling formula 2 .
Tight-binding model. The tight-binding model [28][29][30][31] Hamiltonian can be written in form where N is the number of atoms, E i,α represents on-site energy, c † iα ( c iα ) is the creation (anihilation) operator of the atomic orbital α on the node i. The indices α carry also information about spin (which doubles the number of orbitals). Finally, � i,αβ accounts for the spin-orbit interaction. In the basis of |k, R n α� = e ikR n |R n α� the Hamiltonian matrix elements can be written 32 We perform the calculations using the sp 3 d 5 s * TB Hamiltonian in the nearest neighbors approach, and the t iα,jβ parameters are expressed in terms of the direction cosines 28 . The spin-orbit coupling is accounted for taking the elements � n,αβ between the p-shell orbitals 29 . We take the set of material parameters from Ref. 33 .
To obtain g(k) for the Zeeman doublet in the conduction band one need to perform a similar procedure to the case of the k·p model 6 . The momentum P i,µν (k) and the position R i,µν (k) matrix elements are calculated following Ref. 22 . In the first step, this involves then P i,µν (k) are found using the coefficients c Results. We calculated k-dependent g-factor for the Zeeman doublet in conduction band for the magnetic-field oriented along the [001] direction. To avoid spurious solutions in further k·p calculations for nanostructures 35,36 , we use a slightly reduced value of E p = 21.2 eV (which is close to E p = 21.5 eV recommended in Ref. 37 ). As shown in Fig. 1a, the results obtained from the 8-band k·p and the sp 3 d 5 s * TB model are in a good agreement. At the Ŵ point in the Brillouin zone the g-factor is about −14.4 (the k·p model) or −14.2 (the TB model). Such a strongly negative value is caused by the spin-orbit interaction, and can be predicted from the Roth-Lax-Zwerdling formula 2 . Roth formula however does not contain any k-dependence, and is generally a poor predictor of g-factors in nanostructures 38 . www.nature.com/scientificreports/ For k = 0 , the value of g-factor increases which is related to vanishing of the spin-orbit contributions. At high k , the g-factor tends to its value for free electron g 0 ≈ 2 . However, we note that for larger k values g(k) can (marginally) exceed 2 (which is visible in Fig. 1a for k > 0.9 nm −1 ). This most likely is a result of several approximations used throughout the calculation, including usage of finite basis 38,39 . To further investigate g(k) Fig. 1b,c present the cross-sections of the k·p results for g(k x , k y , 0) and g(k x , 0, k z ) planes. As in the k·p Hamiltonian we neglected terms related to bulk inversion asymmetry, the g(k x , 0, k z ) is equivalent to g(0, k y , k z ) . While, in the xy plane g-factor shows high symmetry, the xz plane exhibit strong anisotropy. This is in agreement to the fact that g-factor is more sensitive to a change of the wave-function in plane perpendicular to the magnetic field than in the field direction 38,40 .
Finally, we find that g(k) can be quite well approximated by the following mesoscopic formula where the parameters α 0 = 16.2 , β 2 xy = 30 nm 2 , β 2 z = 11 nm 2 were fitted to the TB results, and in particular α 0 is chosen such that in a limit of k → 0 , g eff → −14.2 . As shown in Fig. 1a, this formula gives a good agreement to the exact results, especially for smaller k values. β xy = 5.48 nm and β z = 3.32 nm have a unit of length, in an very loose analogy to the of concept magnetic length for Landau levels (equal to 25.6 nm at 1T). Moreover, β xy = β z indicating lack of equivalence of x,y and z as the field is applied and oriented along z ([001]) direction.

Nanostructures
In the first part of this Section we describe the implementations of magnetic field to nanostructure modeling within the k·p and the TB approaches. Then, we compare the numerical results obtained from both methods for an electron confined in a three dimensional InAs box (cube), as a function of box size (edge length) varying from a single lattice constant, i.e. 0.6 nm to over 120 nm and over 60 million atoms involved in the computations. We show, that these results can be qualitatively well reproduced using an effective model based on the bulk g-factor g(k) with virtually no computational cost. We also discuss the effect of symmetry breaking due to underlying crystal lattice, which (in presence of the SO coupling) allows on the mixing of spin configurations in the Zeeman doublet.
Eight band k.p model. In a standard way, the magnetic field enters the k·p Hamiltonian by the substitution k → k + (e/ )A . Since, a straightforward implementation leads to gauge dependent results 41 , the gaugeinvariant scheme was developed 42,43 . Furthermore, the k·p Hamiltonian is supplemented by the magnetic terms where 27 is related to spin. The band angular momentum (which was neglected in the k·p calculations in the previous section) is represented by Finally, the H (r) describes contributions from the remote bands, which are not explicitly included in a given k·p model 19 . As the calculated electron band-edge g-factor ( g = −14.4 ) is already close to the experimental value ( g = −14.7 4 ), we neglect these remote band contributions.
Tight binding model. The magnetic field is implemented to the TB Hamiltonian by the standard Peierls substitution 34,39 with a phase given by Assuming constant magnetic field and a symmetric gauge, one can obtain 44 The contribution from the spin is accounted for via the on-site terms 34 (12)  For nanostructures, in order to avoid any spurious states in the energy range of interest, the dangling bonds of surface atoms are passivated according to well established approach of Ref. 45 . We found that tight-binding results somewhat depend on the choice of the dangling bond shift, however the trends obtained with different shifts are very similar, in agreement with conclusions of Ref. 11 . Strain effects and piezoelectricty are not present in a system due to use of a single chemical compound. For sake of comparison with the k·p and the effective model, any surface reconstruction or presence of image charges are also neglected.
Results: electron in a box. In the following, we calculate the electron g-factor for nano-size InAs cubic box, and by varying box size simultaneously in all three dimensions. We start with a box size of one lattice constant edge length, corresponding to very high quantum confinement, and we progress through intermediate dimensions to very large InAs boxes, with dimensions extending well over 100 nm aiming for near bulk-like properties. The largest of considered box (200 × 200 × 200 in lattice constant units) has edge length equal to 120 nm, corresponding to 64.2 million atoms in the calculation. In this case, the ground electron state energy is 0.422 eV, thus only 3 meV larger than the bulk limit of 0.418 eV. We note as well, in largest considered cases, the tight-binding calculation was performed on parallel computer cluster with 192 computational cores, and MPIparallelized Lanczos solver, taking approximately 12 h to find several lowest electronic states, and the time of the computation scaling proportionally with number of atoms for smaller systems 46 .
The numerical results are presented in Fig. 2, where the values of the g-factor are taken from the differences between the first excited ( E 2 ) and the ground ( E 1 ) state energy levels g = (E 2 − E 1 )/(µ B B) , calculated at B = 1 T. The sign is determined from the spin orientations of states, which are represented by the averages S z . In the tight-binding model, the average z-th spin component for the ψ n state can be calculated as where i describes the atomic sites and indices α , β denote the atomic orbitals.
The results obtained from the TB and k·p models provide qualitatively similar g-factor dependence. In both cases, the increasing size brings g-factor to the bulk value in a similar fashion.
Moreover, in both cases for a very strong confinement the value of the g-factor is close to the free electron case ( g 0 ). This dependence on the system size can be connected to increasing (with confinement) effective energy gap (see Eq. (18) below) 11 and to angular momentum quenching 38 . Yet, there is a notable difference between both approaches, with k·p systematically reporting smaller g-factor magnitudes. However, for a fairness of comparison, quantitative differences of that kind are expected, since we compare 20-band atomistic model (including d-orbitals), and 8-band method based on continuous media approximation. Both method report somewhat different bulk band structures, utilize different of treatment of the cube surface, boundary conditions etc.
We also performed the simulations in a simple effective model which relies on the g eff (k) for bulk InAs (see Eq. (12)) and on the Fourier transform of the wave-function. For sake of simplicity, and computational efficiency, we assume that the electron ground state is build from the bulk conduction band and having the envelope of the three-dimensional infinite well (with L × L × L size) Then, the effective g-factor can be calculated from As shown in Fig. 2, this surprisingly simple, effective model provided results which are in-between k·p approach, and the multi-million atom simulations obtained using the TB model. By performing several numerical tests, we found that the g box (L) differs tight-binding prediction mostly due to oversimplified Eq. (16) assuming hard-wall boundary conditions, whereas approximation of g eff plays a lesser role.
It is known that electron g-factor for any nanostructure can be decomposed into an isotropic part (depending on the effective band gap) and the surface part which depends on the shape of the structure 11 . The first part can be expressed by 11 where E is the effective band gap, and g 0 , E 0 are parameters which can be fitted for a given material. We used this formula for calculations in the box (the red line in Fig. 2). Although for the parameters g 0 = 3.09 and E 0 = 1.69 eV (which we fitted to the TB results) we obtain a very good agreement in middle part of the plot, we cannot simultaneously reproduce the values at the boundaries (in the limit of a very small/large box). This suggest, that in the considered system the part of the g-factor related to the surface plays important role.
Role of atomistic symmetry. So far we have studied nanostructures without analyzing the role of underlying crystal lattice 48,49 . We note that nanostructures analyzed in Fig. 2 were cubic boxes cut from zinc-blende lattice, obtained by terminating box sides with the same ionic species Fig. 3a. Alternatively such geometry can be obtained by choosing the box geometric center to be placed on one of ions. In case of Fig. 3a this center is placed on anion (marked as red spheres), and results discussed so far were obtained for this particular choice.
Such atomic arrangement leads to an overall T d (tetrahedral) symmetry of a nanostructure, which may not be apparent from inspection of Fig. 3a, however it was additionally verified with Jmol 50 and Chemcraft 51 tools allowing for point-group symmetry determination.
Moreover, for the same T d symmetry, one can choose either cation or anion as an origin. Centering the system on the cation effectively corresponds of replacing anions with cations [or red with blue atoms in Fig. 3a]. In both cases (whether centered on anion or cation) the overall symmetry is T d provided that box is terminated consistently with one ionic species only. As we consider magnetic field oriented in the [001] direction, the field effectively lowers the symmetry of the system from the T d to the S 4 point group, therefore the tight-binding results in Fig. 2 corresponded to S 4 symmetry. www.nature.com/scientificreports/ Even though not altering the symmetry, switching indium atomic positions with arsenics, may substantially affect single particles energies, however only for smaller boxes, as shown in Fig. 4. These differences are related to different ionic composition 33 of ground electron states, since nanostructure shown in Fig. 3a will have different number of cation and anion depending on the choice of origin. For the anion-centered case (and anions terminating the surface), with overall number of anions is larger than cations. For the cation-centered case the situation is exactly opposite. Centering the origin on the anion or cation will thus not alter symmetry, but will change the overall ionic "stoichiometry". Nonetheless, despite substantial differences in single particle energies for smaller boxes, anion-cation centering issues have a rather small effect on the final g-factor dependence (inset in Fig. 4). Interestingly, a small difference in g-factor values between anion and cation centered cases can be observed even for larger boxes, again due to difference in stoichiometry, however it decays with system size, and for the edge length of 60 nm (100 × 100 × 100 box; including 4 ×10 6 cations and 4.06 ×10 6 anions) the g-factor varies by about 1% with respect to cation-centered case (with 4.06 ×10 6 cation and 4 ×10 6 anions).
Another possible choice of crystal lattice arrangement with respect to box shape is presented in Fig. 3b, with corresponding tight-binding results shown for comparison in Fig. 4. This particular case is given by a seemingly unimportant shift of the box origin by half of a bond length, i.e. placing it in the exactly between anion and cation. Such choice also leads to termination of box sides in a mixed anion-cation fashion as apparent from Fig. 3b. Moreover, mid-bond case corresponding to exactly equal number of cations and anion, for all considered dimensions. Although there is not strict inversion symmetry in zinc-blende lattice, further replacing anions with cations in Fig. 3b leads to virtually the same single particle spectra (with µ eV differences), therefore it is not considered here.
However, despite ideal stoichiometry, a closer inspection reveals that anion/cation mixing at the boundaries leads to an overall symmetry reduction to C 3v , and with further reduction of symmetry to C 1 , when the magnetic field along z-axis is applied. Yet, as seen in Fig. 4 there is virtually no difference to previous results, with the exception of near-zero g-factor values.
In the previous case of the high symmetry system (whether anion or cation centered), the g-factor value smoothly varied with the box size and crossed 0 (which simply corresponded to a swap of the states in the Zeeman doublet). In contrast, for the low symmetry box the g-factor changes its sign avoiding 0 (while the absolute value remains continuous). This is better illustrated in Fig. 5a where absolute value of g-factor is now presented for mid-bond (low symmetry) case only, and with a magnification of close to g-zero regions. This peculiar feature in g-factor dependence corresponds to an anti-crossing between the states in the Zeeman doublet, which is reflected in the values of S z (Fig. 5b). This behavior is related to the symmetry selection rules. According to the group theory, two states can be mixed (which manifests in an anti-crossing of their energy levels) if they belong to the same irreducible representation 23,52 . In the case of the same-ion-terminated box, the symmetry of the system is described by the S 4 point group (which in presence of the spin-orbit coupling needs to be a double group). Then, the two lowest electron states belong to different irreducible representations and coupling between them is prohibited. In contrast, for the mixed-anion-cation-terminated box, the symmetry is reduced to the C 1 (which including spin is also a double group). In this case, the states in the Zeeman doublet belong to the same representation. Hence, their energies exhibit anti-crossing and spin configurations mix (which is in fact caused by the spin-orbit coupling). This effect will thus happen for any low atomic arrangement, and is not limited to a particular choice presented in Fig. 3b, although application of magnetic field exactly along high symmetry crystal axes, such as [111] (diagonal axis in Fig. 3b) would restore high-symmetry, and allow for electron level crossing rather than anti-crossing. However, no nanostructure can ever be grown to have an ideal atomistic symmetry as systems presented in Fig. 3. In fact, removal of just a single atom can break the overall symmetry and thus will lead to level anti-crossing. The effect of anti-crossing could be also viewed as a presence of off-diagonal terms in the g-tensor. It is known that g-tensors in low symmetry systems (the C 1 group without magnetic field) have 9 independent components 7 , which cannot be reduced to the diagonal form for any magnetic field direction.

Conclusions
In summary, with atomistic tight-binding method, and continuous media approximation k·p approach we have studied Landé g-factor k-dependence for InAs. Based on that, we proposed a mesoscopic model with three effective parameters only, and further we compared this model with multi-million atom tight-binding calculations, for a cubic nanostructure with dimensions varying from single to 120 nanometers and number of atoms reaching 64 million. Despite its simplicity the mesoscopic model shows a good qualitative dependence with actual results ranging between that of the tight-binding and the k·p at virtually no computation cost. Further, we have inspected nanostructure dimensions corresponding to near-zero g-factor, and found that depending on a detail of atomic arrangement electron spin states can undergo an anti-crossing as a function of system size. The effect occurs despite high-shape symmetry, and is not related to cation-anion stoichmetric imbalances, but due to low overall symmetry being a result of cubic shape imposed on underlying zinc-blende lattice. Our results therefore emphasize the key role of symmetry in nanostructures, and show inherent limits to g-factor tuning, especially important for applications involving near-zero g-factor values.