Wigner molecules in phosphorene quantum dots

We use the single-band approximation for Hamiltonian describing the electrons of the conduction band of monolayer phosphorene [35, 64]

	$\displaystyle H_{0}$	$\displaystyle=$	$\displaystyle{\left(-i\hbar\frac{\partial}{\partial x}+e{A_{x}}\right)^{2}}/\,{2m_{x}}+{\left(-i\hbar\frac{\partial}{\partial y}+e{A_{y}}\right)^{2}}/\,{2m_{y}}$		(1)
		$\displaystyle+$	$\displaystyle W(x,y)+{g\mu_{B}B\sigma_{z}}/{2}\,,$

where $W(x,y)$ is the confinement potential. In Eq. (1) we use the effective masses $m_{x}=0.17037m_{0}$ for the armchair crystal direction ($x$) and $m_{y}=0.85327m_{0}$ for the zigzag direction ($y$). The values for the masses were determined in Ref. [35] by fitting the confined energy spectra of the continuum single-band Hamiltonian to the results of the tight-binding method. A detailed comparison of the spectra as obtained by the continuum model to the tight-binding ones is given in Ref. [35] for the harmonic oscillator potential and in Ref. [64] for the annular confinement. In Eq. (1) we take the Landé factor $g=g_{0}=2$ after the k$\cdot p$ theory of Ref. [65]. The values of experimentally extracted $g$-factors vary; in particular an increase with respect to $g_{0}$ was reported [63] at low filling factors, which is attributed [65, 63] to strong electron-electron interaction effects in black phosphorus. The electron-electron interaction in this work is treated in an exact manner. The spin Zeeman term leads to the spin polarization of the confined system. The exact value of the magnetic field producing the spin polarization is affected by the adopted $g$-factor value, but no qualitative effect for the Wigner crystallization of the charge density is expected as long as the spin-orbit coupling is absent. The spectral features of Wigner crystallization for $g=0$ are discussed in the Appendix.

We work with a square mesh with a spacing $\Delta x$ in both the $x$ and $y$ directions. The Hamiltonian acting on the wave function $\Psi_{\mu,\eta}=\Psi(x_{\mu},x_{\eta})=\Psi(\mu\Delta x,\eta\Delta x)$ in the finite-difference approach reads

	$\displaystyle H_{0}\Psi_{\mu,\eta}$	$\displaystyle\equiv$	$\displaystyle\frac{\hbar^{2}}{2m_{x}\Delta x^{2}}\left(2\Psi_{\mu,\eta}-C_{y}\Psi_{\mu,\eta-1}-C_{y}^{*}\psi_{\mu,\eta+1}\right)$		(2)
		$\displaystyle+$	$\displaystyle\frac{\hbar^{2}}{2m_{y}\Delta x^{2}}\left(2\Psi_{\mu,\eta}-C_{x}\Psi_{\mu-1,\eta}-C_{x}^{*}\psi_{\mu+1,\eta}\right)$
		$\displaystyle+$	$\displaystyle W_{\mu,\eta}\Psi_{\mu,\eta}+\frac{g\mu_{B}B}{2}\sigma_{z}\,,$

where $C_{x}=\exp(-i\frac{e}{\hbar}\Delta xA_{x})$ and $C_{y}=\exp(-i\frac{e}{\hbar}\Delta xA_{y})$ introduce the Peierls phases [66] for the description of the orbital effects of the perpendicular magnetic field $(0,0,B)$. For calculation of the phase shifts, we use the symmetric gauge ${\bf A}=(A_{x},A_{y},A_{z})=(-By/2,Bx/2,0)$. Hamiltonian (2) is diagonalized in a finite computational box with the infinite quantum well set at the end of the box (see the Appendix).

For evaluation of a realistic confinement potential $W$ we use a simple model with a phosphorene plane embedded in a Al₂O₃ dielectric that fills the area between two parallel electrodes [Fig. 1(a,b)]. A higher (lower) potential energy for electrons is introduced at the top (bottom) electrode. The bottom electrode is grounded and contains a protrusion that approaches the phosphorene layer. As a result, the electrostatic potential within phosphorene forms a cavity that traps the electrons of the conduction band. The model is a variation [67, 68] of a gated GaAs quantum dot of Ref. [69]. Below we use two models: one with a circular protrusion (Fig. 1(a)) and the other with a rectangular one (Fig. 1(b)). The latter is used in the following to study the case close to the 1D confinement. The confinement potential to be used in the Hamiltonian is given by $W(x,y)=-eV(x,y,z_{p})$, with the electrostatic potential $V$ that we evaluate by solving the Laplace equation $-\nabla^{2}V=0$ and $z_{p}$ is the coordinate of the phosphorene layer. For evaluation of the potential we use the finite element method similar to the one applied in Ref. [68] for a charge-neutral phosphorene plane. The confinement potential at the monolayer is plotted in Fig. 1(c) for the circular protrusion of Fig. 1(a) and in Fig. 1(d) for the rectangular protrusion of Fig. 1(b).

The system of $N$-confined electrons is described with the Hamiltonian

$$H_{N}=\sum_{i=1}^{N}H_{0}(i)+\sum_{j>i}^{N}\frac{e^{2}}{4\pi\epsilon_{0}\epsilon\,r_{ij}}\,.$$

(3)

We take the dielectric constant $\epsilon=9$ assuming that the phosphorene is embedded in Al₂O₃.

In the standard configuration-interaction method [55, 56, 57, 58, 59] the $N$-electron Hamiltonian is diagonalized on the basis of Slater determinants constructed with the single-electron Hamiltonian $H_{0}$ eigenstates. Each of the Slater determinants defines a configuration, e.g. a distribution of electrons over the single-electron states. The number of Slater determinants to be used in the calculation is established by a study of the convergence of the energy estimates. Reaching convergence in the present calculation is challenging because of the strong electron-electron interaction in phosphorene [63]. The energy of the ground-state estimated for $N=4$ at $B=0$ in the circular potential of Fig. 1(c) is plotted with the black line in Fig. 2(a) as a function of the number of the lowest-energy single-electron states $\nu$ that span the Slater determinant basis. The Hamiltonian $H_{N}$ commutes with the operator of the $z$ component of the total spin and also with the parity operator due to the point symmetry of the potential. The symmetries allow for a few-fold reduction of the number of basis elements. In Fig. 2 the four-electron ground state at $B=0$ is the spin singlet $S_{z}=0$ of an even spatial parity. Only Slater determinants of these symmetries contribute to the ground-state wave function. For $\nu=60$ the Slater determinants set counts $\binom{60}{4}=$487 635 elements, of which only about 94 500 determinants correspond to $S_{z}=0$ and either even or odd parity. The right vertical axis in Fig. 2(a) shows the number of Slater determinants of the spin-parity symmetry that are compatible with and contribute to the ground state.

The convergence of the CI method using the $H_{0}$ single-electron eigenstates (black line in Fig. 2(a)) is slow. The electron-electron interaction has a pronounced effect on the electron localization, since the electrons in phosphorene are quite heavy as compared to those in e.g. GaAs, and the deformation of the charge density in terms of the single-electron energy is cheap. This in turn results in a high numerical cost of the convergent calculations that require a large number of single-electron states to be included in the convergent basis. Therefore, due to the strong electron-electron interaction, the set of $H_{0}$ eigenstates is not the best starting point for a convergent CI calculation. The literature indicates a number of methods to speed-up the convergence, including the HF+CI method [60, 61, 62] where the basis for the CI method is based on the Hartree-Fock single-electron spin-orbitals. In the HF+CI approach, the mean-field effects of the electron-electron interaction are accounted for already in the single-electron basis and the CI is responsible only for description of the electron-electron correlation effects that evade the mean-field treatment. The HF charge density in the unrestricted version of the method breaks the symmetry of the confinement potential and its restoration is challenging [70, 71] on its own at the CI stage. Convergence speed-up by the choice of the single-electron basis is achieved [72, 73] with the natural orbitals [74, 75] introduced by Löwdin [76].

In this work, we apply a a simple approach that allows for the convergence speed-up by replacing the potential $W$ in $H_{0}$ by another potential that produces single-electron wave functions of the low-energy spectrum that cover a larger area than the ones for the bare potential $W$. For preparation of the single-electron basis we diagonalize the single-electron Hamiltonian $H_{0}^{\prime}$ with the potential

$$W^{\prime}=W+V_{0}\exp(-(x^{2}+y^{2})/d^{2}).$$

(4)

The $H_{0}^{\prime}$, eigenfunctions are used for the Slater determinants to diagonalize the Hamiltonian $H_{N}$. $V_{0}$ and $d$ of Eq. (4) are used as variational parameters in terms of the $N$-electron energy [77]. The results for 4 electrons and the basis of the $H_{0}^{\prime}$ eigenstates for optimized Gaussian parameters given by the red line in Fig. 2(a) exhibit a substantial convergence speed-up with respect to $H_{0}$ basis. In particular, the basis of $H_{0}^{\prime}$ eigenstates with $\nu=36$ and about 12 thousand elements produce a similar ground-state energy estimate as the $H_{0}$ basis with $\nu=52$ and as much as about 53 thousand Slater determinants.

Figure 2(a) contains also the results obtained with eigenfunctions of Hamiltonian $H_{0}^{\prime\prime}$ (green line) using potential $W^{\prime\prime}(x,y)=W(x)/s$, where $s$ is the scaling factor of the bare potential with its variationally optimal value of $s=2.13$. The scaling enlarges the area covered by the low-energy single-electron wave functions in a manner that becomes equivalent in terms of the 4-electron ground-state energy for the one using $W^{\prime}$ potential for $\nu>32$. The low-energy single-electron wave functions for $W$, $W^{\prime}$ and $W^{\prime\prime}$ potentials are given in the Appendix.

Figure 2(c-k) shows the ground-state charge density in the $x>0$, $y>0$ quarter of the QD as obtained for $\nu=24$ (first row of plots), $\nu=40$ (second row of plots) and $\nu=60$ (third row of plots) with the $H_{0}$ (left column), $H_{0}^{\prime}$ (central column) and $H_{0}^{\prime\prime}$ (right column) eigenfunctions. For $\nu=60$ the results are similar for all the three bases. For lower $\nu$ the results for $H_{0}$ (Fig. 2(c,f)) and $H_{0}^{\prime\prime}$ (Fig. 2(e)) the islands appear closer to the origin than in the convergent result.

In order to illustrate the role of the modified potential in the description of the electron-electron interaction, we plotted in Fig. 2(b) the energy overestimate that is obtained once the basis of Slater determinants is reduced by exclusion of all configurations with more than two electrons [60] above the $32$-nd single-electron energy level. For the Hamiltonian $H_{0}$ the cost of the limited basis is much larger than for $H_{0}^{\prime}$ and grows fast with $\nu$. For $H_{0}^{\prime}$ the overestimate is much lower. The single-electron effects due to $W^{\prime}$ potential are included in the basis, and the double excitations that stay in the basis cover most of the electron-electron correlation effects. A similar result is obtained in the HF+CI method [61, 62].

For analysis of the electron localization, we extract the charge density and the pair correlation function from the $N$-electron wave function $\Psi$. The charge density is obtained as

$$\rho({\bf r})=\langle\Psi|\sum_{i=1}^{N}\delta({\bf r}_{i}-{\bf r})|\Psi\rangle\,.$$

(5)

The pair correlation function extracts the relative localization of the electrons with one of the carrier positions fixed

$$\rho_{12}({\bf r},{\bf r}_{f})=\langle\Psi|\sum_{i,j=1}^{N}\delta({\bf r}_{i}-{\bf r})\delta({\bf r}_{j}-{\bf r}_{f})|\Psi\rangle\,.$$

(6)

In the following we fix the position ${\bf r}_{f}$ of one of the electrons near the local charge density maximum for a discussion of $\rho_{12}$ plots. The spin density can be calculated as

$$\rho^{\alpha_{i}}({\bf r})=\langle\Psi|\sum_{i=1}^{N}\delta({\bf r}_{i}-{\bf r})|\alpha_{i}\rangle\langle\alpha_{i}||\Psi\rangle\,.$$

(7)

where the projection on the spin eigenstates uses $|\alpha_{i}\rangle$ that stands for the single-electron spin eigenstate for $i$-th electron. Similarly, the relative localization including the spin configuration of the electron pair can be included in the $\rho_{12}$ pair-correlation function. In particular, opposite spin distribution can be obtained using spin projections

	$\displaystyle\rho_{12}^{\uparrow\downarrow}({\bf r},{\bf r}_{f})$	$\displaystyle=$	$\displaystyle\langle\Psi|\sum_{i,j=1}^{N}\delta({\bf r}_{i}-{\bf r})\delta({\bf r}_{j}-{\bf r}_{f})\times$		(8)
		$\displaystyle\times$	$\displaystyle(|\alpha_{i}\beta_{f}\rangle\langle\alpha_{i}\beta_{f}|+|\beta_{i}\alpha_{f}\rangle\langle\beta_{i}\alpha_{f}|)\Psi\rangle\,,$

with the spin eigenstates $|\alpha\rangle\neq|\beta\rangle$ .

We discuss first the relatively simple case for $N\leq 4$ (Section III.A.1) where in the low energy states the Wigner molecule formation in the laboratory frame is present $(N=2,N=4)$ or absent $N=3$. Section III.A.2 covers the case for $N=5$ where states of both types are present in the low-energy part of the spectrum.

We discuss two perpendicular orientations of the elongated quantum dot to study the interplay of the potential and effective mass anisotropies.

Figure 2 presents the convergence of the CI method the single-electron Hamiltonians $H_{0}$, $H_{0}^{\prime}$ and $H_{0}^{\prime\prime}$ using the bare QD potential $W$, the potential with the central repulsive Gaussian ($W^{\prime}$) and the reduced potential $W^{\prime\prime}$, respectively. The low-energy single-electron wave functions are plotted in Fig. A12 for the three Hamiltonians starting from the ground-state in the first row, and the subsequent excited states presented in the lower panels of Fig. A12. The seven lowest-energy states have the same character for all Hamiltonians only with states for $H_{0}^{\prime}$ (central column of plots) and $H_{0}^{\prime\prime}$ (right column) covering a larger area than the ones for $H_{0}$, which turns out to speed-up the convergence of the CI method for the system of size increased by the strong electron-electron interaction.

The CI method uses the single-electron wave functions that are obtained with the finite difference approach in a finite computational box with the boundary conditions of vanishing wave function at the edges of the box. The edges of the box form effectively an infinite quantum well that has to be chosen large enough to contain the few-electron system without perturbation to the low-energy states. The influence of the size of the box on the results is given in Fig. A13 for the circular quantum dot. For the circular quantum dot we use the square grid of points of the spacing of $\Delta x=R/nx$ where $R$ is the radius of the computational box and we take $n_{x}=111$. The red (black) line in Fig. A13 shows the ground-state energy for $N=2$ ($N=5$). The growth of the energy for small $R$ is due to the finite-size effect with the quantum well ground state changing as $1/R^{2}$. For calculations for $N=2$ the radius of $R=50$ nm is large enough while for $N=5$ is has to be taken as large as $R=80$ nm.

A striking difference in the energy spectra states with or without the Wigner molecule in the ground state is revealed for the circular potential once the spin Zeeman interaction is removed. This is illustrated in Fig. A14 which reproduces the data from Fig. 3 for $g=0$. The systems with Wigner form of the charge density, i.e. the single-electron islands in the charge density), $N=2$ (Fig. A14(a)) and $N=4$ (Fig. A14(c)) contain a nearly degenerate ground state with energy levels of different parity that interlace in increasing $B$. The separation of the electron charges in the separate, weakly coupled, maxima produces a small energy difference due to the parity, which is similar to the nearly degenerate ground state for identical, weakly coupled quantum dots. The single electron islands are formed by a relatively strong electron-electron interaction and the weakness of the electron tunnelling between the islands can be deduced from the charge density plots of Fig. 4(a,b) and Fig. 4(e,f) for $N=2$ and $N=4$. In both cases, the bunch of energy levels of the ground state is separated by a distinct energy gap from the excited states.

For $N=3$ [Fig. 8(b)] the parity of the ground state changes in growing $B$ as for even $N$ but the spacing between the even and odd parity energy levels near the ground state is much larger than for even $N$ and no energy gap is observed between the nearly degenerate ground state and the rest of the spectrum.