Density functional theory for homogeneous two dimensional materials with magnetic fields

This paper studies DFT models for homogeneous 2D materials in 3D space, under a constant perpendicular magnetic field. We show how to reduce the three--dimensional energy functional to a one--dimensional one, similarly as in our previous work. This is done by minimizing over states invariant under magnetic translations and that commute with the Landau operator. In the reduced model, the Pauli principle no longer appears. It is replaced by a penalization term in the energy.


Introduction
The analysis of quantum properties of two dimensional materials is an active research area in physics and material science.Some 2D materials such as graphene or phosphorene exhibits many interesting physical properties [16,3,10,11] which has many applications such as High Electron Mobility Transistors [12].Some of these properties are not yet fully understood.This is the main motivation to revisit Density Functional Theory (DFT) when applied to quantum two dimensional systems (see [1,9] for previous works).
As in our previous work [9], we study homogeneous two-dimensional slabs, when embedded in three dimensional space, but this time, we include a constant perpendicular magnetic field.We consider a charge distribution µ which is equidistributed in the first two dimensions: µ(x 1 , x 2 , x 3 ) = µ(x 3 ), and with a constant perpendicular magnetic field B = be 3 .
One key result of our previous work was an inequality for the kinetic energy per unit surface for translationally invariant states.Let us quickly summarize the result.Let denote the set of one-body density matrices which commute with all R 2 translations.Here, S 1 (L 2 (R 3 )) stands for locally trace class self-adjoint operators with finite trace per unit surface Tr(γ) < ∞ (see Section 2.3).For R ∈ R 2 ⊂ R 3 , we have denoted by τ R f (x 1 , x 2 , x 3 ) := f (x 1 − R 1 , x 2 − R 2 , x 3 ) the usual translation along the first two dimensions.Let us also introduce the set of reduced states where S 1 (L 2 (R)) refers to the space of trace class self-adjoint operators on L 2 (R).We have proved in [9] that, for any (representable) density ρ = ρ(x 3 ) depending only on the third variable, we have (1.1) inf where ∆ d denotes the Laplacian operator in d-space dimension.The energy appearing in the right hand side leads to one-dimensional reduced models for homogeneous semi-infinite slabs in the context of DFT.One typically obtains a minimization problem of the form where D 1 (•) is the one-dimensional Coulomb interaction energy (see [9, Section 3.1] for a discussion about this term), and E xc some exchange-correlation energy per unit surface.Note that there is no Pauli principle for the operator G; it has been replaced by the penalization term πTr(G 2 ) in the energy, which prevents G from having large eigenvalues.We refer to [9] for details, where we also studied the reduced Hartree-Fock model, which corresponds to the case E xc = 0.
The scope of this paper is to apply a similar reduction when taking into account magnetic effects.Without considering the spin (we refer to Section 4 for the case with spin), the kinetic energy per unit surface is of the form 1 2 Tr (−i∇ 3 + A) 2 γ , where A = b(0, x 1 , 0) is a vector potential so that curl A = B = be 3 .We chose a gauge which is not symmetric, but which will simplify some computations.The Laplacian operator −∆ 3 has been replaced by the Landau operator In analogy with [9], we only consider states commuting with the so-called magnetic translations m R and with the Landau operator.We refer to Section 2 for the definition of these operators, and for the justification of this choice.We denote the set of such states by (1.3) In Theorem 2.7, we prove that any γ ∈ P A has a simple decomposition, in terms of the different projectors on the Landau levels.Using this decomposition, we prove in Theorem 3.1 that, in the magnetic case, we have an equality similar to (1.1), which reads with the penalization term F defined by where {x} denotes the fractional part of x.The function F is studied in Proposition 3.2.It is a piece-wise linear function, reflecting the contribution of the different Landau levels.The function F is not new, and already appears in the context of the two-dimensional Thomas Fermi (TF) theory under constant magnetic fields (see [14] and related references).The (spinless) TF kinetic energy takes the form For b = 0, we recover the usual two-dimensional TF kinetic energy ´πρ 2 .It is different from the three-dimensional TF kinetic energy of a gas under a constant magnetic field, which has been derived and studied in [6,19], and is obtained by assuming that the electron density is constant, hence also invariant under the third-direction translation.
Equation (1.4) allows the reduction of DFT models for two-dimensional homogeneous slabs under constant magnetic field.In fact, one obtains a one-dimensional problem of the form (compare with (1.2)) Note that the exchange-correlation function may depend on the external magnetic field b (see [6, Eqn.(4.1)-(4. 2)] for the expression of the exchange energy for the Landau gas).In Section 3.3, we study the corresponding reduced Hartree-Fock model, where This article is structured as follows.In Section 2, we start by introducing the Landau operator and studying its spectral decomposition, then we define magnetic translations {m R } R∈R 2 with some of their properties, and we characterize the states in P A .Using this characterization, we explain in section 3 how to reduce the kinetic energy, we give some properties of the penalization term F and we study the corresponding reduced Hartree-Fock model.Finally, we show in Section 4 how to extend our results to systems with spin.
Acknowledgments.The research leading to these results has received funding from OCP grant AS70 "Towards phosphorene based materials and devices".

States commuting with magnetic translations
In this section, we prove that states γ ∈ P A have a particular structure., where H n (x) = (−1) n e x 2 d n dx n e −x 2 refers to the n-th Hermite polynomial and In particular, the spectral projector Its density ρ Pn (x) := P n (x; x) = b 2π is constant and independent of n.Proof.First, we remark that L A 2 commutes with all translations in the x 2direction.We introduce the Fourier transform with respect to the x 2 -variable and its inverse We have The operator h b,k is a translation of the harmonic oscillator h with ε n = b(2n + 1) and ϕ n as defined in (2.1).We deduce that the spectral decomposition of h b,k is which proves (2.2).Using (2.7), we see that W is an eigenfunction of L A 2 , corresponding to eigenvalue ε n , if and only if it is of the form where g ∈ L 2 (R).Thus To compute the kernel P n (x, y) of the projector on E n , we use the Moyal identity, that we recall here.
Proof.We first prove the result for f 1 , g 1 , f 2 , g 2 ∈ C ∞ 0 (R) and conclude by density.Parseval identity gives In particular, if Hence Using that m |ψ m ψ m | = I L 2 , we obtain P n (x; y) = ϕ n,x , ϕ n,y , which, given that ϕ n is real-valued, gives (2.5).The density of P n is thus 2.2.Magnetic translations.The Landau operator does not commute with the usual translations, however it commutes with the magnetic translations, that we define now.We write The operators p A,1 and p A,2 do not commute, and do not commute with L A 2 .Actually, we have

However, introducing the dual momentum operators
) .Note that we have added a phase factor in order to match the usual convention.Using the Baker-Campbell-Haussdorf formula and the fact that [ p A,1 , p A,2 ] = ib commutes with all operators, we obtain the explicit expression where τ R f (x) := f (x − R) is the usual translation operator.By construction, the magnetic translations commute with L A 2 and P n , but they do not commute among them.Actually, we have An important feature of magnetic translations is that they form an irreducible family on each eigenspace E n , in the sense of [2,Definition 2.3.7].
Proposition 2.4.The set of magnetic translation operators (m R ) R is an irreducible family of operators on each E n , in the sense that Proof.Assume otherwise, and let Ψ so that Then there is Using the Moyal identity, and the fact that Applying the inverse Fourier transform to k → f (k)g(k − R 1 ) = 0 shows that f (k)g(k − R 1 ) = 0 a.e. for all R 1 ∈ R. Squaring and integrating in R 1 gives f = 0, a contradiction.

Diagonalisation of states commuting with magnetic translations.
In what follows, we are interested in one-body density matrices which commute with all magnetic translations.In the case without magnetic field, if a state commutes with all usual translations τ R , then it commutes with the Laplacian operator.In the magnetic case, there are operators which commute with all magnetic translations (m R ) R∈R 2 , but which do not commute with the Landau operator (we give an example of such an operator in Remark 2.6 below).So, we rather consider one-body density matrices which commute with m R , and with the Landau operator.It turns out that such operators have an explicit and simple characterization.
We first enunciate our result in dimension two before turning to the three dimensional case.
If η is a locally trace class operator, then its density is constant Proof.Since η commutes with L A 2 , it commutes with any spectral projector P n , hence leaves invariant E n = Ran(P n ), for all n ∈ N. The operator η n := P n ηP n ∈ S(E n ) commutes with all magnetic translations.Since the family {m R } R is irreducible, it implies that η n is proportional to the identity on E n , hence is of the form η n = λ n I En .This is a kind of Schur's Lemma, see [ with all magnetic translations.However, we have, using the decomposition (2.7) that This happens if and only if ζ is an eigenstate of the harmonic oscillator, that is ζ = ϕ n for some n ∈ N 0 .
The three-dimensional analogue of the previous Proposition reads as follows.
Theorem 2.7.Let γ be a bounded operator on L 2 (R 3 ) commuting with all magnetic translations m R and with L A 2 ⊗ I.Then, there exists a family (γ n ) n∈N of bounded operators on L 2 (R), with γ n ≤ γ , and so that If γ is a locally trace class operator, then its density depends only on x 3 Proof.Let us consider two fixed test functions φ, ψ ∈ L 2 (R), and define the operator The conditions on γ imply that η φ,ψ is a bounded self-adjoint operator that commutes with all m R and with L A .Thus, using Proposition 2.5, η φ,ψ can be decomposed as Since η is a bounded operator, we have for any normalized We deduce that the map (φ, ψ) → λ n (φ, ψ) is sesquilinear and bounded, with bound smaller than η op .The result then follows by taking γ n the bounded self-adjoint operator on L 2 (R) defined by φ, γ n ψ := λ n (φ, ψ).
Finally, for η of the form (2.11), we obtain, using ρ Pn = b 2π , For an operator η of the form (2.11), the trace per unit-surface, defined as the limit takes the simpler form where we have used that the density of any Landau level is ρ Pn = b 2π .

Reduction of the kinetic energy, and applications
We now exploit the particular structure of states γ ∈ P A to deduce their kinetic energy.
3.1.Reduction of the kinetic energy.Recall that the set P A has been defined in (1.3) as the set of one-body density matrices γ, satisfying the Pauli principle 0 ≤ γ ≤ 1, which commute with all magnetic translations, and with the Landau operator L A 2 ⊗ I.We also recall that the set of reduced states G is defined by The main result of this section is the following.Theorem 3.1.For any γ ∈ P A , there is an operator G ∈ G satisfying ρ G = ρ γ and where Conversely, for any G ∈ G, there is γ ∈ P A so that ρ γ = ρ G , and for which there is equality in (3.1).In particular, for any (representable) density ρ, Proof.According to Theorem (2.7), any γ ∈ P A can be decomposed as For a state of the form (3.4), we define the operator G γ ∈ S(L 2 (R)) by Since γ n ≥ 0, we have G γ ≥ 0 as well, so G γ ∈ G. Also, since ρ Pn (x) = b 2π , we deduce that ρ G = ρ γ .
Recalling that , and using Proposition 2.1, we obtain that The first term cannot be expressed directly as a function of G, but we have an inequality for this term.Since G is a positive operator with finite trace, it is compact, and admits a spectral decomposition of the form G = j g j |ψ j ψ j | with g j > 0 and j g j < ∞.Evaluating the trace of γ n in the {ψ j } basis, and changing the order of the sums (all terms are positive), we obtain Since 0 ≤ γ n ≤ 1, the quantity m j (n) := ψ j , γ n ψ j satisfies 0 ≤ m j (n) ≤ 1.
In addition, we have b 2π n m j (n) = b 2π ψ j , n γ n ψ j = ψ j , Gψ j = g j .So we have the inequality (3.5) Since the ε n are ranked in increasing order, we can apply the bathtub principle [13,Theorem 1.4].The optimal m for the above minimization is given by . We now calculate the infimum in the RHS of (3.5) using the explicit formula of the optimal function m * .Recalling that ε n = b(2n + 1) and denoting by with the function F defined in (3.2).Summing in j and gathering the terms gives the inequality which proves the first part of the Theorem.Conversely, given G = j g j |ψ j ψ j | ∈ G, we consider the state and m * j defined as in (3.6).The operator γ * belongs to P A , satisfies G γ * = G, and gives an equality in (3.1).

Some properties of the function F. Let us collect some useful properties of the function F . A plot of F is displayed in
with equality in the left for 2πg b ∈ N, and equality in the right for 2πg b ∈ N+ 1 2 , and F (b, g) → πg 2 as b → 0. For any b ≥ 0, the map g → F (b, g) − πg 2 is b 2π periodic and the map g ∈ R + → F (b, g) is piece-wise linear, increasing and convex.Finally, for all 0 ≤ g < b 2π , we have F (b, g) = 1 2 bg Proof.The first part is straightforward from the definition (3.2).To see that it is convex, piece-wise linear and increasing, we use the alternative form where we have denoted by x := 2πg b .When 0 ≤ g < b 2π , which corresponds to 0 ≤ x < 1, F (b, g) = 1  2 bg.Remark 3.3.The left inequality of (3.7) implies Tr(F (b, G)) ≥ πTr(G 2 ), hence, together with (1.1), that In particular, the kinetic energy is higher with the magnetic field.This is a kind of diamagnetic inequality for 2D materials.
Remark 3.4.The fact that F (b, g) → πg 2 as b → 0, for all g ∈ R + , means that, our reduction approach in this manuscript is coherent with the one without magnetic field already treated in [9].
Remark 3.5.Splitting F into F (b, G) = πg 2 + F (b, g), we see that the effect of adding a magnetic field B = (0, 0, b) is a periodic perturbation of the energy with no magnetic field.

3.3.
Reduced DFT models.In the context of DFT, the previous result suggests modelling the electronic state in a homogeneous slab of charge distribution µ(x) = µ(x 3 ) under a constant magnetic field B = b(0, 0, x 3 ) by a reduced state G ∈ G whose energy per unit surface is given by (3.9) Here, E xc models is an exchange-correlation energy per unit surface, and D 1 is the one-dimensional Hartree term.This last term has been extensively studied in our previous work [9, Section 3.1], and is defined as follows.For We have proved in [9,Proposition 3.3] that the elements f ∈ C have null integral ´f = 0, that the map C f → D 1 (f ) is convex, and that with the mean-field potential The function Φ f is continuous, and is the (unique) solution to In practice, we restrict the minimization problem to the G so that ρ G −µ ∈ C.This implies in particular the neutrality condition Tr(G) = ´ρG = ´µ = ν.
On the other hand, if G is a trace class operator, then ρ G ∈ L 1 (we assume that µ ∈ L 1 as well), and if Note that there is no Pauli principle on G for admissible states G.It has been replaced by a penalization term +F (b, G) in the energy.
Remark 3.6.The energy (3.9) is obtained when minimizing a 3-dimensional DFT model over transitionally invariant states.In particular, this model does not include possible spatial symmetry breaking along the first 2 variables.Such phenomena are known to exist in two-dimensional electron gas under magnetic field due to the de Haas-van Alphen effect [5].In some real-life systems e.g.Br 2 , magnetic domains form, sometimes called Condon domains [4,15].Our simple model is unable to capture these effects.

3.4.
The reduced Hartree-Fock case.Let us illustrate the previous discussion in the particular case of the reduced Hartree-Fock (rHF) model, in which E xc = 0 in (3.9).We let 0 ≤ µ ∈ L 1 (R) be a nuclear density describing a homogeneous 2D material and denote by ν = ´R µ the total charge per unit surface.
We denote by the corresponding rHF energy per unit surface, and study the minimization problem Following the exact same lines as [9, Theorem 2.7], one has the following.
Theorem 3.7.The problem I rHF b admits a minimizer, and all minimizers share the same density.
We skip the proof for brevity.The uniqueness of the density comes from the fact that the problem is strictly convex in ρ G .However, unlike the case without magnetic field, F (b, •) is not strictly convex for b > 0. It is unclear to us whether the minimizer of I rHF b is unique.
We would like to write the Euler-Lagrange equations for a minimizer G * .Recall that g → F (b, g) is continuous and convex (but is not smooth).We denote by f b := ∂ g F (b, •) its subdifferential, a set-valued function defined by In our case, since the function F b is piece-wise linear, f b is explicit.From (3.8), we obtain (in the following lines, {a} denotes the singleton a) Its inverse map, noted h b , is the set-valued function so that y ∈ f b (x) iff x ∈ h b (y).One finds, for y > 0, We extend the definition of h b by setting h b (y) = 0 for y < 0.
In order to work with functions, it is useful to introduce the maps ] for all y ∈ R. The Euler-Lagrange equations for G * takes the following form (see end of the section for the proof). .Then there is λ ∈ R so that where ρ * = ρ G * is the associated density of G * , and Φ * is the mean-field potential, defined as the unique solution of the last equation.
The first equation can also be written as and means that if G * = j g j |ψ j ψ j | is the spectral decomposition of the optimizer, then ψ j is also an eigenfunction of H * for an eigenvalue ε j so that Conversely, if ε < λ is an eigenvalue of H * , then ε = ε j for some j.
In practice, for numerical purpose, one rather considers an approximation F δ b of F , which is smooth, strictly convex, and so that F δ n − F ∞ < δ.In this case, one can repeat the arguments in [9, Theorem 2.7], and the first line of (3.12) becomes Remark 3.9 (Strong magnetic fields).In the case where b > 2πν, any G ∈ G ν are positive and satisfies Tr(G) = ν, hence all eigenvalues of G are smaller than ν.In particular, This minimizer is therefore independent of b > 2πν, reflecting the fact that all electrons lie in the lowest Landau level.Following the previous lines, we deduce that G * is a rank-1 operator, of the form Proof of Proposition 3.8.Let G * be a minimizer of I rHF b , and let ρ * and Φ * be the corresponding density and mean-field potential, and set H * := − 1 2 ∆ 1 + Φ * .Recall that ρ * (hence Φ * and H * ) is uniquely defined.First, we claim that G * commutes with H * .This is a standard result in the case where the map F is smooth (say of class C 1 ), using that In our case however, the map F is only piece-wise smooth, and we need a direct proof.Let A be a finite-rank symmetric operator on L 2 (R), and set G t,A := e −itA G * e itA .Since G t,A is a unitary transformation of G * , we have In particular, Together with the fact that and the definition of H * , we deduce that Since the minimum of E rHF b is obtained for t = 0, the linear term in t must vanish, that is: where we have used cyclicity of the trace, and the fact that A is finite rank (so all operators are trace-class).Since this is true for all finite rank symmetric operators A, we deduce as wanted that [H * , G * ] = 0.
Recall that G * is positive compact (even trace class).So there is M := rank(G * ) ∈ N ∪ {∞} and an orthonormal family {ψ j } j∈[1,M ] so that g j |ψ j ψ j |, with Hψ j = ε j ψ j , and with g 1 ≥ g 2 ≥ • • • ≥ g M > 0. The orthonormal family {ψ } 1≤j≤M spans Ran(G * ).For two indices (i, j) ∈ [1, M ], we consider the operator For t small enough (|t| < min{g i , g j }), the operator G (i,j) t is positive, and with Tr(G In addition, we have By optimality of G * , this quantity is always positive.Taking the limit t → 0 + , and recalling the definition of f ± b in (3.11), we deduce that for this eigenvalue, and there is equality.In other words, we have proved that It remains to prove the result on Ker(G * ).Let ψ ∈ Ker(G * ).This time, we consider the perturbed state which is in G ν for all 0 ≤ t ≤ g j .Taking the limit t → 0 + , and reasoning as before, we get where we took the supremum in 1 ≤ j ≤ M in the last inequality.We deduce that, for all ψ ∈ Ker(G * ),

Models with spin
In this section, we explain how to extend our results to the case where the spin is taken into account.In this case, the density matrix is an operator γ ∈ S(L 2 (R 3 , C 2 )) satisfying the Pauli-principle 0 ≤ γ ≤ 1.Such an operator can be decomposed as a 2 × 2 matrix of the form The spin-density 2 × 2 matrix of γ is R γ (x) = γ(x, x) (see [7] for details), and the total density is ρ γ (x) = R ↑↑ γ (x) + R ↓↓ γ (x) = γ ↑↑ (x, x) + γ ↓↓ (x, x).The kinetic energy operator is now the Pauli operator , where σ contains the Pauli matrices.In the case of constant magnetic field B = (0, 0, b) with the gauge A = b(0, x 1 , 0), this operator becomes In what follows, we denote by B := b 0 0 −b .This term corresponds to the Zeeman term.There are several ways to read this operator.Indeed, we have In the decomposition (4.1) (resp.(4.2)), we split the Pauli operator as .) First splitting: keeping the spin structure.This splitting is useful if one wants to keep track of the spin-density structure of R γ as a function of x 3 .This happens for instance when using exchange-correlations functionals which depend explicitly on the spin, as in the LSDA model [17,8].In this case, we consider one-body density matrices γ of the form Remark 4.1.One can prove that the set of such states is smaller than the set of states γ commuting with P A 2 and the magnetic translations m R .This comes from the fact that the family (m R ) is not irreducible on E n := Ran P n .
The conclusion of Theorem 3.1 still holds in the spin setting, and the proof is similar with only minor modifications: we now consider the G matrix defined by and the optimal m * j functions are given by m * j (0) = { We obtain that for any representable density ρ. with the new functional A plot of F spin (b = 1, g) is displayed in Figure 1.The constant π 2 in the π 2 g 2 term is the two-dimensional Thomas-Fermi constant, when the spin of the electron is included.Again, we see that magnetic effects gives a correction to the Thomas-Fermi approximation, but this time, including the Zeeman term, the magnetic energy is always lower than the no-magnetic one, and we have

Remark 2 . 2 . 2 −x 2 x 1 ,
The definition (2.4) is slightly different from the classical Wigner transform (see for example [18, Chapter 2]) which is rather adapted to study Landau operator with the gauge Ã = b for b = 1.A gauge transformation links the two transforms.

Figure 1 below. Proposition 3 . 2 .
The function F in (3.2) is continuous and satisfies

for n ≥ 1 .
This is similar to what we have studied in the present article, but the γ n operators now act on L 2 (R, C 2 ).Introducing the operator G := b 2π γ n and the setG spin := G ∈ L 2 (R, C 2 ), G ≥ 0 , we obtain as before that, for all representable spin-density 2 × 2 matrix R, we haveinf γ of the form (4.3) ∆ 1 G) + Tr(F (b, G)) + b ˆR(R ↑↑ − R ↓↓ ).The last term is the Zeeman term, where we have used thatTr(BG) = bTr(G ↑↑ − G ↓↓ ) = b ´(R ↑↑ − R ↓↓ ).Second splitting: loosing the spin structure.If one is only interested in keeping the total density ρ γ instead of the full spin-density 2 × 2 matrix R γ , one should instead use the spectral decomposition of the operator P A 2 := L A 2 + B. This one is easily deduced from the one of L A 2 in Proposition 2.1.Using that ε n = b(2n + 1), we obtain P A 2 = ∞ n=0 ε n P n with ε n = 2nb, The lowest Landau level has now energy ε 0 = 0, and is only occupied by spin-down electrons.The corresponding eigenspace has density ρ P 0 = b 2π .The other Landau levels are «doubly» occupied, with density ρ Pn = b π .In this case, we consider density matrices γ of the form (4.4) γ = ∞ n=0 P n ⊗ γ n , γ n ∈ S(L 2 (R, C)), 0 ≤ γ n ≤ 1.
2, Proposition 2.3.8].Remark 2.6.The hypothesis ηL A 2 = L A 2 η is not a consequence of the commutation with the operators (m R ) R .Indeed, consider for a normalized ζ ∈ L 2 (R), the projector P ζ onto the vectorial space