Nagaoka ferromagnetism in an array of phosphorene quantum dots

We consider an array of four quantum dots defined in phosphorene containing three excess electrons, i.e., in the conditions of near half filling when itinerant Nagaoka ferromagnetism is expected to appear in a square array with isotropic interdot hopping. The interdot hopping in the array arranged in a square inherits the anisotropy from the form of the phosphorene conduction band. We apply the configuration interaction method for discussion of the appearance and stability of the spin-polarized ground state and discuss the compensation of the effective mass anisotropy by the geometry of the quantum dot array. Our study shows strong stability of Nagaoka ferromagnetism for optimized geometry of the array, with the Nagaoka gap as large as ∼ 230 µeV. A phase diagram for the ground-state spin ordering versus the geometric parameters of the array is presented. We study the suppression of the ferromagnetism in a transition of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\times 2$$\end{document}2×2 array to a quasi-1D chain and indicate that the shift of one of the quantum dots away from the array center is enough to transform the system to a quantum dot chain. A shift in the zigzag crystal direction induces the low-spin ground state more effectively than a shift along the armchair direction. We also discuss the robustness of the spin ordering against detuning one of the dots. The ferromagnetic ground-state survives as long as the detuning is not large enough to trap one of the electrons within a single quantum dot (for positive detuning) or remove one of the quantum dots of the accessible energy range (for negative detuning).


Single electron
An electron in the phosphorene monolayer can be described with the single-band continuum Hamiltonian 24,27 , with effective masses m x = 0.17037m 0 along the armchair direction (see Fig. 1) and m y = 0.85327m 0 along the zigzag direction, derived using a tight-binding Hamiltonian in Ref. 24 and V(x, y) is the confinement potential.The effective masses were obtained by fitting the energy spectra of the single-band continuum Hamiltonian to the tight-binding model for the harmonic oscillator 24 and the ring-like confinement potential 27 .The magnetic field is applied perpendicular to the plane of 2D phosphorene crystal for non-zero Zeeman term and is introduced with symmetric gauge A = (A x , A y , A z ) = (−yB/2, xB/2, 0) .In Eq. (1), we use the Landé g-factor, g = g 0 = 2 .The value is taken from Ref. 28 where k • p theory was used to extract the Landé g-factors.The purpose of including the Zeeman term in the Hamiltonian is to remove the degeneracy of spin eigenstates.A range of g values were reported [29][30][31][32] .The specific value of g changes the magnitude of the magnetic field required to lift the degeneracy.
(1) www.nature.com/scientificreports/ We introduce the external magnetic field for a residual lifting of the spin degeneracy which is useful for the characterization of the spectrum.In the discussion of Nagaoka ferromagnetism, for convenience, we set a residual magnetic field B = 1 mT to lift the degeneracies with respect to the z component of the spin σ z .The Nagaoka gap measured experimentally for the GaAs system was of the order of a few µeVs 11 .We expect our system to have a much larger gap for the strongest ferromagnetic state.The magnetic field produces splitting between states of the z component of the total spin σ z = −3/2 and σ z = −1/2 of about 0.22 µeV.As we will see in further results, this energy will be small compared to the largest Nagaoka gap, and will affect the ground state only when the Nagaoka gap is of the same order (near the phase transition) and not the results near the point of largest gap.The energy gap is calculated as �E = E 3/2 − E 1/2 where E 3/2 is the energy of the spin quantum number S = 3/2 state and E 1/2 , the energy of the spin S = 1/2 state.The Hamiltonian in Eq. ( 1) is solved using the finite difference method with gauge-invariant discretization 33 .We use a square mesh with spacing x in both directions so that the action of Hamiltonian on the wave-function � α,β = �(x α , x β ) = �(α�x, β�x) is given by with Peierls phases 33 C x = exp(−i e �xA x ) and C y = exp(−i e �xA y ) to account for orbital effects of the mag- netic field.The Hamiltonian in Eq. ( 2) is diagonalized in a finite computational box of size 50 nm × 50 nm and the spacing x of 0.4 nm.The numerical diagonalization is carried out using the Lanczos algorithm as in Ref. 24 .

Model potential
Inside the computational square box of side length 50 nm, the model quantum dots are arranged in a rectangle with distance between the centers of dots 2µ x in the armchair direction (x) and 2µ y in the zigzag crystal direction (y).The spacings influence both the interdot tunneling rates (hopping energy) and the interdot electron-electron interaction.The applied model confinement potential is such that the centers of the dots are located at (±µ x , ±µ y ) .Each dot is determined by a Gaussian potential with the size given by the parameter s and the depth of the potential by V d .Throughout the work, the size of the dots is taken to be s = 7 nm.The chosen soft Gaussian potential is a good approximation for the electrostatic quantum dots 34 and the construction of four such electrostatic quantum dots setup should be within the experimental reach with the help of a pair of flat gate electrodes with circular intrusion near the center of the dots [35][36][37] .

Three-electron Hamiltonian and electron density
The energy spectrum for three interacting electrons in our system is calculated in a continuum approach using the following Hamiltonian where H 0 is a single electron Hamiltonian in Eq. (1) and the dielectric constant has the value ǫ = 9 for all calcula- tions in the work.The present study calls for an exact treatment of the electron-electron correlation due to the strong electron-electron interaction in phosphorene 29 .We employ the configuration interaction approach 38,39 .In this method, the Hamiltonian H 3e is diagonalized in the basis of three-electron Slater determinants con- structed with the lowest 52 eigenfunctions of the single-electron Hamiltonian H 0 .The Slater determinants are built by filling the three electrons in the single-electron states in all available combinations.The three-electron ground state and excited states are then written as the linear combination of these configurations of electrons.The number of single-electron eigenfunctions spanning the basis of the Slater determinants was chosen on the basis of convergence of the three-electron energy eigenvalues.For 52 single-electron eigenfunctions, the basis of three-electron Slater determinants contains 52 3 = 22100 elements.The Hamiltonian matrix is then numerically diagonalized to obtain three-electron eigenfunctions and corresponding energies.
We investigate the electron distribution in the dots by extracting the electron density from the three electron wave-functions , (2) where is the three-electron ground-state wave function.

Results
We discuss our results in the following order: First, we discuss the parameters of the system for the energy spectrum and the electron density, including the spin ground-state configuration.Next, we calculate a complete phase diagram for the nonferromagnetic and Nagaoka ferromagnetic phases.This gives us insight on the configuration for which the Nagaoka ferromagnetism is the most robust.We then discuss a similar phase diagram but with a shallower potential.We investigate the case where we gradually change the configuration of dots to form a pseudo-1D chain of dots to understand the competition between the ferromagnetic and nonferromagnetic states and the effect of anisotropy on this competition.Lastly, we observe the effect on the Nagaoka phase of a changed potential of one of the quantum dots.

Nagaoka spin polarization
The anisotropy of the effective mass causes the x and y axes of the system to be inequivalent.Therefore, we vary both parameters, µ x and µ y to search for the high-spin configuration in the ground state.The potential and square root of the density of the ground state for different combinations of the parameters are plotted in Fig. 2. The potential depth considered for the following calculations is V d = 125 meV.The densities are calculated using Eq.
( 5) www.nature.com/scientificreports/ (5) for the ground state.For spacings µ x = 9 nm, µ y = 7 nm (Fig. 2a) the electrons are located near the center of four dots, the charge densities of the array are distinctly separated with little or no density overlap (Fig. 2b).Note that even though the dots are far away, the square root density in the center of dots is not circular, but is more spread in the x direction (Fig. 2a), which is a consequence of the mass anisotropy.For the rest of the plots with µ y ≤ 6 nm the potential [see Fig. 2c,e,g] resembles two dumbbell-like shapes due to the electron tunneling in the zigzag direction, for which the assumed distance between the dots is smaller than in the armchair direction.The trace of electron tunneling between the dots can be observed in Fig. 2d for µ x = 8 nm, µ y = 6 nm.As we move the dots closer together in the y direction, we can observe in Fig. 2f that the square root densities of the top-bottom pair of dots now overlap considerably.Finally, for dots arranged in a rectangle with µ x = 6.8 nm, µ y = 5.2 nm, the square root electron density is highest along the edges of the rectangle.The energy spectra for the same parameters are plotted in Fig. 3.For µ x = 9 nm and µ y = 7 nm, the ground state at B = 0 is a four-fold degenerate state with spin eigenvalue S = 3/2 but the spin doublet 1/2 is above with the energy difference of 40 µ eV (Nagaoka gap �E = −40 µeV) .The low-energy spectra are illustrated in Fig. 3a.Although the tunneling coupling between the dots is weak, it is already large enough to promote the spin-polarized quadruplet with the total spin eigenvalue of S = 3/2 and the z component eigenvalues σ z = −3/2, −1/2, 1/2, 3/2 in the ground state.An asymptotic case of large interdot distances corresponds to tun- neling and interaction being negligible so that the quadruplet becomes degenerate with the spin S = 1/2 doublet.
Reducing the spacing between the dots to µ x = 8 nm and µ y = 6 nm, the energy gap between the lowest spin S = 1/2 and spin S = 3/2 states becomes 100 µ eV with the high-spin ground state [see Fig. 3b].
The energy spectrum for µ x = 9 nm and µ y = 5 nm is plotted in Fig. 3c.In this case, the electrons interact weakly in the x-direction and strongly in the y-direction.We have effectively two extended quantum dots with an electron in one of the other two dots.Since the ground state of the two electrons for negligible magnetic field is a singlet state 40 , the top-bottom pair of dots will contain two electrons of opposite spins, resulting in the ground state being a spin S = 1/2 state due to the spin of the solitary electron.The ground state becomes spin-polarized by the Zeeman interaction only at a high magnetic field of about 1.9 T.
The spin singlet state is removed from the ground state when the tunneling in the x direction is enhanced for a reduced distance in x direction to µ x = 6.8 nm with µ y = 5.2 nm we obtain the ground state, which is again spin polarized with a large energy gap of �E = −230 µ eV in Fig. 3d.This is the maximal gap that we obtain for the applied single-dot potential depth.
The energy gap E for various geometry of the quantum dot array is given in a phase diagram that displays the nonferromagnetic and Nagaoka ferromagnetic phases in Fig. 4. The diagram is calculated for a small magnetic field of 1 mT.For the case with µ y ≥ 7 nm, the four dots are located far away from each other, rendering tun- neling negligible, and hence the lowest spin S = 1/2 and spin S = 3/2 states are nearly degenerate, as seen on the right side of the figure.In the case where the dots are located much closer to each other in y direction (left side of the diagram), the ground state tends to be the spin-1/2 state as in Fig. 3c.This is a quantitative result for the strong tunnel coupling forming two double-dot subsystems.The region in between corresponds to the extended ferromagnetic system with a spin-polarized ground state.In this region, the tunneling between neighbouring dots is enough for the 'hole' to move around the four dots.The largest gap is �E = −230 µ eV and is located at µ x = 6.8 nm, µ y = 5.2 nm.
To understand in detail the effects of the interactions depending on the location of the dots, we study the cross sections of the phase diagram to see the energy gap as a function of the parameter µ y for various (µ x − µ y ) in Fig. 5a.The top-most (bluest) line shows the change in the energy gap E with µ y for µ x − µ y = 1 nm.It barely goes below the zero line, indicating a very weak Nagaoka ferromagnetic phase near µ y = 6 nm.As the difference µ x − µ y increases, the next four lines from the top slowly start to go much lower than the zero line, indicating that the system is in a stronger ferromagnetic phase and that it would take more energy to invert one of the spins and break the phase.The lowest line is the one that goes the lowest under the zero line, which is for the difference µ x − µ y = 1.6 nm.The second line from the bottom (greenest) is for the difference of 1.9 nm and is now higher than the difference of 1.6 nm.At this point, the electron density begins to take the shape of two dumbbells with decreasing x-direction overlap.
Fig. 5b shows the energy gap as a function of µ x for different y spacings.The plot has few interesting features that give insight on the system.First, every plot has a peak at low values of µ x slightly below the fixed µ y distance.The initial growth of the energy gap at the left side of the plot is due to an extension of the region accessible for electrons that has a larger influence on the spin-unpolarized state than on the spin-polarized state for which the electrons cannot occupy the same location anyway.When the system is separated to single-electron islands maxima the lines dip to a certain value before tending to a constant value as µ x increases.This constant itself tends to zero for the lines as the parameter µ y becomes large, indicating a vanishing tunnel coupling.In order from topmost (red) line in Fig. 5b, for µ y = 4.5 nm the energy gap is always positive and for this y spacing, the system never attains the ferromagnetic ordering.For other parameters, a ferromagnetic ground state is found for a range of µ x .
The maximal Nagaoka gap is expected to shift on the µ x , µ y plane for a varied potential depth that affects the strength of the tunnel coupling.The phase diagram for the potential depth of V d = 60 meVis plotted in Fig. 6.The results are qualitatively similar to the case of V d = 125 meV presented above.The largest Nagaoka gap �E = −211 µ eV occurs for the parameters µ x = 7.52 nm and µ y = 5.35 nm.Because of the shallower potential, the electrons are less localized within the separate dots.This results in a larger tunnel coupling and the Nagaoka phase is achieved for larger inter-dot distances compared to the case of V d = 125 meV.Note that the parameter µ x changed from 6.8 nm to 7.5 nm, while the parameter µ y changed only about 0.15 nm.This can be attributed to the heavier electron effective mass in y direction and the tunneling energy changing more strongly with the interdot distance.

Shift from the rectangular geometry
In order to investigate the robustness of the spin-polarized ground state we deform the rectangular arrangement of dots so that electrons are trapped in a finite pseudo-1D chain.Starting from the initial state with the largest energy gap, i.e. µ x = 6.8 nm, µ y = 5.2 nm, we shift the location of the top-right dot in both x and y directions separately.The energy gap as a function of the shift �µ is plotted in Fig. 7.The anisotropy in the effective mass makes the nature of Nagaoka transition different in armchair (x) and zigzag (y) direction.For a shift in armchair direction, the spin-1/2 states have much lower energy than the high-spin state and the electrons from a 1D pseudo chain structure as seen in inset (f) of Fig. 7.The electron occupancy of the lower right dot is low for the shift of �µ y = 3.5 nm in the zigzag direction.One of the electrons get fixed in the shifted dot, which is far away from the rest, which minimizes the Coulomb interaction energy.The two remaining electrons occupy the deeper left dumbbell instead of the shallower lower right dot [see inset (a) and (b) of Fig. 7].No tunneling is possible between the rightmost dots and only a trace tunneling in the topmost dots.A similar effect is seen in the x-direction shift, but the geometry of the dots is such that the electrons form a chain with no tunneling between the two upper dots.Ref. 24 showed that when the dot arrays are deformed to form a quantum dot chain, the spin polarization in the ground state is excluded by the Lieb-Mattis theorem 40 , which restricts the ground state solutions of such a 1D chain to low spin values.However, the anisotropy of the masses leads to an asymmetry in the tunneling of electrons between the dots.In the case of a shift in the zigzag direction, the system can  (b) as a function of parameter µ x for various µ y 's.The lines shift towards the negative sides as the system parameters approach the largest Nagaoka gap point, after which the lines in both subfigures approach zero.
be divided into two parts: the left double-dot subsystem and two single dots.As in the case of Fig. 3c, the double dot holds the singlet state of two electrons, and the third electron lowers the Coulomb repulsion by occupying the shifted dot.Tunnel coupling in the y direction decreases rapidly due to the higher effective mass, resulting in a transition occurring at a shift of about ∼ 1.5 nm.In contrast, lifting the Nagaoka ground-state polarization requires a shift of approximately 1.85 nm in the x-direction.

Potential detuning
In addition to moving the dots, it is possible to assess the tolerance of the ferromagnetic state to disorder by changing the potential depth of one of the dots.We vary the potential depth V d of just the top right dot by an amount dV.Fig. 8 shows the Nagaoka gap as a function of the change dV in the range −20 meV to 20 meV.The changes in the electron density as the potential changes are also shown in the insets of Fig. 8.The most evident feature of the plot in Fig. 8 is the asymmetry in the energy gap variation for negative and positive potential changes.More precisely, for the positive change dV the Nagaoka ferromagnetic transitions to the low-spin state at dV ≈ 7.0 meV, while the same transition occurs at a bit smaller change of dV ≈ −5.4 meV for the negative change.The transition to the low-spin ground state for positive detuning is due to the localization of an electron in the detuned dot and the transition for negative detuning is due to the delocalization of an electron from the detuned dot which is excluded from the array by the energy mismatch leaving the three dots with an exact half-filling.When the top right dot is made deeper [insets (c) and (d)] the electron occupancy of this dot becomes larger than 1.When the dot is made shallower [insets (a) and (b)] the dot is emptied.In both cases the conditions for observation of the itinerant ferromagnetism are lifted, and the ground state acquires the low spin.The color scale for the potential and square root of density is shown.Unlike the case in Section "Shift from the rectangular geometry", the anisotropy of the curve is inherent and is not unique to phosphorene.

Figure 1 .
Figure 1.(a) Crystal structure of monolayer phosphorene indicating the zigzag direction (y-axis) and the armchair direction (x-axis).(b) Top view of the phosphorene crystal.

Figure 2 .
Figure2.The potential (left column) and the square root of electron density (right column) with parameters, µ x = 9 nm, µ y = 7 nm (a, b); µ x = 8 nm, µ y = 6 nm (c, d); µ x = 9 nm, µ y = 5 nm (e, f); µ x = 6.8 nm, µ y = 5.2 nm (g, h).For all the cases, the magnetic field is 1 mT and the potential depth V d = 125 meV.The color scale for each plot is to the right of each plot.Subfigures (g, h) indicate the potential and square root of the density for the configuration of dots with the largest Nagaoka gap.

Figure 3 .
Figure 3.The energy spectrum as a function of magnetic field B for parameter sets of µ x = 9 nm (a), µ y = 7 nm, (b) µ x = 8 nm, µ y = 6 nm, (c) µ x = 9 nm, µ y = 5 nm, (d) µ x = 6.8 nm, µ y = 5.2 nm.Potential depth is set V d = 125 meV for all the chosen sets in (a-d).The colour of each line denotes the spin eigenvalue of the energy state, and the scale is shown to right of each plot.The spectra show the shifting of ground states and the changing Nagaoka gap as the configurations of the dots change.Subfigure (d) shows the spectrum for the configuration with the largest Nagaoka gap.

Figure 4 .
Figure 4.The energy gap, E as a function of the parameters µ y and (µ x − µ y ) with potential depth, V d = 125 meV.The energy gap is calculated for every point at the center of the hexagons.The green and grey regions in the diagram indicate the configurations for which the ground state is a low-spin state.The red and yellow regions show the configurations with a spin polarized ground state.The largest energy gap �E = −230 µ eV occurs at µ y = 5.2 nm and µ x − µ y = 1.6 nm and is indicated with the most saturated red color.

Figure 5 .
Figure 5. (a)The Nagaoka gap as a function of the parameters µ y for different values of (µ x − µ y ) . (b) as a function of parameter µ x for various µ y 's.The lines shift towards the negative sides as the system parameters approach the largest Nagaoka gap point, after which the lines in both subfigures approach zero.

Figure 6 .Figure 7 .Figure 8 .
Figure 6.The energy gap, E as a function of the parameters µ y and (µ x − µ y ) with potential depth, V d = 60 meV.The energy gap is calculated for every point at the center of the hexagons.The largest energy gap �E = −211 µ eV occurs at µ y = 5.35 nm and µ x − µ y = 2.21 nm.The diagram is similar as Fig. 4, but the lower potential depth has resulted in a major shift of the largest Nagaoka gap up on the µ x − µ y scale and a minor shift on µ y scale.