Stability of the Dirac cone in artificial graphene formed in quantum wells: a computational many-electron study

We carry out a comprehensive computational study on the stability of the Dirac cone in artificial graphene realized in nanopatterned quantum wells. Our real-space approach allows us to vary the size, shape, and positioning of the quantum dots in the hexagonal lattice. We compare the (noninteracting) single-particle calculations to density-functional studies within both local-density approximation and meta-generalized-gradient approximation. Furthermore, the density-functional results are compared against numerically precise path-integral quantum Monte Carlo calculations. As a whole, our results indicate high stability of the Dirac bands against external parameters, which is reassuring for further experimental investigations.


Introduction
Artificial graphene (AG) is a generic term for physical systems that employ a man-made honeycomb lattice for charge carriers [1]. AG has been realized in several systems such as nanopatterned two-dimensional (2D) electron gas (2DEG) in a semiconductor heterostructure, [2][3][4][5] copper surfaces covered by molecules [6], and in trapped cold atoms in an optical lattice [7,8]. Although these systems lack many important properties of graphene such as the plasticity, transparency, and strength, the geometry of the lattice itself can be controlled almost at will. This enables tunable band structures that can be used to design new Dirac materials.
Following the pioneering works of Park and Louie [2] and Gibertini et al [3] AG realized on high-mobility 2DEG has attracted many further investigations. The quantum dot lattice in a GaAs heterostructure has a tunable lattice constant from tens to hundreds of nanometers. Hence, the electron-electron (e-e) interactions in the lattice have an important role that leads to collective excitations according to the Mott-Hubbard model [4]. In a previous study of some of the present authors [5] it was found that the Dirac point is relatively stable against e-e interactions. However, the general stability of graphene-like properties in nanopatterned AG towards, e.g., lattice distortion and the shape of the QDs is still unknown. This is a critical issue for the usability of semiconductor AG in applications. Moreover, strong e-e interactions may enhance the sensitivity to irregularities in the QD lattice.
In this paper we show that the band structure of AG realized in nanopatterned 2DEG is very stable with respect to the shape and size of the QD potential and to external lattice distortion-at least within experimental variations. This is surprising in view of, e.g., recent studies on molecular graphene showing rather complex band structures [8]. In conjunction with these findings we bring the theoretical modeling of AG systems into the next level. In particular, we employ realistic real-space modeling with density-functional theory [9] (DFT) applied within both the 2D local-density approximation [10] (LDA) and the 2D meta-generalized-gradient approximation [11] (mGGA). These computations are further supplemented by 2D path-integral Monte Carlo (PIMC) simulations that provide us with numerically exact densities for comparison. The mGGA results confirm the existence of the Dirac cone in an ideal AG previously studied with tight-binding calculations [3] and with the LDA [5]. The mGGA yields a larger gap between the lowest and excited bands in accordance with manybody calculations on other graphene-like systems [12][13][14]. This improvement is due to the better description of the electronic exchange in our mGGA, particularly the correct asymptotic behavior and the correct one-particle limit of the exhange potential [11]. However, in the parameter range corresponding to the transition between a Dirac system and a metallic state, we report and analyze cases that are problematic for the mGGA. In this respect, there is still the need for an all-around DFT functional for correlated 2D systems. Nevertheless, the present study confirms the Dirac properties of AG through accurate many-electron calculations and suggests that those properties are robust against several experimental parameters.

Numerical framework
We consider electrons confined in circularly symmetric GaAs/AlGaAs QDs that are arranged in a hexagonal (honeycomb) lattice. We use the effective mass approximation with the material parameters * = m m 0.067 0 for the effective mass and   = 12.4 0 for the dielectric constant. Each QD is occupied with one electron, the potential depth is = V 0.6 0 meV, the lattice constant (distance between the QDs) is a=150 nm, and the QD radius R and the shape of the potential are varied-see figure 1 and text below. We point out that the lattice constant a determines the relative strength of e-e interactions as the Coulomb energy scales as µ r 1 in contrast with the µ r 1 2 scaling of the kinetic energy. The value a=150 nm chosen here agrees with the experiments of Gibertini et al [3] and corresponds to a relatively strongly interacting system. It can be expected that forthcoming experiments aim at reducing the lattice constant [15].
Our model is strictly 2D so that the degrees of freedom in the direction perpendicular to the 2D plane has been omitted. However, the e-e interaction is of the Coulomb form, so that we mimic an experimental situation of a quasi-2D system. We use a rectangular 2D grid with periodic boundary conditions. The rectangular supercell includes four QDs, and the relevant path Γ-M-K-Γ can be extracted by the 'unfolded' path as described and visualized in reference [5].
We perform four types of calculations: (i) independent (noninteracting) particles (IPs), (ii) DFT with 2D LDA, (iii) DFT with 2D mGGA, and (iv) PIMC. In the LDA we employ the parametrized form for the correlation by Attaccalite et al [10]. In our mGGA [11] we use a 2D version of a Becke-Johnson-type model [16,17] for the exchange potential that produces the set of correct properties such as the asymptotic behavior and the oneparticle limit. For the correlation in our mGGA we use the form of the LDA. Both the IP and DFT calculations are carried out with the Octopus code, [18][19][20][21] where the convergence parameters are mesh spacing and number of k-points in the irreducible Brillouin zone.
Our 2D PIMC simulations require a larger supercell than the DFT calculations due to minimization of the finite-size effects: we see that a supercell of eight QDs is already large enough for good accuracy. The fermion sign problem is dealt with the restricted (or fixed-node) path-integral method using free particle nodes [22]. Our simulations at T=0.1 K together with the Trotter number M=3200 [23][24][25] lead to an accurate description of the ground-state electron density.
The real-space grid spacing is 2.94 nm for IP and LDA, 1.96 nm for mGGA. In the k-space we use a rectangular uniform grid in the irreducible Brillouin zone (rectangle Γ-X-S-Y) with 9×25 k-points for the IP and LDA calculations and with 14×40 k-points for the mGGA calculations. The convergence threshold for the Kohn-Sham calculations is an error of the eigenvalues lower than´-1.0 10 3 meV.

Shape of the quantum dots
First we examine the effect of the shape of the QD confinement on the band structure. The confining potential is given by where α determines the softness of the potential. We compare the Gaussian case with a = 2 to the limit of a hard-wall potential with a  ¥. The dot radius is set to = R 52.5 nm = a 0.35 , and the other parameter values are given in section 2.
The potential profiles and the resulting band structures are shown in figure 1. The (red) dashed lines, (blue) thin solid lines, and (black) thick solid lines show the IP, LDA, and mGGA results, respectively. We first notice that the band structures of the hard-wall (a = ¥) and soft-wall (a = 2) AG are astonishingly similar. Both cases show a distinctive Dirac point, regardless of the type of calculation (IP, LDA, or mGGA). In the hard-wall case, the IP and LDA results are consistent with previous results in [5]. In the soft-wall case, the higher bands are slightly pushed to lower energies. This is expected from the fact that the Gaussian potential spreads out on the edge of the QD, thus increasing the density of states (DOSs).
It can be deduced from figure 1 that with the chosen parameters the e-e interactions play a rather minor role in ground-state properties. Interestingly, LDA and mGGA have opposite effects in comparison with the noninteracting IP result. In the mGGA-expecting to represent here the most accurate calculation-the higher band is shifted up in energy. A qualitative similar correction to the LDA has been found in GW calculations for the band structure of MoS 2 [12] and graphane [13,14]. Nevertheless, as a whole our examination shows that the actual shape of the confining potential does not play a major role for the ground-state properties of AG. On the basis of the vast literature for GaAs QDs it can be expected that the confinement at the bottom of the QD is close to harmonic, [26] which is the case also for a Gaussian model considered here.

Size of the quantum dots
Next we set a = 8 in equation (1) corresponding to intermediate softness of the QD confinement in comparison with the two potentials in the previous section ( figure 2). The lattice constant a is also kept constant, but the radius R of the potential is changed. We focus here on the regime, where the QDs become too small for confining electrons. The resulting delocalization of the electronic density is then expected to lead to metallic behavior.
Starting from the IP picture, the left panel of figure 2 shows that the system becomes metallic at = ¼ R a 0.1 0.2. Due to the noninteracting picture, this result applies to any lattice constant a. In the interacting picture, however, a determines the relative strength of e-e interactions (see section 2). Thus, the transition point from a graphene-like band structure to a metal in terms of R/a depends on a.
According to the two rightmost panels of figure 2 there is a major difference between LDA and mGGA for the metallic transition: in the former case this occurs atR a 0.2, whereas the mGGA does not show any transition in the examined range. At small values for R/a the LDA density is spread rather homogeneously in the system, whereas the mGGA preserves the density localization and thus the Dirac cone.
To assess the accuracy of LDA or mGGA, especially at small values of R/a where the differences are considerable, we benchmark the results against PIMC calculations. Figure 3 shows the electron densities with = R a 0.2 (a) and = R a 0.4 (b), while the lattice constant is kept at a=150 nm. Here we show a one-dimensional cut of the full 2D density for better visibility of the differences; the inset in figure 3(a) shows the 1D cut along which the density is plotted. The trends are opposite in the two cases: with small QD radius the PIMC density (thick solid line) is rather homogeneously distributed in the system and only minor localization at the dots is visible. In contrast, both LDA (dashed-dotted line) and especially mGGA (thin solid line) show much stronger localization. In the case of = R a 0.4, however, PIMC shows the strongest localization and mGGA is close, whereas in the LDA the localization is slightly weaker.
On the basis of earlier investigations [11] the 2D mGGA used here is expected to be more accurate than the LDA. As shown in figure 3 this assumption holds only in the strong-localization regime with relative large values of R/a (0.3), which in fact corresponds to the experimental situation in [3]. The discrepancy at small QD radii is due to the exchange potential v x : the overall mGGA exchange potential is deeper and broader than that of the LDA, leading to pronounced localization of the density. Despite the fact that the v x in our mGGA is relatively close to the exact exchange potential with, e.g., ther 1 asymptotic behavior, we assume that the missing correlation beyond the LDA is problematic in this regime. This calls for further examinations and developments of mutually compatible 2D exchange and correlation functionals within DFT.  = R a 0.4, while a=150 nm. The smaller radius pushes the bands to higher energies, which eventually leads to delocalization. This is pronounced in the PIMC result in (a).

Displacement of quantum dots
Finally, we examine the effect of QD lattice distortion on the DOSs. We apply random perturbations to the QD positions (displacements) in the lattice within a magnitude of 0%-2% in the supercell. The angles of the displacements are randomly set within 0°-360°. Figure 4(a) shows the IP, LDA, and mGGA results for the DOS calculated with a Lorentzian broadening of g = 0.006 meV. The results with perturbation (gray lines) display a sum of 25 calculations with random sets of displacements. In all cases, the results of the perturbed calculations are very close to the unperturbed situation: the maximum shifts of the order of 10 −3 meV, and we always find a distinctive 'Dirac valley' surrounded by two peaks in the DOS, where the one with a higher energy has a longer tail [3,5,27,28]. This demonstrates the high stability of AG against displacements in the QD lattice, which is encouraging for further experimental studies exploiting Dirac physics in nanopatterned AG. In molecular graphene [6,8]-representing an alternative experimental approach-this challenge of positioning is not present due to the precise adsorption on an atomic level.
In the above examination, the displacements are also periodic across the supercells, each consisting of four QDs. This raises the question whether the apparent stability of the DOS in figure 4(a) results from the periodic calculation. Therefore, in figure 4(b) we compare the mGGA results for the DOS using the original supercell of four QDs against using doubly sized supercells of eight QDs. The perturbations to all the QDs are random. The agreement between the results confirm that the periodic displacements do not considerably affect the examination of the stability.

Conclusions
In summary, we have computationally examined the stability of the graphene-like band structure and especially the Dirac cone in AG formed in the 2D electron gas in a nanopatterned quantum well. We have employed several computational methods within a realistic real-space model of the system: the IP approximation, DFT within the 2D LDA, the 2D mGGA, and a PIMC method. We have used the experimental parameter values with a large lattice constant (150 nm), so that the e-e interactions are relatively strong.
We have found that the Dirac cone is very stable against the shape and size of the confining potential for the quantum dots, as well as against variations in the lattice geometry. In the mGGA the higher bands are shifted down in energy in accordance with previous studies comparing the LDA and more accurate calculations. Generally, the softening of the potential has very little effect, but if the radius of the quantum dots is considerably decreased we find a transition to a metallic state as the electron density is spread outside the dots. In this regime, mGGA shows over-localization in comparison with the PIMC, although at larger radii mGGA is more accurate than the LDA.
Regarding the lattice variation, we have found that the displacement of quantum dots in the lattice has only a minor effect on the band structure. Therefore, it might be expected that there is a rather large error threshold in nanopatterining experiments to observe and exploit Dirac physics in AG. However, further examinations are needed to systematically study the effects of the fabrication and actual nanopatterining-including size variations of the quantum dots and also smaller lattice constants that are becoming realizable.