Nagaoka spin-valley ordering in silicene quantum dots

We study a cluster of quantum dots defined within silicene that host confined electron states with spin and valley degrees of freedom. Atomistic tight-binding and continuum Dirac approximation are applied for few-electron system in quest for spontaneous valley polarization driven by inter-dot tunneling and electron-electron interaction, i.e. a valley counterpart of itinerary Nagaoka ferromagnetic ordering recently identified in GaAs square cluster of quantum dots with three excess electrons [P. Dehollain, {\it et al.}, Nature {\bf 579}, 528 (2020)]. We find that for Hamiltonian without intrinsic-spin orbit coupling -- similar to the one of graphene with staggered potential -- the valley polarization in the ground-state can be observed in a range of inter-dot spacing provided that the spin of the system is frozen by external magnetic field. The inter-valley scattering effects are negligible for cluster geometry that supports the valley polarized ground-state. In presence of a strong intrinsic spin-orbit coupling that is characteristic to graphene no external magnetic field is necessary for observation of ground-state that is polarized in both spin and valley. The effective magnetic field due to the spin-orbit interaction produces a perfect anticorrelation of the spin and valley isospin components in the degenerate ground-state.


I. INTRODUCTION
The Hubbard model for cubic lattice with on-site Coulomb interaction dominating over the inter-site hopping produces spin polarized ground-state near halffilling [1]. In the theory of itinerary ferromagnetism [2] the phenomenon is referred to as Nagaoka ferromagnetism. Semiconductor quantum dots were pointed out as a possible two-dimensional realization of the Hubbard model and artificial molecules or clusters formed by multiple quantum dots were studied in the context of a spin polarization driven by inter-dot tunneling and electronelectron interaction [3,4]. Spin-ordered ground state was recently experimentally identified in electrostatic quantum dots defined in GaAs [5] in a three-electron system for quantum dots arranged in a square cluster, a case previously theoretically studied in Ref. [3].
In graphene [6] and graphene-like Xene materials [7] the electron states near the charge neutrality point are additionally characterized by the valley isospin due to the presence of two non-equivalent Dirac points in the Brillouin zone. Electrostatic confinement in graphene is excluded by the Klein tunneling effect [8]. However in bilayer graphene [9] or buckled silicene [10] the energy gap and thus the electrostatic confinement can be formed by perpendicular electric field. The valley degree of freedom of confined states can be resolved by the transport spectroscopy [11,12] in a similar way as the spin is [13]. A possible Nagaoka ordering of the valley would contribute to the toolbox of valleytronics [14,15] for confined states in two-dimensional materials.
In this paper we look for the counterpart of the Nagaoka ferromagnetism in the valley degree of freedom in a 2D system. We focus our attention on the groundstate valley polarization for a three-electron system in a square cluster of quantum dots defined in silicene [7], i.e. a counterpart of the case studied experimentally in GaAs system [5]. With respect to the III-V quantum dots the silicene besides the valley degree of freedom hosts a strong intrinsic spin-orbit coupling [16,17] which as we show below can play a role in the Nagaoka ordering. Without the spin-orbit coupling term the tight-binding Hamiltonian is identical with the one for the monolayer graphene with staggered potential [18][19][20] up to the numerical value of the inter-atomic hopping energy. For that reason below we solve both the problems with and without the spin-orbit coupling. We demonstrate that in the absence of the spin-orbit coupling the spin degree of freedom of the three-electron system needs to be frozen for the valley ordering to be observed. The intrinsic spinorbit coupling splits the fourfold degeneracy of the confined single-electron ground state with respect to the spin and valley forming spin-valley doublets in a manner similar to the one closely studied for carbon nanotubes [21]. We find that in presence of the intrinsic spin-orbit coupling Nagaoka ordering in both valley and spin appear simultaneously.

A. Atomistic tight-binding Hamiltonian
For the tight-binding model we define a flake of buckled silicene [22,23] with ions of the A sublattice at positions r A k = k 1 a 1 + k 2 a 2 with the crystal lattice vectors a 1 = a 1 2 , √ 3 2 , 0 , and a 2 = a (1, 0, 0), with the silicene lattice constant a = 3.89 Å. The B sublattice is shifted by a base vector r B k = r A k + (0, d, δ) where d = 2.25 Å is the in-plane nearest-neighbor distance and the vertical shift arXiv:2012.12524v1 [cond-mat.mes-hall] 23 Dec 2020 of the sublattices is denoted by δ = 0.46 Å. Calculations are performed for a hexagonal flake with armchair edges and side length of about 30 nm with approximately 72 000 p z spin-orbitals.
We determine the eigenstates of atomistic tightbinding Hamiltonian [16,17,23], where the first sum describes the nearest neighbor hopping with the energy t = 1.6 eV [16,17]. In Eq. (2) p kl = exp i e ´ r l r k A · dl stands for the Peierls phase. The integral in the exponent of p kl accounts for the Aharonov-Bohm phase shifts that the wave functions acquire from the vector potential A in hopping. We consider the magnetic field with both perpendicular B z and in-plane B x components B = (B x , 0, B z ) with the vector potential A = (−B z y/2, B z x/2 − B x z, 0). Due to the 2D nature of the material the in-plane component does not produce noticeable orbital effects. The in-plane field is introduced in order to manipulate spins of the confined states via the spin Zeeman effect included in the last term in Eq. (2), with µ B as the Bohr magneton and g = 2 as the Landé factor. The third sum in Eq. (2) introduces the intrinsic spin-orbit interaction [24] with the coupling constant t SO = 3.9 meV [16,17] and ν kl = +1 (−1) for the path of the next-nearest neighbor hopping from ion l to k via the common neighbor that turns counterclockwise (clockwise). The second sum in Eq. (2) introduces the external potential with V k standing for the potential on r k ion.
In silicene the bias between sublattices opens the energy gap in the band structure [10] that allows for formation of the confinement potential. For the confinement potential we assume that the bias is independent of the electron position within the plane, where the sum over i runs over 4 quantum dots with centers g i and r ik = |r k − g i |. We take R = 4.2 nm for the dot radius, and V g = 0.3 eV for the depth of potential cavities. In Eq. (2) the electrostatic potential is lowered on both sublattices near a center of each quantum dot. This type of potential variation -with nearly equal bias and a minimum on both sublattices -can be achieved with a pair of flat gate electrodes with one that contains a circular intrusion near the quantum dot center [25]. Note, that type I quantum dots can also be produced with flat gates provided that they contain apertures near the confinement area [26]. Three-electron charge densities for the centers of the quantum dots g i forming a square of side length X are given in Fig. 1 with the A (B) sublattice placed on the left (right) column of plots, and the side of the square increasing from top to bottom.

B. Continuum Hamiltonian
The atomistic calculations are supported with modelling of the continuum approximation that explicitly resolves the valley degree of freedom. We work with a fourcomponent wave function spanned on sublattice and spin subspaces ψ = ψ A↑ ψ B↑ ψ A↓ ψ B↓ T , and the energy operator [16,27] where σ and τ are the Pauli matrices in the spin and sublattice subspaces, respectively. I sublattice and I spin are the identity matrices, the Fermi velocity is v F = 3dt/2 . The wave vector operators are defined as k = −i∇ + e A, and η = ±1 is the valley index.
Hamiltonian (3) is diagonalized by the finite element method. The computational box is divided into typically about 2200 triangular elements with 18000 nodes supporting Lagrange interpolating polynomials of the second degree [28] as shape functions covering the spin and sublattice spaces. In Eq. (3) the last expression is an artificial Wilson term [29] that is applied to remove the spurious states [29][30][31][32] due to the fermion doubling problem from the low-energy spectrum. We take the Wilson parameter W D = 36 meV nm 2 which increases the energy of the fast oscillating states with a negligible influence on the actual solutions of the Dirac equation that are smooth near the charge neutrality point.

C. Calculations for three electrons
We look for signatures of Nagaoka ordering effects in the system of three excess electrons in a quadruple quantum dot. The electron-electron interaction for gapless graphene flakes can lead to generation of electron and hole pairs so that only the confined charge and not the number of electrons is a good quantum number [33]. In the calculations that follow for the three-electron system the typical total interaction energy is about 60 meV, i.e. 20 meV per electron pair, i.e. much lower than the potential bias between the sublattices. Since the interaction energy is lower than the energy gap and the considered quantum dot does not support confinement of holes we neglect the effects of pair-generation by Coulomb interaction [33] and assume that the number of conduction band electrons is fixed [34]. In both the atomistic method and in the continuum approximation we diagonalize the Hamiltonian in the basis of three-electron wave functions constructed by the lowest-energy 48 confined eigenstates of the conduction band that produces the basis of 17 296 Slater determinants.
The Hamiltonian for the system of interacting electrons is where d † i is the electron creation operator for the energy level E i . The two-electron Coulomb matrix elements that with κ = e 2 /(4π 0 ). We take 0 = 4.5 for the dielectric constant, that corresponds to SiO 2 or thin layers of Al 2 O 3 [35,36] applied as a matrix embedding the silicene monolayer.

Coulomb integrals in the continuum approach
The integration of the Coulomb matrix elements is performed in a somewhat different manner for the continuum approximation and in the atomistic approach. For the continuum approach we calculate the Coulomb integrals using the formula where the four components of each function are calculated by the finite-element method. The Kronecker deltas with the valley indices that stand before the integral imply the neglect of the inter-valley scattering effects by the short-range component of the Coulomb interaction [37][38][39]. The tight-binding approach accounts for the intervalley scattering and the comparison of the two methods that follows provides estimation of the inter-valley scattering effects.

Coulomb integrals in the atomistic approach
In the tight-binding method the single-electron wave functions ψ are expanded in the basis of 3p z spin-orbitals of Si ions, The Coulomb matrix elements are summed over the ions, We apply the two-center approximation [34] for the integral. p a z ( The single site integral (a = b) is calculated for 3p z Si atomic orbitals with p z (r) = N z 1 − Zr

III. RESULTS AND DISCUSSION
A. Single-dot single-electron results Before we proceed to the discussion of the results for three electrons in the multiple quantum dot we first present the spectrum for a single electron in a single quantum dot. In the absence of spin-orbit interaction and without external magnetic field the single-electron confined ground state is four-fold degenerate with respect to both spin and valley [ Fig. 2(a)]. The splitting of energy levels by the perpendicular magnetic field B z is stronger with respect to the valley than to the spin [ Fig. 2(a)]. The valley degeneracy is preserved for the in-plane field B x and the energy levels are split only with respect to the spin [ Fig. 2(b)].
The intrinsic spin-orbit coupling is diagonal in zcomponent of the spin. The coupling introduces an effective magnetic field perpendicular to the plane of confinement with orientation that is opposite in the sense of eigenvalue sign to the valley isospin η (cf. the term with t SO coupling constant in Eq. (3)). This effective field For a circular quantum dot the single-electron ground state corresponds to angular momentum quantum number 0 for the wave function component on sublattice A and ±1 on sublattice B [25]. In consequence, the charge density in a single quantum dot corresponds to a maximum and a zero of the charge density in the centers of the quantum dots on A and B sublattice, respectively. The single-electron properties are consistent with the charge density distribution for three interacting electrons with well separated wave functions -see Fig. 1(e,f) for the centers of potential minima distributed on a square with side length of X = 11.7 nm. The wave function component on the B sublattice is less strongly localized and thus it mediates the inter-dot tunneling in a stronger extent.

Nagaoka valley ordering
The three-electron spectrum is given in Fig. 3(a) as calculated with the tight-binding approach and in Fig.  3(b) as obtained with the finite-element method. The results of the atomistic and continuum approaches agree very well up to a relative shift of the entire spectra on the energy scale of a few meV. In the results of the atomistic approach we plot the energy levels with colors indicating the total spin z or x component. On the energy levels calculated with finite-element modeling we mark by the color of the lines the total valley isospin component for the three electrons.
The three-electron ground-state at B z = 0 is four-fold degenerate. The degenerate energy levels correspond to the eigenvalue of S z = 1 2 3 i=1 σ i z z-component of the total spin ± 1 2 and the z-component of the total valley isospin of Σ KK = 1 2 3 i=1 η i z equal to ± 1 2 . The groundstate is therefore not polarized neither in spin nor in valley and no Nagaoka ordering is found. The first excited state is sixteen-fold degenerate. The energy levels that are degenerate at B z = 0 correspond to the total valley index Σ KK = − 3 2 , − 1 2 , 1 2 , 3 2 . For each of these values an energy level with each values of the z-components of the spin S z = − 3 2 , − 1 2 , 1 2 , 3 2 is found. Nagaoka ferromagnetism was observed [5] in a quadruple quantum dot defined in GaAs, with three electrons per 8 available spin-orbitals. In silicene system we have 16 available spin-valley-orbitals due to the additional valley degree of freedom. One can try to eliminate the spin degree of freedom and reduce the number of equivalent states to 8 by applying a strong magnetic field to freeze the spin degree of freedom. The application of the perpendicular magnetic field lifts the valley degeneracy [ Fig.  2(a)] that would eventually lead to the valley polarization by the external field. Here, we want to preserve the valley degeneracy in order to study the electron-electron interaction triggering the valley polarization. For that reason we choose to apply the in-plane field that interacts only with the spin and not the valley of confined states [ Fig.  2(b)]. The results are given in Fig. 3(c,d) with B z = 10 mT applied to slightly split the energy levels. The states that are spin polarized in the −x direction for all valley configurations are promoted to the lower energy by the B x field.
The structure of the low-energy spectrum spinpolarized by strong B x field is revealed once a weak B z field is additionally applied. In Figure 4 we set B x = 20 T and calculate the energy levels as functions of B z . In the applied range of B z ≤ 0.3 T the spins remain nearly perfectly polarized in the −x direction. The groundstate at B z = 0 is fourfold degenerate with the valley isospin component that takes values − 3 2 , − 1 2 , 1 2 and 3 2 . This structure of the ground-state energy level is a valleyordered counterpart of spin-polarized three-electron state with the total spin quantum number S = 3 2 and the spin component quantum numbers S z = − 3 2 , − 1 2 , 1 2 and + 3 2 . In the excited part of the spectrum in Fig. 4 we find a number of valley non-polarized states forming doublets at B z = 0 with the valley isospin equal to ± 1 2 . By analogy to the spin degree of freedom we attribute the total valley isospin quantum number V = 3 2 to the four-fold degenerate state with the isospin component Σ KK changing from −V to V with steps of 1. The doublets thus correspond to V = 1 2 with the components Σ KK = ± 1 2 . The Nagaoka polarization of the ground-state appears due to the inter-dot electron tunneling and thus it is determined by the system geometry. In order to study the valley polarization we kept B x = 20T and varied the positions of the centers of the dots. We define ∆E 31 as the energy difference between the lowest valley polarized state with V = 3 2 and the lowest non-polarized state with V = 1 2 state. ∆E 31 < 0 corresponds to the Nagaoka ordered valley in the ground-state. The black line in Fig.  5 shows the result for the continuum Hamiltonian as a function of X -the side length of the square on which the centers of the dots are placed [ Fig. 1]. The valleypolarized ground-state is found for X 10.5nm. The polarized ground-state is most stable for X = 11.7nm that was selected for plots in Fig. 1(e,f), Fig. 3 and Fig. 4. For low values of X the four quantum dots become more strongly tunnel coupled and they eventually are transformed to a single quantum dot for which the ground-state is not polarized. In the limit of large X the inter-dot tunneling is negligible. Once the electrons are separated and their hopping between the dots is removed the valley isospin has no influence on the energy for B z = 0, hence the degeneracy of V = 3 2 and V = 1 2 states at large X. The black lines shows the energy difference between the lowest valley polarized and non-polarized states as obtained with the continuum approach for Bx = 20 T and tSO = 0 as a function of the side length of a square on which the centers of the four quantum dots are placed [see Fig. 1]. The red line shows the results obtained for X = 11.7 nm, that corresponds to the most stable valley-polarized ground state as a function of a shift ∆x of the y position of the quantum dot localized in the first quadrant of the Cartesian coordinate. The insets show the potential on the A sublattice for X = 11.7 nm and ∆X = 0 (left) and ∆X = 4 nm (right). Same scale for VA is applied for both sublattices. The frame in the insets has the length of 60 nm.
Ref. [5] found that the ground-state spin ordering vanishes when the array of dots is deformed to approach the limit of a quantum dot chain. For a chain of dots [41] the Lieb-Mattis theorem [42,43] excludes spin polarization of the ground state. Here we looked for a similar effect in the valley degree of freedom. We fixed the value of X to 11.7 nm and moved the quantum dot localized in the first quadrant of the coordinate system (x > 0, y > 0, see the insets to Fig. 5) shifting its center by ∆X in the y direction. The results are displayed in Fig. 5 by the red line. The shift by ∆X 0.75 nm makes the valley polarized and non-polarized states degenerate. For ∆X 2.1 nm the value of ∆E 31 is more or less inverted from the ∆X = 0 case. The energy spectra for X = 11.7 nm and ∆X = 2.1 nm are plotted in Fig. 6. The ground-state here is a V = 1 2 doublet and the valley polarized V = 3 2 quadruplet appears as an excited state nearly degenerate at B z = 0 with an excited doublet. The energy spacing between the quadruplet and the excited doublet is larger for the atomistic model than in the continuum approach [see Fig. 6(a) and Fig. 6(b)].

Inter-valley scattering effects
Results for the energy spectra presented so far were very similar in the continuum and tight-binding approaches. In the continuum approach that we apply here all the inter-valley scattering effects are neglected. The inter-valley scattering can appear as (i) a single-electron effect triggered by the armchair edge [40] of the flake and (ii) as an interaction effect due to the short-range component of the Coulomb potential [37][38][39].
The effect (i) can be observed for a smaller flake, when the tails of quantum-dot-confined wave functions tunnel to the armchair edge that induces inter-valley mixing. Figure 7 shows a zoom of the low-part of the spectra for parameters applied in Fig. 4 with the side length of the flake of 20.5 nm [ Fig. 7(a)], 23 nm [ Fig. 7(b)] and 25 nm [ Fig. 7(c)]. The valley mixing due to the edge effect lifts the degeneracy at B z = 0 and opens avoided crossings between energy levels that in the continuum approach correspond to different valley isospin quantum numbers [compare with Fig. 4(b)]. For strong inter-valley mixing the dependence of energy levels of B z deviates from linear, in particular near B z = 0.    The effect (ii) is not triggered when the electrons occupy separate quantum dots. In this case only the long range tail of the Coulomb potential is resolved by carriers. The applied method for the continuum approach explicitly neglects the inter-valley scattering. The scattering due to interaction can be introduced to the continuum Hamiltonian, although the procedure is far from straightforward [37][38][39]. The atomistic approach applied here inherently accounts for the inter-valley scattering. As we find above, the continuum Hamiltonian produces results that are very close to the tight-binding spectra for the geometry of the system for which Nagaoka valley ordering is found. In Fig. 8 we plotted the results for quantum dot centers placed at the corners of the square of side length X = 10 nm. For this parameters the valley ordering is no longer observed in the ground-state [see Fig. 5]. In Fig. 8 we can see that the agreement between the two methods is no longer perfect. In particular, the valleyordered quadruplet is found here only by the continuum approach but in the atomistic model the quadruplet is split into two doublets. The two-electron levels splitting by the inter-valley scattering due to the electron-electron interaction was discussed in detail in Ref. [25]   electron pair in a single quantum dot (see Fig 4(a) and Fig. 4(b) in [25]). In the parameter range where Nagaoka ordering is found the inter-valley scattering by the Coulomb interaction is negligible or missing.
C. Nagaoka polarization in presence of the spin-orbit interaction Figure 9 shows the spectra for X = 11.7 nm and t SO = 3.9 meV in the absence of the in-plane field B x = 0. We notice that with the spin-orbit interaction the agreement between the continuum and atomistic spectra is still very good. Moreover, the pattern of energy levels is similar to the one found in Fig. 4 for t SO = 0 and B x = 20 T. For t SO = 3.9 meV and B = 0 a single dot hosts a two-fold degenerate ground state [ Fig. 2(c,d)] with opposite spin and valley isospin components ησ z = −1. This pattern of energy levels replaces the valley-degenerate ground-state with η = ±1 found for t SO = 0, B z = 0 and strong in-plane field that freezes the spin found for Fig. 2(b). In these two above cases the Coulomb integrals between the single-electron states in the low-energy part of the spectrum are similar, hence the agreement of the three-electron spectra in Fig. 9 and Fig. 4. In presence of the intrinsic spin-orbit interaction there is a perfect anticorrelation between the S z and Σ KK quantum numbers [ Fig. 9(a,b)]. All the information on the eigenstates is therefore redundantly included in the spin and valley sets of quantum numbers. The above conclusions hold also for the deformed system with ∆X = 2.1 nm that are plotted in Fig. 10.

IV. SUMMARY AND CONCLUSIONS
We studied a system of three-electrons in a square cluster of quantum dots defined within material that provides valley degree of freedom to the confined singleelectron states using the continuum approach that allows for identification of the valley isospin in the atomistic tight-binding spectra. We found that the Nagaoka-type polarization of the valley in a system without the intrinsic spin-orbit coupling is found in conditions when the spin degree of freedom is frozen by in-plane magnetic field. Non-polarized ground state is promoted when the spatial symmetry of the cluster is lifted by a shift of one of the quantum dots that transforms the cluster toward a chain of quantum dots. In presence of the intrinsic spin-orbit coupling the spin-valley polarization is observed along with the perfect anticorrelation of the spin and valley isospin components already in the absence of external magnetic field. The pattern of the energy levels near the ground-state for systems with and without the spin-orbit coupling is very similar provided that a strong in-plane magnetic field is applied to the latter. The spontaneous ground-state valley polarization in the system can be harnessed for studies of valley manipulation in multiple quantum dots.