Exact Coulomb cutoff technique for supercell calculations in two dimensions

We present a reciprocal space technique for the calculation of the Coulomb integral in two dimensions in systems with reduced periodicity, i.e., finite systems, or systems that are periodic only in one dimension. The technique consists in cutting off the long-range part of the interaction by modifying the expression for the Coulomb operator in reciprocal space. The physical result amounts in an effective screening of the spurious interactions originated by the presence of ghost periodic replicas of the system. This work extends a previous report [C. A. Rozzi et al., Phys. Rev. B 73, 205119 (2006)], where three-dimensional systems were considered. We show that the use of the cutoffs dramatically enhances the accuracy of the calculations for a given supercell size, and it allows to describe two-dimensional systems of reduced periodicity with substantially less computational effort. In particular, we consider semiconductor quantum-dot arrays having potential applications in quantum information technology.


I. INTRODUCTION
The interest in low-dimensional electronic structures has increased steadily during the past few decades. This is mostly due to breakthroughs in semiconductor technology in the 1970s and early 1980s. At present, low-dimensional systems form a significant portion of the whole field of condensed matter physics. Some examples are layered semiconductor devices such as metal-oxide-semiconductor field effect transistors, quantum Hall systems, spintronic devices, and quantum wells, wires, and dots 1 (QDs).
Following the advances in constructing techniques for single two-dimensional (2D) QDs with tunable atom-like properties, it has become possible to couple QDs to form artificial "molecules".
Coupled QDs have significant potential for solid-state quantum computation through, e.g., coherent manipulation of spins. 2 Furthermore, extended lattices or arrays of QDs have been fabricated, 3 which, in addition to the proposed applications in quantum information, 4 show interesting magnetic phase transitions 5,6 which may be exploited in quantum transport and spintronics.
From the theoretical point of view, dealing with systems of arbitrary periodic dimensionality may be complicated. In the simplest case of a fully periodic lattice of elements, periodic boundary conditions are applied at every cell border, and Bloch's theorem describes the discrete-translation invariant form of the orbitals. On the other hand, in all the cases in which the system has reduced periodicity, the use of a supercell with periodic boundary conditions in all the directions becomes problematic. In fact, the response function of a periodic lattice is generally very different from the response of a system with reduced periodicity (such as an isolated system, a chain, a slab, etc.) and the convergence of the fully periodic quantity to the reduced-periodic ones as a function of the supercell size is very slow. A large supercell is a numerical disadvantage, but it is specially necessary to avoid the influence of the periodic images if long range operators are used.
The main issue here is indeed the computation of a long-range operator, i.e., the Hartree (or "Coulomb") potential: which is ubiquitous in Science, and which we study here in the context of 2D electronic structure calculations. In particular, we will exemplify our approach by utilizing density-functional theory (DFT), although the method we propose can also be useful in different fields. In Fourier space, the convolution integral (1) is transformed into a trivial product. This fact adds to the other undoubted advantages of the supercell approach, such as the natural inclusion of the periodic boundary con-ditions, and the existence of very efficient fast Fourier transform algorithms.
The attempt to retain these advantages, i.e., to compute the Hartree integral in reciprocal space, has led to the creation of several cutoff schemes for finite systems, whose main intent is to provide an effective truncated Coulomb interaction such that the system becomes unaware of the existence of its periodic replicas. 7,8,9 More recently, an exact scheme was proposed to achieve a broader goal, 10 namely to truncate the Coulomb interaction in a 3D system in the dimensions along which the system is confined, leaving it long-ranged in the dimensions in which the system is periodic.
In this paper we show that a similar scheme can be drawn in a 2D space, allowing us to correct the spurious supercell effect when treating finite systems and one-dimensional (1D) chains. In particular we focus on the case of a single infinite chain of 2D few-electron QDs, that, in a classical supercell approach, would mistakenly appear as periodic in both directions, while it should be treated as a truly periodic system only along the x direction. We study the magnetic ground state

II. METHOD
We follow closely the procedure described in Ref. [10]. The Hartree integral is the convolution of the charge density n with the Coulomb interaction potential v. In 2D, We consider the charge density to be in a unit cell. This unit cell may then be replicated in both directions to fill all 2D space (the two dimensional periodic case, or "2D/2D"), replicated in only one direction ("2D/1D"), or not replicated at all (finite case, or "2D/0D"). If we move to reciprocal space, the unit cell is always replicated periodically in both directions, which is undesired in both the 2D/1D and 2D/0D cases. The Hartree integral in reciprocal space reduces to the simple multiplication: where G = (G x , G y ) are the reciprocal vectors. The Fourier transform of v(x, y) can be readily This expression implies full 2D periodicity, and therefore it contains spurious terms if the periodicity is reduced. Our aim is to modify the expression of the Coulomb interaction v(G) →ṽ(G) to an effectively truncated interaction:ṽ for some suitable region D , such that it avoids the interaction of the "real" cells with the spurious replicas, while maintaining all interactions between points in the real cells. In order to achieve this goal, it will be necessary to increase the size of the original unit cell.
For finite systems (2D/0D) the cutoff region can be conveniently defined as r < R, for some cutoff radius R. We can then easily perform the Fourier integral in polar coordinates: where J 0 is the Bessel function of the first kind, and 1 F 2 is the generalized hypergeometric function.
Note that the G = 0 case is finite and continuous, and poses no difficulty:ṽ 0D (G = 0) = 2πR. The value of R must be sufficient to contain all possible interactions in the original unit cell where the charge is contained; if we imagine it to be a square of radius L, then R = √ 2L. If we now enlarge this cell (padding the density with zeros) to L ′ = (1 + √ 2)L, the spurious replicas will not interact thanks to the interaction cutoff, and the Hartree integration will be exact.
The 2D/1D case is more subtle. We assume the charge density to be contained in a strip defined by |y| < R/2: the system is a chain of unit cells along the x axis. We define the cutoff region D as |y| < R -and the unit cell is also enlarged in the y direction to |y| < R. It is easy to see that in this manner, the "ghost" replicas do not interact with the original chain, but the real interactions are preserved. The reciprocal space expression for the truncated Coulomb potential is: = 4 where K 0 is the modified Bessel function of the second kind. However, in this case, since lim G x →0 + K 0 (G x y) = +∞, the integral is undefined on the whole line G x = 0.
In the fully periodic (2D/2D) case, we also have a singular point at G = 0 [see Eq. (4)]. If we assume charge neutrality, this singularity poses no problem. In Eq. (3) the term v(G = 0) multiplies n(G = 0), which is zero because this is in fact the charge neutrality condition. In the 2D/1D case, we must estimate how these divergent terms affect Eq. (3), or, more precisely, the back Fourier transform In order to see how the infinities appear, we consider the integral in Eq. (7) for G x = 0, but using first a finite integration domain also in the x direction, −h < x < h. This integral is convergent.
If we perform the integration and retain only the terms that do not vanish in the limit h → ∞ we The first term, which we call v ∞ (0, G y ), diverges as h → ∞. However, it can easily be seen that it can be ignored if we assume charge neutrality. We perform the G y integration in Eq. (9) for Now we have: In the integral |y ′ | < R/2, since the charge is contained in that region. This is the region of interest, and therefore we are also interested in looking at the potential only for |y| < R/2. As a consequence, |y ′ − y| < R, and we can conclude Z dG y v ∞ (0, G y )n(0, G y )e −iG y y = 2π log(2h) This is obviously zero if we assume charge neutrality. Therefore we have prooved that we can safely ignore the diverging terms, and retain only the regular ones. With this in mind, the case of G = 0 is now trivial:ṽ 1D (0, 0) = −4R(log R − 1).

III. RESULTS
To test the cutoff method, we consider 1D QD arrays similar to those in Ref. [6]. Each rectangular unit cell contains two QDs, and each QD has N electrons bound by a Gaussian positive background charge density. The total Gaussian background charge density has the form where r = (x, y), R = (na x , 0) with n = 0, 1, 2, . . ., and r s = 2 is the average density at the center of the QD. Note the use of eff. a.u. throughout. 11 We solve the Kohn-Sham equations within spin-DFT on a 2D grid with and without the cutoff method described above. In the former case, the system is periodic only in x direction, whereas in the latter case it is periodic in both x and y directions. For the exchange and correlation we use the 2D local-spin density approximation (LSDA) with the parametrization of the correlation by Attaccalite et al. 12 This parametrization has been shown to be more accurate than the form of Tanatar and Ceperley 13 in the partially spin-polarized regime. 14 All the numerical calculations are done using the octopus 15 code.
In Fig. 1 we show the Fermi energy E F of a QD chain with three electrons per dot as a function of the lattice constant a y perpendicular to the chain. The lattice constant in the x direction is fixed to a x = 5.06. Using the cutoff method leads to a very fast convergence, so that at a y = 3a x the Fermi energy is converged to six digits. In contrast, without cutoff, i.e., the system being periodic in both directions, a much larger supercell is needed in order to achieve comparable accuracy. This is due to the long-range Coulomb interaction between parallel QD chains. Therefore, for accurate calculations the cutoff scheme is essential in reducing the cell size and thus the computational cost.
Next we focus our attention to the physical effects caused by multiple parallel QD chains in comparison with a single QD chain. In the former case no cutoff is used so that for a y = a x we have parallel chains located at y = na x , n = ±1, ±2, . . .. In the latter case (single chain) we use the cutoff method with a y = 3a x to guarantee a high precision (see Fig. 1).
In Fig. 2 we plot the band structure for a QD chain with three electrons per dot at a x = 5.06 (left panel), corresponding to the system considered in Fig. 1. We find that the presence of parallel chains leads to a rather minor effect on the bands. The qualitative shape of the structure is very similar, and the shift in the bands due to the interchain effects is of the same order of magnitude as the shift in E F . Fig. 2 can be also directly compared with the results of Ref. [6]. Overall, we find an excellent agreement, also in the case of a larger lattice constant a x = 8.55 (right panel). Here the spin degeneracy is lifted due to the exchange effect which leads to a magnetic ground state. 6 Next we consider the interchain effects on phase transitions in QD chains. As noted by Ceperley. 13 However, further studies are required regarding the validity of the LSDA in the highcorrelation regime and/or in the presence of magnetic fields. These would also be ideal systems to test the 2D density functionals recently developed. 16,17

IV. SUMMARY
The reciprocal space is the natural venue for treating periodic systems. If the periodicity is not complete (the system is not periodic in all the space dimensions) or even absent, the computations require the use of a large supercell in order to avoid the spurious interactions due to "ghost" system replicas; this is specially manifested in the computation of the Coulomb, or Hartree, integral, which is a trivial multiplication if we use plane waves. For 3D systems, it has been shown that the use of a screened Coulomb interaction greatly improves the accuracy in the calculation of ground states quantities, and substantially simplifies the evaluation of excited-state properties of reducedperiodicity systems. In this work we have shown that the same ideas can be applied to the case of the 2D electron gas; we have provided the relevant formulae for finite systems in 2D, and for systems that are periodic in only one dimension.
Moreover, the corrective cutoffs are exact and rather straightforward to apply. We have demonstrated their effectiveness by computing band structures and Fermi energies of one-dimensional periodic arrays of quantum dots formed by the confinement of two-dimensional electron gas. We expect that our technique will be of great interest for studying these novel systems of "artificial molecules" and "crystals", although the ubiquity of the problem we have faced will probably make our technique relevant in other fields, too.