# Density Functional Theory for two-dimensional homogeneous materials with magnetic fields

David Gontier, Salma Lahbabi, Abdallah Maichine

► To cite this version:

David Gontier, Salma Lahbabi, Abdallah Maichine. Density Functional Theory for two-dimensional homogeneous materials with magnetic fields. Journal of Functional Analysis, 2023, 285 (9), pp.110100. 10.1016/j.jfa.2023.110100. hal-03881872

HAL Id: hal-03881872

https://hal.science/hal-03881872v1

Submitted on 11 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# DENSITY FUNCTIONAL THEORY FOR TWO-DIMENSIONAL HOMOGENEOUS MATERIALS WITH MAGNETIC FIELDS

DAVID GONTIER, SALMA LAHBABI, AND ABDALLAH MAICHINE

ABSTRACT. 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.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>Acknowledgments</td><td>4</td></tr><tr><td>2. States commuting with magnetic translations</td><td>4</td></tr><tr><td>2.1. Two dimensional Landau operator</td><td>4</td></tr><tr><td>2.2. Magnetic translations</td><td>6</td></tr><tr><td>2.3. Diagonalisation of states commuting with magnetic translations</td><td>7</td></tr><tr><td>3. Reduction of the kinetic energy, and applications</td><td>9</td></tr><tr><td>3.1. Reduction of the kinetic energy</td><td>9</td></tr><tr><td>3.2. Some properties of the function F</td><td>11</td></tr><tr><td>3.3. Reduced DFT models</td><td>12</td></tr><tr><td>3.4. The reduced Hartree-Fock case</td><td>13</td></tr><tr><td>4. Models with spin</td><td>17</td></tr><tr><td>References</td><td>20</td></tr></table>

## 1. 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

<footnote>Date: July 18, 2023.</footnote>
(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 $\mu$ which is equidistributed in the first two dimensions: $\mu(x_1, x_2, x_3) = \mu(x_3)$, and with a constant perpendicular magnetic field $\mathbf{B} = b\mathbf{e}_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

$$ \mathcal{P} := \{ \gamma \in \mathcal{S}^1(L^2(\mathbb{R}^3)), \ 0 \le \gamma \le 1, \ \forall \mathbf{R} \in \mathbb{R}^2, \ \tau_{\mathbf{R}}\gamma = \gamma\tau_{\mathbf{R}} \} $$

denote the set of one-body density matrices which commute with all $\mathbb{R}^2$ translations. Here, $\mathcal{S}^1(L^2(\mathbb{R}^3))$ stands for locally trace class self-adjoint operators with finite trace per unit surface $\mathrm{Tr}(\gamma) < \infty$ (see Section 2.3). For $\mathbf{R} \in \mathbb{R}^2 \subset \mathbb{R}^3$, we have denoted by $\tau_{\mathbf{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

$$ \mathcal{G} := \{ G \in \mathcal{S}^1(L^2(\mathbb{R})), \ G \ge 0 \}, $$

where $\mathcal{S}^1(L^2(\mathbb{R}))$ refers to the space of trace class self-adjoint operators on $L^2(\mathbb{R})$. We have proved in [9] that, for any (representable) density $\rho = \rho(x_3)$ depending only on the third variable, we have

$$ (1.1) \quad \inf_{\substack{\gamma \in \mathcal{P} \\ \rho\gamma = \rho}} \left\{ \frac{1}{2} \mathrm{Tr}(-\Delta_3 \gamma) \right\} = \inf_{\substack{G \in \mathcal{G} \\ \rho_G = \rho}} \left\{ \frac{1}{2} \mathrm{Tr}(-\Delta_1 G) + \pi \mathrm{Tr}(G^2) \right\}, $$

where $\Delta_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

$$ (1.2) \quad \inf \left\{ \frac{1}{2} \mathrm{Tr}(-\Delta_1 G) + \pi \mathrm{Tr}(G^2) + \frac{1}{2} \mathcal{D}_1(\rho_G - \mu) + E^{\mathrm{xc}}(\rho_G), \ G \in \mathcal{G} \right\}, $$

where $\mathcal{D}_1(\cdot)$ is the one-dimensional Coulomb interaction energy (see [9, Section 3.1] for a discussion about this term), and $E^{\mathrm{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 $\pi \mathrm{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^{\mathrm{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

$$ \frac{1}{2} \mathrm{Tr}((-i\nabla_3 + \mathbf{A})^2 \gamma), $$
where $\mathbf{A} = b(0, x_1, 0)$ is a vector potential so that $\text{curl } \mathbf{A} = \mathbf{B} = be_3$. We chose a gauge which is not symmetric, but which will simplify some computations. The Laplacian operator $-\Delta_3$ has been replaced by the Landau operator

$$ \mathbf{L}_3^\mathbf{A} := (-i\nabla_3 + \mathbf{A})^2 = \mathbf{L}_2^\mathbf{A} - \partial_{x_3x_3}^2, \quad \text{with} \quad \mathbf{L}_2^\mathbf{A} = -\partial_{x_1x_1}^2 + (-i\partial_{x_2} + bx_1)^2. $$

In analogy with [9], we only consider states commuting with the so-called magnetic translations $\mathbf{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)

$$ \mathcal{P}^{\mathbf{A}} := \{ \gamma \in \mathcal{G}^1(L^2(\mathbb{R}^3)), 0 \le \gamma \le 1, \forall \mathbf{R} \in \mathbb{R}^2, \mathbf{m}_R \gamma = \gamma \mathbf{m}_R, \mathbf{L}_2^{\mathbf{A}} \gamma = \gamma \mathbf{L}_2^{\mathbf{A}} \}. $$

In Theorem 2.7, we prove that any $\gamma \in \mathcal{P}^{\mathbf{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

$$ (1.4) \qquad \inf_{\substack{\gamma \in \mathcal{P}^{\mathbf{A}} \\ \rho\gamma = \rho}} \left\{ \frac{1}{2} \mathrm{Tr}(\mathbf{L}_3^{\mathbf{A}}\gamma) \right\} = \inf_{\substack{G \in \mathcal{G} \\ \rho_G = \rho}} \left\{ \frac{1}{2} \mathrm{Tr}(-\Delta_1 G) + \mathrm{Tr}(F(b, G)) \right\} $$

with the penalization term $F$ defined by

$$ F(b, g) := \pi g^2 + \frac{b^2}{4\pi} \left\{ \frac{2\pi g}{b} \right\} \left( 1 - \left\{ \frac{2\pi g}{b} \right\} \right), $$

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

$$ E_{\text{kin}}^{\text{TF}}(b, \rho) := \int_{\mathbb{R}^2} F(b, \rho(\mathbf{x})) \mathrm{d}\mathbf{x}. $$

For $b = 0$, we recover the usual two-dimensional TF kinetic energy $\int \pi \rho^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))

$$ \inf \left\{ \frac{1}{2} \mathrm{Tr}(-\Delta_1 G) + \mathrm{Tr}(F(b, G)) + \frac{1}{2} \mathcal{D}_1(\rho_G - \mu) + E^{\mathrm{xc}}(b, \rho_G), \ G \in \mathcal{G} \right\}. $$
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 $E^{xc} = 0$.

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 \in R^2}$ with some of their properties, and we characterize the states in $\mathcal{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”.

## 2. STATES COMMUTING WITH MAGNETIC TRANSLATIONS

In this section, we prove that states $\gamma \in \mathcal{P}^A$ have a particular structure.

**2.1. Two dimensional Landau operator.** We start by recalling some classical facts about the (two-dimensional) Landau operator $\mathbf{L}_2^A = -\partial_{x_1x_1}^2 + (-i\partial_{x_2} + bx_1)^2$. In what follows, we assume $b \neq 0$. For $n \in \mathbb{N}_0 = \{0, 1, 2, \dots\}$, we introduce the function $\varphi_n: \mathbb{R} \to \mathbb{R}$ defined by

$$ (2.1) \qquad \varphi_n(x) := a_n |b|^{1/4} \mathcal{H}_n(\sqrt{|b|}x) e^{-\frac{1}{2}|b|x^2}, $$

where $\mathcal{H}_n(x) = (-1)^n e^{x^2 \frac{d^n}{dx^n}} e^{-x^2}$ refers to the n-th Hermite polynomial and $a_n = (2^n n!)^{-1/2} / \pi^{1/4}$ is a normalization constant so that $\|\varphi_n\|_{L^2(\mathbb{R})} = 1$.

**Proposition 2.1.** The operator $\mathbf{L}_2^A$ has purely discrete spectrum

$$ (2.2) \qquad \sigma(\mathbf{L}_2^A) = b(2\mathbb{N}_0 + 1). $$

The eigenvalue $\varepsilon_n := b(2n + 1)$ is of infinite multiplicity, with eigenspace

$$ (2.3) \qquad E_n := \ker(\mathbf{L}_2^A - \varepsilon_n) = \{W(\varphi_n, g), g \in L^2(\mathbb{R})\}, $$

where $W$ is a Wigner type transform defined on $L^2(\mathbb{R}) \times L^2(\mathbb{R})$ by

$$ (2.4) \qquad W(\varphi, g)(\mathbf{x}) := \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} e^{-ikx_2} \varphi\left(x_1 - \frac{k}{b}\right) g(k)dk. $$

In particular, the spectral projector $\mathbf{P}_n$ onto $\ker(\mathbf{L}_2^A - \varepsilon_n)$ has kernel

$$ (2.5) \qquad \mathbf{P}_n(\mathbf{x}; \mathbf{y}) = \frac{1}{2\pi} \int_{\mathbb{R}} e^{-ik(x_2 - y_2)} \varphi_n\left(x_1 - \frac{k}{b}\right) \varphi_n\left(y_1 - \frac{k}{b}\right) dk. $$

Its density $\rho_{\mathbf{P}_n}(\mathbf{x}) := \mathbf{P}_n(\mathbf{x}; \mathbf{x}) = \frac{b}{2\pi}$ is constant and independent of n.
**Remark 2.2.** 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 $\tilde{\mathbf{A}} = \frac{b}{2} \begin{pmatrix} -x_2 \\ x_1 \end{pmatrix}$, for $b = 1$. A gauge transformation links the two transforms.

*Proof.* First, we remark that $\mathbf{L}_2^A$ commutes with all translations in the $x_2$-direction. We introduce the Fourier transform with respect to the $x_2$-variable

$$ (2.6) \qquad \mathcal{F}[f](x_1, k) := \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} f(x_1, x_2) e^{ikx_2} dx_2 $$

and its inverse

$$ \mathcal{F}^{-1}[\phi](x_1, x_2) := \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} \phi(x_1, k) e^{-ikx_2} dk. $$

We have

$$ (2.7) \qquad \mathcal{F}\mathbf{L}_2^A \mathcal{F}^{-1} = \int_{\mathbb{R}}^{\oplus} \mathfrak{h}_{b,k} dk, \quad \text{with} \quad \mathfrak{h}_{b,k} := -\partial_{xx}^2 + (bx - k)^2. $$

The operator $\mathfrak{h}_{b,k}$ is a translation of the harmonic oscillator $\mathfrak{h} := -\partial_{xx}^2 + b^2 x^2$, whose spectral decomposition is

$$ \mathfrak{h} = \sum_{n=0}^{\infty} \varepsilon_n |\varphi_n\rangle\langle\varphi_n|, $$

with $\varepsilon_n = b(2n + 1)$ and $\varphi_n$ as defined in (2.1). We deduce that the spectral decomposition of $\mathfrak{h}_{b,k}$ is

$$ \mathfrak{h}_{b,k} := \sum_{n=0}^{\infty} \varepsilon_n |\varphi_n(\cdot - \frac{k}{b})\rangle\langle\varphi_n(\cdot - \frac{k}{b})|, $$

which proves (2.2). Using (2.7), we see that $W$ is an eigenfunction of $\mathbf{L}_2^A$, corresponding to eigenvalue $\varepsilon_n$, if and only if it is of the form

$$ W = \mathcal{F}^{-1} \left[ (x_1, k) \mapsto g(k) \varphi_n \left( x_1 - \frac{k}{b} \right) \right] = W(\varphi_n, g), $$

where $g \in L^2(\mathbb{R})$. Thus

$$ E_n := \ker(\mathbf{L}_2^A - \varepsilon_n) = \{W(\varphi_n, g), g \in L^2(\mathbb{R})\}. $$

To compute the kernel $\mathbf{P}_n(\mathbf{x}, \mathbf{y})$ of the projector on $E_n$, we use the Moyal identity, that we recall here.

**Proposition 2.3.** Let $f_1, g_1, f_2, g_2$ in $L^2(\mathbb{R})$. Then, $W(f_1, g_1), W(f_2, g_2) \in L^2(\mathbb{R}^2)$ and (Moyal identity)

$$ (2.8) \qquad \langle W(f_1, g_1), W(f_2, g_2) \rangle_{L^2(\mathbb{R}^2)} = \langle f_1, f_2 \rangle_{L^2(\mathbb{R})} \langle g_1, g_2 \rangle_{L^2(\mathbb{R})}. $$
Proof. We first prove the result for $f_1, g_1, f_2, g_2 \in C_0^\infty(\mathbb{R})$ and conclude by density. Parseval identity gives

$$ \begin{aligned} \langle W(f_1, g_1), W(f_2, g_2) \rangle_{L^2(\mathbb{R}^2)} &= \iint_{\mathbb{R}^2} (\overline{f_1}f_2) \left(x_1 - \frac{k}{b}\right) (\overline{g_1}g_2)(k) dx_1 dk \\ &= \langle f_1, f_2 \rangle_{L^2(\mathbb{R})} \langle g_1, g_2 \rangle_{L^2(\mathbb{R})}. \end{aligned} $$

□

In particular, if $(\psi_m)_{m \in \mathbb{N}}$ is any basis of $L^2(\mathbb{R})$, then $\{W(\varphi_n, \psi_m)\}_{m \in \mathbb{N}}$ is a basis of $E_n$. Notice that, for any fixed $\mathbf{x} \in \mathbb{R}^2$, we have

$$ W(\varphi, g)(\mathbf{x}) = \langle \varphi_{\mathbf{x}}, g \rangle_{L^2(\mathbb{R})}, \quad \text{with} \quad \varphi_{\mathbf{x}}(k) := \frac{1}{\sqrt{2\pi}} e^{ikx_2} \bar{\varphi} \left( x_1 - \frac{k}{b} \right). $$

Hence

$$ \mathbf{P}_n(\mathbf{x}; \mathbf{y}) = \sum_{m \in \mathbb{N}} W(\varphi_n, \psi_m)(\mathbf{x}) \overline{W(\varphi_n, \psi_m)}(\mathbf{y}) = \sum_{m \in \mathbb{N}} \langle \varphi_{n, \mathbf{x}}, \psi_m \rangle \langle \psi_m, \varphi_{n, \mathbf{y}} \rangle. $$

Using that $\sum_m |\psi_m\rangle\langle\psi_m| = \mathbb{I}_{L^2}$, we obtain $\mathbf{P}_n(\mathbf{x}; \mathbf{y}) = \langle\varphi_{n,\mathbf{x}}, \varphi_{n,\mathbf{y}}\rangle$, which, given that $\varphi_n$ is real-valued, gives (2.5). The density of $\mathbf{P}_n$ is thus

$$ \rho_{\mathbf{P}_n}(\mathbf{x}) = \mathbf{P}_n(\mathbf{x}; \mathbf{x}) = \frac{1}{2\pi} \int_{\mathbb{R}} \left| \varphi_n \left( x_1 - \frac{k}{b} \right) \right|^2 dk = \frac{b}{2\pi}. $$

□

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

$$ \mathbf{L}_2^\mathbf{A} := p_{\mathbf{A},1}^2 + p_{\mathbf{A},2}^2, \quad \text{with} \quad p_{\mathbf{A},1} := -i\partial_{x_1}, \quad p_{\mathbf{A},2} := -i\partial_{x_2} + bx_1. $$

The operators $p_{\mathbf{A},1}$ and $p_{\mathbf{A},2}$ do not commute, and do not commute with $\mathbf{L}_2^\mathbf{A}$. Actually, we have

$$ [p_{\mathbf{A},1}, p_{\mathbf{A},2}] = -ib, \quad [p_{\mathbf{A},1}, \mathbf{L}_2^\mathbf{A}] = -2ibp_{\mathbf{A},2}, \quad [p_{\mathbf{A},2}, \mathbf{L}_2^\mathbf{A}] = 2ibp_{\mathbf{A},1}. $$

However, introducing the dual momentum operators

$$ \tilde{p}_{\mathbf{A},1} := -i\partial_{x_1} + bx_2, \quad \tilde{p}_{\mathbf{A},2} := -i\partial_{x_2}, $$

we can check that $[\tilde{p}_{\mathbf{A},1}, \mathbf{L}_2^\mathbf{A}] = [\tilde{p}_{\mathbf{A},1}, \mathbf{L}_2^\mathbf{A}] = 0$. The magnetic translation $\mathbf{m}_R, \mathbf{R} \in \mathbb{R}^2$, is the unitary operator

$$ \begin{aligned} \mathbf{m}_R &= \exp(-\frac{i}{2}bR_1R_2) \exp(-i\tilde{\mathbf{p}}_A \cdot \mathbf{R}) \\ &= \exp(-\frac{i}{2}bR_1R_2) \exp(-i(\tilde{\mathbf{p}}_{A,1}R_1 + \tilde{\mathbf{p}}_{A,2}R_2)). \end{aligned} $$

Note that we have added a phase factor in order to match the usual convention. Using the Baker-Campbell-Hausssdorf formula and the fact that $[\tilde{\mathbf{p}}_{\mathbf{A},1}, \tilde{\mathbf{p}}_{\mathbf{A},2}] = ib$ commutes with all operators, we obtain the explicit expression

$$ \mathbf{m}_R = \exp(-ibR_1x_2)\tau_R, \quad \text{that is} \quad (\mathbf{m}_R f)(\mathbf{x}) = \exp(-ibR_1x_2)f(\mathbf{x} - \mathbf{R}), $$
where $\tau_R f(x) := f(x - R)$ is the usual translation operator. By construction, the magnetic translations commute with $L_2^A$ and $P_n$, but they do not commute among them. Actually, we have

$$m_R m_{\tilde{R}} = e^{ibR_2 \tilde{R}_1} m_{R+\tilde{R}} \quad \text{and} \quad m_R^* = m_R^{-1} = e^{ibR_1 R_2} m_{-R}.$$

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

$$(2.9) \qquad \forall \Psi \in E_n \setminus \{0\}, \quad E_n = \operatorname{span}\{m_R \Psi : R \in \mathbb{R}^2\}.$$

*Proof.* Assume otherwise, and let $\Psi$ so that

$$\tilde{E}_n(\Psi) := \operatorname{span}\{m_R \Psi : R \in \mathbb{R}^2\} \subsetneq E_n.$$

Then there is $\Phi \in E_n \setminus \{0\}$ so that $\Phi \perp \tilde{E}_n(\Psi)$. Let $f, g \in L^2(\mathbb{R})$ so that $\Psi = W(\varphi_n, g)$ and $\Phi = W(\varphi_n, f)$. Using the Moyal identity, and the fact that

$$m_R W(\varphi_n, g) = W(\varphi_n, t_R g), \quad \text{with} \quad t_R g : k \mapsto e^{-ibR_1 R_2} e^{ikR_2} g(k - bR_1), \\ \text{the condition } \langle \Phi, m_R \Psi \rangle = 0 \text{ for all } R \in \mathbb{R}^2 \text{ reads}$$

$$\forall \mathbf{R} \in \mathbb{R}^2, \quad \langle f, t_{\mathbf{R}} g \rangle = 0, \quad \text{hence} \quad \int_{\mathbb{R}} \bar{f}(k) e^{ikR_2} g(k - bR_1) dk = 0.$$

Applying the inverse Fourier transform to $k \mapsto \bar{f}(k)g(k - R_1) = 0$ shows that $\bar{f}(k)g(k - R_1) = 0$ a.e. for all $R_1 \in \mathbb{R}$. Squaring and integrating in $R_1$ gives $f = 0$, a contradiction. $\square$

**2.3. 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 $\tau_R$, then it commutes with the Laplacian operator. In the magnetic case, there are operators which commute with all magnetic translations $(m_R)_{R \in \mathbb{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.

**Proposition 2.5.** Let $\eta \in S(L^2(\mathbb{R}^2))$ be such that $\eta m_R = m_R \eta$ for all $R \in \mathbb{R}^2$ and $\eta L_A = L_A \eta$. Then, there is a family of real numbers $(\lambda_n)_{n \in \mathbb{N}_0}$ so that

$$(2.10) \qquad \eta = \sum_{n \in \mathbb{N}_0} \lambda_n P_n.$$
If $\eta$ is a locally trace class operator, then its density is constant

$$ \rho_{\eta} = \frac{b}{2\pi} \sum_{n \in \mathbb{N}_0} \lambda_n. $$

**Proof.** Since $\eta$ commutes with $\mathbf{L}_2^A$, it commutes with any spectral projector $\mathbf{P}_n$, hence leaves invariant $E_n = \text{Ran}(\mathbf{P}_n)$, for all $n \in \mathbb{N}$. The operator $\eta_n := \mathbf{P}_n\eta\mathbf{P}_n \in S(E_n)$ commutes with all magnetic translations. Since the family $\{\mathbf{m}_R\}_R$ is irreducible, it implies that $\eta_n$ is proportional to the identity on $E_n$, hence is of the form $\eta_n = \lambda_n \mathbb{I}_{E_n}$. This is a kind of Schur's Lemma, see [2, Proposition 2.3.8]. $\square$

**Remark 2.6.** The hypothesis $\eta\mathbf{L}_2^A = \mathbf{L}_2^A\eta$ is not a consequence of the commutation with the operators $(\mathbf{m}_R)_R$. Indeed, consider for a normalized $\zeta \in L^2(\mathbb{R})$, the projector $P_\zeta$ onto the vectorial space

$$ E_{\zeta} := \{W(\zeta, f), f \in L^2(\mathbb{R})\}. $$

Since $\mathbf{m}_R W(\zeta, f) = W(\zeta, \mathbf{t}_R f) \in E_\zeta$, the set $E_\zeta$ is invariant by $\mathbf{m}_R$, hence $P_\zeta$ commutes with all magnetic translations. However, we have, using the decomposition (2.7) that

$$ \begin{aligned} \mathbf{L}_2^A [W(\zeta, f)] &= \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} e^{-ikx_2} \left[ (-\partial_{x_1 x_1}^2 + (bx_1 - k)^2) \zeta \left(x_1 - \frac{k}{b}\right) \right] f(k) dk \\ &= W(\tilde{\zeta}, f) \end{aligned} $$

with $\tilde{\zeta} := (-\partial_{xx}^2 + b^2x^2) \zeta$. So $\mathbf{L}_2^A E_\zeta = E_{\tilde{\zeta}}$, and $\mathbf{L}_2^A$ leaves $E_\zeta$ invariant iff $\tilde{\zeta}$ is collinear to $\zeta$. This happens if and only if $\zeta$ is an eigenstate of the harmonic oscillator, that is $\zeta = \varphi_n$ for some $n \in \mathbb{N}_0$.

The three-dimensional analogue of the previous Proposition reads as follows.

**Theorem 2.7.** Let $\gamma$ be a bounded operator on $L^2(\mathbb{R}^3)$ commuting with all magnetic translations $\mathbf{m}_R$ and with $\mathbf{L}_2^A \otimes \mathbb{I}$. Then, there exists a family $(\gamma_n)_{n \in \mathbb{N}}$ of bounded operators on $L^2(\mathbb{R})$, with $\|\gamma_n\| \le \|\gamma\|$, and so that

$$ (2.11) \qquad \gamma = \sum_{n=0}^{\infty} \mathbf{P}_n \otimes \gamma_n. $$

If $\gamma$ is a locally trace class operator, then its density depends only on $x_3$

$$ \rho_{\gamma}(\mathbf{x}) = \rho_{\gamma}(x_3) = \frac{b}{2\pi} \sum_{n=0}^{\infty} \rho_{\gamma_n}(x_3). $$

**Proof.** Let us consider two fixed test functions $\phi, \psi \in L^2(\mathbb{R})$, and define the operator $\eta_{\phi, \psi}$ on $L^2(\mathbb{R}^2)$ by

$$ \forall f, g \in L^2(\mathbb{R}^2), \quad \langle f, \eta_{\phi, \psi}, g \rangle_{L^2(\mathbb{R}^2)} := \langle f \otimes \phi, \gamma(g \otimes \psi) \rangle_{L^2(\mathbb{R}^3)}. $$
The conditions on $\gamma$ imply that $\eta_{\phi,\psi}$ is a bounded self-adjoint operator that commutes with all $m_R$ and with $L_A$. Thus, using Proposition 2.5, $\eta_{\phi,\psi}$ can be decomposed as

$$ \eta_{\phi,\psi} = \sum_{n \in \mathbb{N}_0} \lambda_n(\phi, \psi) \mathbf{P}_n. $$

Since $\eta$ is a bounded operator, we have for any normalized $\Phi_n \in L^2(\mathbb{R}^2)$ in the range of $\mathbf{P}_n$,

$$ |\lambda_n(\phi, \psi)| = |\langle \Phi_n \otimes \phi, \gamma(\Phi_n \otimes \psi) \rangle_{L^2(\mathbb{R}^3)}| \le \|\gamma\|_{\text{op}} \|\psi\|_{L^2(\mathbb{R})} \|\phi\|_{L^2(\mathbb{R})}. $$

We deduce that the map $(\phi, \psi) \mapsto \lambda_n(\phi, \psi)$ is sesquilinear and bounded, with bound smaller than $\|\eta\|_{\text{op}}$. The result then follows by taking $\gamma_n$ the bounded self-adjoint operator on $L^2(\mathbb{R})$ defined by

$$ \langle \phi, \gamma_n \psi \rangle := \lambda_n(\phi, \psi). $$

Finally, for $\eta$ of the form (2.11), we obtain, using $\rho_{\mathbf{P}_n} = \frac{b}{2\pi}$,

$$ \rho_{\eta}(\mathbf{x}) = \eta(\mathbf{x}, \mathbf{x}) = \sum_{n \in \mathbb{N}} \rho_{\mathbf{P}_n}(x_1, x_2) \rho_{\gamma_n}(x_3) = \frac{b}{2\pi} \sum_{n \in \mathbb{N}} \rho_{\gamma_n}(x_3). $$

□

For an operator $\eta$ of the form (2.11), the trace per unit-surface, defined as the limit

$$ \underline{\text{Tr}}(\eta) = \lim_{L \to \infty} \frac{1}{L^2} \text{Tr}(\mathbf{1}_{\Gamma_L} \eta \mathbf{1}_{\Gamma_L}), \quad \Gamma_L = [-\frac{L}{2}, \frac{L}{2}]^2 \times \mathbb{R} $$

takes the simpler form

$$ \underline{\text{Tr}}(\eta) = \frac{b}{2\pi} \sum_{n=0}^{\infty} \text{Tr}_1(\eta_n), $$

where we have used that the density of any Landau level is $\rho_{\mathbf{P}_n} = \frac{b}{2\pi}$.

### 3. REDUCTION OF THE KINETIC ENERGY, AND APPLICATIONS

We now exploit the particular structure of states $\gamma \in \mathcal{P}^A$ to deduce their kinetic energy.

**3.1. Reduction of the kinetic energy.** Recall that the set $\mathcal{P}^A$ has been defined in (1.3) as the set of one-body density matrices $\gamma$, satisfying the Pauli principle $0 \le \gamma \le 1$, which commute with all magnetic translations, and with the Landau operator $L_2^A \otimes \mathbb{I}$. We also recall that the set of reduced states $\mathcal{G}$ is defined by

$$ \mathcal{G} := \{G \in S(L^2(\mathbb{R})), \ G \ge 0\}. $$

The main result of this section is the following.
**Theorem 3.1.** For any $\gamma \in \mathcal{P}^A$, there is an operator $G \in \mathcal{G}$ satisfying $\rho_G = \rho_\gamma$ and

$$ (3.1) \qquad \frac{1}{2} \underline{\mathrm{Tr}}(\mathbf{L}_3^A \gamma) \geq \frac{1}{2} \mathrm{Tr}(-\Delta_1 G) + \mathrm{Tr}(F(b, G)), $$

where

$$ (3.2) \qquad F(b,g) := \pi g^2 + \frac{b^2}{4\pi} \left\{ \frac{2\pi g}{b} \right\} \left( 1 - \left\{ \frac{2\pi g}{b} \right\} \right). $$

Conversely, for any $G \in \mathcal{G}$, there is $\gamma \in \mathcal{P}^A$ so that $\rho_\gamma = \rho_G$, and for which there is equality in (3.1). In particular, for any (representable) density $\rho$,

$$ (3.3) \qquad \inf_{\substack{\gamma \in \mathcal{P}^A \\ \rho_\gamma = \rho}} \left\{ \frac{1}{2} \underline{\mathrm{Tr}}(\mathbf{L}_3^A \gamma) \right\} = \inf_{\substack{G \in \mathcal{G} \\ \rho_G = \rho}} \left\{ \frac{1}{2} \mathrm{Tr}(-\Delta_1 G) + \mathrm{Tr}(F(b, G)) \right\}. $$

*Proof.* According to Theorem (2.7), any $\gamma \in \mathcal{P}^A$ can be decomposed as

$$ (3.4) \qquad \gamma = \sum_{n=0}^{\infty} \mathbf{P}_n \otimes \gamma_n \quad \text{with} \quad \gamma_n \in \mathcal{S}(L^2(\mathbb{R})), \quad 0 \le \gamma_n \le 1. $$

For a state of the form (3.4), we define the operator $G_\gamma \in \mathcal{S}(L^2(\mathbb{R}))$ by

$$ G_\gamma := \frac{b}{2\pi} \sum_{n=0}^{\infty} \gamma_n. $$

Since $\gamma_n \ge 0$, we have $G_\gamma \ge 0$ as well, so $G_\gamma \in \mathcal{G}$. Also, since $\rho_{\mathbf{P}_n}(\mathbf{x}) = \frac{b}{2\pi}$, we deduce that $\rho_G = \rho_\gamma$.

Recalling that $\mathbf{L}_3^A := \mathbf{L}_2^A \otimes \mathbb{I}_{L^2(\mathbb{R})} + \mathbb{I}_{L^2(\mathbb{R}^2)} \otimes (-\Delta_1)$, and using Proposition 2.1, we obtain that

$$ \begin{align*} \frac{1}{2} \underline{\mathrm{Tr}} (\mathbf{L}_3^A \gamma) &= \frac{1}{2} \sum_{n=0}^{\infty} \underline{\mathrm{Tr}} (\varepsilon_n \mathbf{P}_n \otimes \gamma_n) + \frac{1}{2} \sum_{n=0}^{\infty} \underline{\mathrm{Tr}} (\mathbf{P}_n \otimes (-\Delta_1 \gamma_n)) \\ &= \frac{b}{4\pi} \sum_{n=0}^{\infty} \varepsilon_n \mathrm{Tr}(\gamma_n) + \frac{b}{4\pi} \sum_{n=0}^{\infty} \mathrm{Tr}(-\Delta_1 \gamma_n) \\ &= \frac{b}{4\pi} \sum_{n=0}^{\infty} \varepsilon_n \mathrm{Tr}(\gamma_n) + \frac{1}{2} \mathrm{Tr}(-\Delta_1 G). \end{align*} $$

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 = \sum_j g_j |\psi_j\rangle\langle\psi_j|$ with $g_j > 0$ and $\sum_j g_j < \infty$. Evaluating the trace of $\gamma_n$ in the $\{\psi_j\}$ basis, and changing the order of the sums (all terms are positive), we obtain

$$ \sum_{n=0}^{\infty} \varepsilon_n \mathrm{Tr}(\gamma_n) = \sum_{j=1}^{\infty} \left( \sum_{n=0}^{\infty} \varepsilon_n \langle \psi_j, \gamma_n \psi_j \rangle \right). $$
Since $0 \le \gamma_n \le 1$, the quantity $m_j(n) := \langle \psi_j, \gamma_n \psi_j \rangle$ satisfies $0 \le m_j(n) \le 1$. In addition, we have $\frac{b}{2\pi} \sum_n m_j(n) = \frac{b}{2\pi} \langle \psi_j, \sum_n \gamma_n \psi_j \rangle = \langle \psi_j, G\psi_j \rangle = g_j$. So we have the inequality (3.5)

$$ \sum_{n=0}^{\infty} \varepsilon_n \mathrm{Tr}(\gamma_n) \ge \sum_{j=1}^{\infty} \inf_m \left\{ \sum_{n=0}^{\infty} \varepsilon_n m(n), 0 \le m(n) \le 1, \sum_{n=0}^{\infty} m(n) = \frac{2\pi g_j}{b} \right\}. $$

Since the $\varepsilon_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

$$ (3.6) \qquad m_j^*(n) = \begin{cases} 1 & \text{for all } 0 \le n \le \lfloor \frac{2\pi g_j}{b} \rfloor - 1 \\ \{\frac{2\pi g_j}{b}\} & \text{for } n = \lfloor \frac{2\pi g_j}{b} \rfloor \\ 0 & \text{otherwise.} \end{cases}. $$

We now calculate the infimum in the RHS of (3.5) using the explicit formula of the optimal function $m^*$. Recalling that $\varepsilon_n = b(2n + 1)$ and denoting by $x := \frac{2\pi g_j}{b}$, we obtain

$$ \begin{aligned} \sum_{n=0}^{\infty} \varepsilon_n m_j^*(n) &= b \sum_{n=0}^{\lfloor x \rfloor - 1} (2n + 1) + b(2\lfloor x \rfloor + 1) \{x\} \\ &= b(x^2 + \{x\}(1 - \{x\})) = \frac{4\pi}{b} F(b, g_j), \end{aligned} $$

with the function $F$ defined in (3.2). Summing in $j$ and gathering the terms gives the inequality

$$ \frac{1}{2} \mathrm{Tr}(\mathbf{L}_3^A \gamma) \ge \mathrm{Tr}(F(b, G)) + \frac{1}{2} \mathrm{Tr}(-\Delta_1 G), $$

which proves the first part of the Theorem. Conversely, given $G = \sum_j g_j |\psi_j\rangle\langle\psi_j|$ $\in \mathcal{G}$, we consider the state

$$ \gamma^* := \sum_{n=0}^{\infty} \mathbf{P}_n \otimes \gamma_n^*, \quad \text{with} \quad \gamma_n^* := \sum_{j=1}^{\infty} m_j^*(n) |\psi_j\rangle\langle\psi_j| $$

and $m_j^*$ defined as in (3.6). The operator $\gamma^*$ belongs to $\mathcal{P}^A$, satisfies $G_{\gamma^*} = G$, and gives an equality in (3.1). $\square$

**3.2. Some properties of the function F.** Let us collect some useful properties of the function $F$. A plot of $F$ is displayed in Figure 1 below.

**Proposition 3.2.** The function $F$ in (3.2) is continuous and satisfies

$$ (3.7) \qquad \pi g^2 \le F(b, g) \le \pi g^2 + \frac{b^2}{16\pi}, $$

with equality in the left for $\frac{2\pi g}{b} \in \mathbb{N}$, and equality in the right for $\frac{2\pi g}{b} \in \mathbb{N}+\frac{1}{2}$, and $F(b, g) \to \pi g^2$ as $b \to 0$. For any $b \ge 0$, the map $g \mapsto F(b, g) - \pi g^2$ is
$\frac{b}{2\pi}$ periodic and the map $g \in \mathbb{R}_+ \mapsto F(b, g)$ is piece-wise linear, increasing and convex. Finally, for all $0 \le g < \frac{b}{2\pi}$, we have $F(b, g) = \frac{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

$$ (3.8) \qquad F(b,g) = \frac{b^2}{4\pi} (x(1+2|x|) - |x| - |x|^2), $$

where we have denoted by $x := \frac{2\pi g}{b}$. When $0 \le g < \frac{b}{2\pi}$, which corresponds to $0 \le x < 1$, $F(b, g) = \frac{1}{2}bg$. $\square$

**Remark 3.3.** The left inequality of (3.7) implies $\mathrm{Tr}(F(b, G)) \ge \pi \mathrm{Tr}(G^2)$, hence, together with (1.1), that

$$ \inf_{\substack{\gamma \in \mathcal{P}^A \\ \rho\gamma = \rho}} \left\{ \frac{1}{2} \mathrm{Tr} \left( \mathbf{L}_3^A \gamma \right) \right\} \ge \inf_{\substack{\gamma \in \mathcal{P} \\ \rho\gamma = \rho}} \left\{ \frac{1}{2} \mathrm{Tr} \left( -\Delta_3 \gamma \right) \right\}. $$

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) \to \pi g^2$ as $b \to 0$, for all $g \in \mathbb{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) = \pi g^2 + \tilde{F}(b, g)$, we see that the effect of adding a magnetic field $\mathbf{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 $\mu(\mathbf{x}) = \mu(x_3)$ under a constant magnetic field $\mathbf{B} = b(0, 0, x_3)$ by a reduced state $G \in \mathcal{G}$ whose energy per unit surface is given by

$$ (3.9) \qquad \mathcal{E}(G) = \frac{1}{2}\mathrm{Tr}(-\Delta_1 G) + \mathrm{Tr}(F_b(G)) + \frac{1}{2}\mathcal{D}_1(\rho_G - \mu) + E^{\mathrm{xc}}(\rho_G). $$

Here, $E^{\mathrm{xc}}$ models is an exchange-correlation energy per unit surface, and $\mathcal{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 $f \in \mathcal{C} := \{f \in L^1(\mathbb{R}), W_f \in L^2(\mathbb{R})\}$, where $W_f(x) := \int_{-\infty}^x f$ is a primitive of $f$, we have

$$ \mathcal{D}_1(f) := 4\pi \int_{\mathbb{R}} |W_f|^2(x)dx. $$

We have proved in [9, Proposition 3.3] that the elements $f \in \mathcal{C}$ have null integral $\int f = 0$, that the map $\mathcal{C} \ni f \mapsto \mathcal{D}_1(f)$ is convex, and that

$$ \mathcal{D}_1(f) = 4\pi \iint_{(\mathbb{R}_+)^2 \times (\mathbb{R}_-)^2} \min\{|x|, |y|\} f(x)f(y)dxdy = \int_{\mathbb{R}} \Phi_f(x)f(x)dx, $$
with the mean-field potential

$$ \Phi_f(x) := 4\pi \int_{\mathbb{R}^{\pm}} \min\{|x|, |y|\} f(y) dy, \quad x \in \mathbb{R}^{\pm}. $$

The function $\Phi_f$ is continuous, and is the (unique) solution to

$$ -\Phi_f''(x) = 4\pi f, \quad \Phi_f'(x) \xrightarrow{x \to \pm\infty} 0, \quad \Phi_f(0) = 0. $$

In practice, we restrict the minimization problem to the $G$ so that $\rho_G - \mu \in \mathcal{C}$. This implies in particular the neutrality condition $\mathrm{Tr}(G) = \int \rho_G = \int \mu = \nu$. On the other hand, if $G$ is a trace class operator, then $\rho_G \in L^1$ (we assume that $\mu \in L^1$ as well), and if $\rho_G - \mu$ has null integral, then $\mathcal{D}_1(\rho_G - \mu) < \infty$ iff $\rho_G - \mu \in \mathcal{C}$.

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 \le \mu \in L^1(\mathbb{R})$ be a nuclear density describing a homogeneous 2D material and denote by $\nu = \int_{\mathbb{R}} \mu$ the total charge per unit surface.

We denote by

$$ (3.10) \qquad \mathcal{E}_b^{\mathrm{rHF}}(G) := \frac{1}{2} \mathrm{Tr}(-\Delta_1 G) + \mathrm{Tr}(F(b, G)) + \frac{1}{2} \mathcal{D}_1(\rho_G - \mu) $$

the corresponding rHF energy per unit surface, and study the minimization problem

$$ \mathcal{I}_b^{\mathrm{rHF}} := \inf \{ \mathcal{E}_b^{\mathrm{rHF}}(G), G \in \mathcal{G}^{\nu} \}, \quad \text{with} \quad \mathcal{G}^{\nu} := \{ G \in \mathcal{G}, \mathrm{Tr}(G) = \nu \}. $$

Following the exact same lines as [9, Theorem 2.7], one has the following.

**Theorem 3.7.** The problem $\mathcal{I}_b^{\mathrm{rHF}}$ 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 $\rho_G$. However, unlike the case without magnetic field, $F(b, \cdot)$ is not strictly convex for $b > 0$. It is unclear to us whether the minimizer of $\mathcal{I}_b^{\mathrm{rHF}}$ is unique.

We would like to write the Euler-Lagrange equations for a minimizer $G_*$. Recall that $g \mapsto F(b, g)$ is continuous and convex (but is not smooth). We
denote by $f_b := \partial_g F(b, \cdot)$ its subdifferential, a set-valued function defined by

$$f_b(g) := \{a \in \mathbb{R}, \forall g' \in \mathbb{R}, F(b, g') - F(b, g) \geq a(g' - g)\}.$$

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$)

$$f_b(g) = \begin{cases} \left\{ \frac{b^2}{4\pi} (2n + 1) \right\} & \text{if } \exists n \in \mathbb{N}_0, \quad n < \frac{2\pi g}{b} < n + 1 \\ \left[ \frac{b^2}{4\pi} (2n - 1), \frac{b^2}{4\pi} (2n + 1) \right] & \text{if } \frac{2\pi g}{b} = n \in \mathbb{N}_0. \end{cases}$$

Its inverse map, noted $h_b$, is the set-valued function so that $y \in f_b(x)$ iff $x \in h_b(y)$. One finds, for $y > 0$,

$$h_b(y) = \begin{cases} \left[ n \frac{b}{2\pi}, (n+1) \frac{b}{2\pi} \right] & \text{if } n := \frac{1}{2} \left( \frac{4\pi}{b^2} y - 1 \right) \in \mathbb{N}_0, \\ \left\{ n \frac{b}{2\pi} \right\} & \text{if there is } n \in \mathbb{N}_0 \text{ so that } n - 1 < \frac{1}{2} \left( \frac{4\pi}{b^2} y - 1 \right) < n. \end{cases}$$

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

$$(3.11) \qquad f_b^{\pm}(g) := \lim_{t \to 0^{\pm}} \frac{1}{t} (F(b, g+t) - F(b, g))$$

so that $f_b(g) = [f_b^-(g), f_b^+(g)]$ for all $g \in \mathbb{R}^+$. Of course, if $g \notin \frac{b}{2\pi}\mathbb{N}_0$ is a regular point, then $f_b(y) = f_b^+(y) = f_b^-(y)$. We define the maps $h_b^{\pm}$ similarly, so that $h_b(y) = [h_b^-(y), h_b^+(y)]$ for all $y \in \mathbb{R}$.

The Euler-Lagrange equations for $G_*$ takes the following form (see end of the section for the proof).

**Proposition 3.8 (Euler-Lagrange equations).** Let $G_*$ be a minimizer of $\mathcal{I}_b^{\text{RF}}$. Then there is $\lambda \in \mathbb{R}$ so that

$$(3.12) \qquad \begin{cases} h_b^-(\lambda - H_*) \leq G_* \leq h_b^+(\lambda - H_*) \\ H_* := -\frac{1}{2}\Delta_1 + \Phi_* \\ -\Phi_*'' = 4\pi(\rho_* - \mu), \quad \Phi_*'(x) \xrightarrow{x \to \pm\infty} 0, \quad \Phi_*(0) = 0, \end{cases}$$

where $\rho_* = \rho_{G_*}$ is the associated density of $G_*$, and $\Phi_*$ is the mean-field potential, defined as the unique solution of the last equation.

The first equation can also be written as

$$G_* \in h_b(\lambda - H_*),$$

and means that if $G_* = \sum_j g_j |\psi_j\rangle\langle\psi_j|$ is the spectral decomposition of the optimizer, then $\psi_j$ is also an eigenfunction of $H_*$ for an eigenvalue $\varepsilon_j$ so that

$$g_j \in h_b(\lambda - \varepsilon_j), \quad \text{or equivalently} \quad \varepsilon_j \in \lambda - f_g(g_j).$$

Conversely, if $\varepsilon < \lambda$ is an eigenvalue of $H_*$, then $\varepsilon = \varepsilon_j$ for some $j$.

In practice, for numerical purpose, one rather considers an approximation $F_b^\delta$ of $F$, which is smooth, strictly convex, and so that $\|F_n^\delta - F\|_\infty < \delta$. In
this case, one can repeat the arguments in [9, Theorem 2.7], and the first line of (3.12) becomes

$$ (F_b^\delta)'(G_*) = \lambda - H_*, \quad \text{or, equivalently} \quad G_* = [(F_b^\delta)'']^{-1}(\lambda - H_*). $$

**Remark 3.9** (Strong magnetic fields). In the case where $b > 2\pi\nu$, any $G \in \mathcal{G}^\nu$ are positive and satisfies $\mathrm{Tr}(G) = \nu$, hence all eigenvalues of $G$ are smaller than $\nu$. In particular,

$$ F_b(G) = \frac{1}{2}bG, \quad \text{hence} \quad \mathrm{Tr}(F_b(G)) = \frac{1}{2}b\nu $$

is a constant, independent of $G \in \mathcal{G}^\nu$. In this case, $G_*$ is also the minimizer of

$$ \inf \left\{ \frac{1}{2} \mathrm{Tr}(-\Delta_1 G) + \frac{1}{2} \mathcal{D}_1(\rho_G - \mu), \ G \in \mathcal{G}^\nu \right\}. $$

This minimizer is therefore independent of $b > 2\pi\nu$, 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 $G_* = \nu|\psi_*\rangle\langle\psi_*|$, with $\psi_*$ minimizing

$$ \inf \left\{ \frac{\nu}{2} \int_{\mathbb{R}} |\nabla \psi_*|^2 + \frac{1}{2} \mathcal{D}_1 (\nu |\psi_*|^2 - \mu), \ \psi_* \in L^2(\mathbb{R}), \ \|\psi_*\| = 1 \right\}. $$

*Proof of Proposition 3.8.* Let $G_*$ be a minimizer of $\mathcal{I}_b^{\mathrm{rHF}}$, and let $\rho_*$ and $\Phi_*$ be the corresponding density and mean-field potential, and set $H_* := -\frac{1}{2}\Delta_1 + \Phi_*$. Recall that $\rho_*$ (hence $\Phi_*$ 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

$$ \mathrm{Tr}(F(G_* + H)) = \mathrm{Tr}(F(G_*)) + \mathrm{Tr}(F'(G_*)H) + o(H). $$

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(\mathbb{R})$, and set $G_{t,A} := e^{-itA}G_*e^{itA}$. Since $G_{t,A}$ is a unitary transformation of $G_*$, we have $\mathrm{Tr}(G_{t,A}) = \mathrm{Tr}(G_*) = \nu$, that is $G \in \mathcal{G}^\nu$ and $\mathrm{Tr}(F(b, G_{t,A})) = \mathrm{Tr}(F(b, G_*))$. In particular,

$$ \mathcal{E}_b^{\mathrm{rHF}}(G_{t,A}) = \mathcal{E}_b^{\mathrm{rHF}}(G_*) + \frac{1}{2} (\mathrm{Tr}(-\Delta_1[G_{t,A} - G_*]) + \mathcal{D}_1(\rho_{G_{t,A}} - \mu) - \mathcal{D}_1(\rho_* - \mu)). $$

Together with the fact that

$$ G_{t,A} = G_* + it[G_*, A] + o(t), $$

and the definition of $H_*$, we deduce that

$$ \mathcal{E}_b^{\mathrm{rHF}}(G_{t,A}) = \mathcal{E}_b^{\mathrm{rHF}}(G_*) + it\mathrm{Tr}(H_*[G_*, A]) + o(t). $$
Since the minimum of $\mathcal{E}_b^{\text{rHF}}$ is obtained for $t=0$, the linear term in $t$ must vanish, that is:

$$ 
\begin{aligned}
0 &= \operatorname{Tr}(H_*[G_*, A]) = \operatorname{Tr}(H_*G_*A - H_*AG_*) \\
  &= \operatorname{Tr}(H_*G_*A - G_*H_*A) = \operatorname{Tr}([H_*, G_*]A),
\end{aligned}
$$

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 := \text{rank}(G_*) \in \mathbb{N} \cup \{\infty\}$ and an orthonormal family $\{\psi_j\}_{j \in [1, M]}$ so that

$$ G_* = \sum_{j=1}^{M} g_j |\psi_j\rangle\langle\psi_j|, \quad \text{with} \quad H\psi_j = \varepsilon_j\psi_j, $$

and with $g_1 \ge g_2 \ge \cdots \ge g_M > 0$. The orthonormal family $\{\psi_j\}_{1 \le j \le M}$ spans $\text{Ran}(G_*)$. For two indices $(i, j) \in [1, M]$, we consider the operator

$$ G_t^{(i,j)} := G_* + t (|\psi_i\rangle\langle\psi_i| - |\psi_j\rangle\langle\psi_j|). $$

For $t$ small enough ($|t| < \min\{g_i, g_j\}$), the operator $G_t^{(i,j)}$ is positive, and with $\operatorname{Tr}(G_t^{(i,j)}) = \operatorname{Tr}(G_*) = \nu$, hence $G_t^{(i,j)} \in \mathcal{G}^\nu$. In addition, we have

$$ 
\begin{aligned}
\mathcal{E}_b^{\text{rHF}}(G_t^{(i,j)}) - \mathcal{E}_b^{\text{rHF}}(G_*) = (F(b, g_i + t) - F(b, g_i)) &+ (F(b, g_j - t) - F(b, g_j)) \\
&+ t \underbrace{\operatorname{Tr}(H_*(|\psi_i\rangle\langle\psi_i| - |\psi_j\rangle\langle\psi_j|))}_{=\varepsilon_i - \varepsilon_j} + o(t).
\end{aligned}
$$

By optimality of $G_*$, this quantity is always positive. Taking the limit $t \to 0^+$, and recalling the definition of $f_b^\pm$ in (3.11), we deduce that

$$ f_b^+(g_i) + \varepsilon_i \ge f_b^-(g_j) + \varepsilon_j. $$

This inequality is valid for all $1 \le i, j \le M$, so

$$ \inf_{1 \le i \le M} (f_b^+(g_i) + \varepsilon_i) \ge \sup_{1 \le j \le M} (f_b^-(g_j) + \varepsilon_j) =: \lambda $$

For all $1 \le j \le M$, we have

$$ f_b^-(g_j) + \varepsilon_j \le \lambda \le f_b^+(g_j) + \varepsilon_j, \quad \text{hence} \quad \lambda - \varepsilon_j \in [f_b^-(g_j), f_b^+(g_j)]. $$

This is also $\lambda - \varepsilon_j \in f_b(g_j)$, or equivalently $g_j \in h_b(\lambda - \varepsilon_j)$. Note that if there is an eigenvalue $g_j \notin \frac{b}{2\pi}\mathbb{Z}$, then $f_b^-(g_j) = f_b^+(g_j)$ for this eigenvalue, and there is equality. In other words, we have proved that

$$ G_* \in h_b(\lambda - H_*) \quad \text{on} \quad \text{Ran}(G_*) $$

It remains to prove the result on $\operatorname{Ker}(G_*)$. Let $\psi \in \operatorname{Ker}(G_*)$. This time, we consider the perturbed state

$$ G_t^{(j)} := G_* + t (|\psi\rangle\langle\psi| - |\psi_j\rangle\langle\psi_j|), $$
which is in $G^\nu$ for all $0 \le t \le g_j$. Taking the limit $t \to 0^+$, and reasoning as before, we get

$$f_b^+(0) + \langle \psi, H_* \psi \rangle \ge f_b^-(g_j) + \varepsilon_j, \quad \text{hence} \quad f_b^+(0) + \langle \psi, H_* \psi \rangle \ge \lambda,$$

where we took the supremum in $1 \le j \le M$ in the last inequality. We deduce that, for all $\psi \in \text{Ker}(G_*)$,

$$f_b^+(\langle \psi, G_* \psi \rangle) = f_b^+(0) \ge \lambda - \langle \psi, H_* \psi \rangle,$$

so

$$f_b^+(G_*) \ge \lambda - H_*, \quad \text{on } \text{Ker}(G_*).$$

Together with the fact that $f_b^-(G_*) = 0$ on $\text{Ker}(G_*)$, we obtain $G_* \in h_b(\lambda - H_*)$ on $\text{Ker}(G_*)$ as well. $\square$

## 4. 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 $\gamma \in S(L^2(\mathbb{R}^3, \mathbb{C}^2))$ satisfying the Pauli-principle $0 \le \gamma \le 1$. Such an operator can be decomposed as a $2 \times 2$ matrix of the form

$$\gamma = \begin{pmatrix} \gamma^{\uparrow\uparrow} & \gamma^{\uparrow\downarrow} \\ \gamma^{\downarrow\uparrow} & \gamma^{\downarrow\downarrow} \end{pmatrix}.$$

The spin-density $2 \times 2$ matrix of $\gamma$ is $R_\gamma(\mathbf{x}) = \gamma(\mathbf{x}, \mathbf{x})$ (see [7] for details), and the total density is $\rho_\gamma(\mathbf{x}) = R_\gamma^{\uparrow\uparrow}(\mathbf{x}) + R_\gamma^{\downarrow\downarrow}(\mathbf{x}) = \gamma^{\uparrow\uparrow}(\mathbf{x}, \mathbf{x}) + \gamma^{\downarrow\downarrow}(\mathbf{x}, \mathbf{x})$.

The kinetic energy operator is now the Pauli operator

$$\mathbf{P}_3^A := [\sigma \cdot (-i\nabla + \mathbf{A})]^2,$$

where $\sigma$ contains the Pauli matrices. In the case of constant magnetic field $\mathbf{B} = (0, 0, b)$ with the gauge $\mathbf{A} = b(0, x_1, 0)$, this operator becomes

$$\mathbf{P}_3^A = \mathbf{L}_3^A \mathbb{I}_2 + \begin{pmatrix} b & 0 \\ 0 & -b \end{pmatrix} = \mathbf{L}_2^A \mathbb{I}_2 + \begin{pmatrix} b & 0 \\ 0 & -b \end{pmatrix} - \partial_{x_3 x_3}^2 \mathbb{I}_2.$$

In what follows, we denote by $\mathcal{B} := \begin{pmatrix} b & 0 \\ 0 & -b \end{pmatrix}$. This term corresponds to the Zeeman term. There are several ways to read this operator. Indeed, we have

$$L^2(\mathbb{R}^3, \mathbb{C}^2) \simeq L^2(\mathbb{R}^3) \otimes \mathbb{C}^2 \simeq L^2(\mathbb{R}^2) \otimes L^2(\mathbb{R}) \otimes \mathbb{C}^2$$

$$(4.1) \qquad \simeq L^2(\mathbb{R}^2) \otimes L^2(\mathbb{R}, \mathbb{C}^2)$$

$$(4.2) \qquad \simeq L^2(\mathbb{R}^2, \mathbb{C}^2) \otimes L^2(\mathbb{R}).$$

In the decomposition (4.1) (resp. (4.2)), we split the Pauli operator as

$$\mathbf{P}_3^A = \mathbf{L}_2^A \otimes \mathbb{I} + \mathbb{I} \otimes (\mathcal{B} - \partial_{x_3 x_3}^2), \quad (\text{resp. } \mathbf{P}_3^A = (\mathbf{L}_2^A + \mathcal{B}) \otimes \mathbb{I} + \mathbb{I} \otimes (-\partial_{x_3 x_3}^2)).$$

First splitting: keeping the spin structure. This splitting is useful if one wants to keep track of the spin-density structure of $R_\gamma$ 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 $\gamma$ of the form

$$ (4.3) \qquad \gamma = \sum_{n=0}^{\infty} \mathbf{P}_n \otimes \gamma_n, \quad \gamma_n \in S(L^2(\mathbb{R}, \mathbb{C}^2)), \quad 0 \le \gamma_n \le 1. $$

This is similar to what we have studied in the present article, but the $\gamma_n$ operators now act on $L^2(\mathbb{R}, \mathbb{C}^2)$. Introducing the operator $G := \frac{b}{2\pi} \sum \gamma_n$ and the set

$$ G^{\text{spin}} := \{G \in L^2(\mathbb{R}, \mathbb{C}^2), G \ge 0\}, $$

we obtain as before that, for all representable spin-density $2 \times 2$ matrix $R$, we have

$$ \begin{aligned} & \inf_{\substack{\gamma \text{ of the form } (4.3) \\ R_{\gamma}=R}} \left\{ \frac{1}{2} \text{Tr}(\mathbf{P}_3^A \gamma) \right\} \\ &= \inf_{\substack{G \in G^{\text{spin}} \\ R_G=R}} \left\{ \frac{1}{2} \text{Tr}(-\Delta_1 G) + \text{Tr}(F(b, G)) \right\} + b \int_{\mathbb{R}} (R^{\uparrow\uparrow} - R^{\downarrow\downarrow}). \end{aligned} $$

The last term is the Zeeman term, where we have used that $\text{Tr}(BG) = b\text{Tr}(G^{\uparrow\uparrow} - G^{\downarrow\downarrow}) = b \int(R^{\uparrow\uparrow} - R^{\downarrow\downarrow})$.

**Second splitting: loosing the spin structure.** If one is only interested in keeping the total density $\rho_\gamma$ instead of the full spin-density $2 \times 2$ matrix $R_\gamma$, one should instead use the spectral decomposition of the operator $\mathbf{P}_2^A := \mathbf{L}_2^A + B$. This one is easily deduced from the one of $\mathbf{L}_2^A$ in Proposition 2.1. Using that $\varepsilon_n = b(2n + 1)$, we obtain $\mathbf{P}_2^A = \sum_{n=0}^\infty \tilde{\varepsilon}_n \tilde{\mathbf{P}}_n$ with

$$ \tilde{\varepsilon}_n = 2nb, \quad \tilde{\mathbf{P}}_0 = \begin{pmatrix} 0 & 0 \\ 0 & \mathbf{P}_0 \end{pmatrix}, \quad \tilde{\mathbf{P}}_n = \begin{pmatrix} \mathbf{P}_{n-1} & 0 \\ 0 & \mathbf{P}_n \end{pmatrix} \quad \text{for } n \ge 1. $$

The lowest Landau level has now energy $\tilde{\varepsilon}_0 = 0$, and is only occupied by spin-down electrons. The corresponding eigenspace has density $\rho_{\tilde{\mathbf{P}}_0} = \frac{b}{2\pi}$. The other Landau levels are «doubly» occupied, with density $\rho_{\tilde{\mathbf{P}}_n} = \frac{b}{\pi}$.

In this case, we consider density matrices $\gamma$ of the form

$$ (4.4) \qquad \gamma = \sum_{n=0}^{\infty} \tilde{\mathbf{P}}_n \otimes \tilde{\gamma}_n, \quad \tilde{\gamma}_n \in S(L^2(\mathbb{R}, \mathbb{C})), \quad 0 \le \tilde{\gamma}_n \le 1. $$

**Remark 4.1.** One can prove that the set of such states is smaller than the set of states $\gamma$ commuting with $\mathbf{P}_2^A$ and the magnetic translations $\mathbf{m}_R$. This comes from the fact that the family $(\mathbf{m}_R)$ is not irreducible on $\tilde{E}_n := \text{Ran}\tilde{\mathbf{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

$$ G_{\gamma} = \frac{b}{2\pi} \left( \tilde{\gamma}_0 + 2 \sum_{n=1}^{\infty} \tilde{\gamma}_n \right) \quad \in \mathcal{G}, $$
and the optimal $m_j^*(0) = \{\frac{2\pi g_j}{b}\}$ and $m_j^*(n) = 0$ if $\frac{2\pi g_j}{b} < 1$, and, if $\frac{2\pi g_j}{b} \ge 1$,

$$m_j^*(n) = \begin{cases} 1 & \text{for } n \le \lfloor \frac{\pi}{b}g_j + \frac{1}{2} \rfloor - 1 \\ \{\frac{\pi}{b}g_j + \frac{1}{2}\} & \text{for } n = \lfloor \frac{\pi}{b}g_j + \frac{1}{2} \rfloor \\ 0 & \text{otherwise.} \end{cases}$$

Denoting by $y = \frac{\pi}{b}g_j + \frac{1}{2}$, this gives the energy

$$\frac{b}{2\pi} \sum_{n=0}^{\infty} \tilde{\epsilon}_n m_j^*(n) = \frac{b^2}{2\pi} \sum_{n=1}^{[y]-1} 2n + \frac{b^2}{\pi} [y] \{y\} = \frac{b^2}{2\pi} (y^2 - y + \{y\}(1 - \{y\})).$$

We obtain that for any representable density $\rho$.

$$\inf_{\substack{\gamma \text{ of the form (4.4)} \\ \rho\gamma = \rho}} \left\{ \frac{1}{2} \text{Tr}(\mathbf{P}_3^\mathbf{A} \gamma) \right\} = \inf_{\substack{G \in \mathcal{G} \\ \rho_G = \rho}} \left\{ \frac{1}{2} \text{Tr}(-\Delta_1 G) + \text{Tr}(F^{\text{spin}}(b, G)) \right\},$$

with the new functional

$$F^{\text{spin}}(b, g) = \frac{\pi}{2}g^2 - \frac{b^2}{8\pi} + \frac{b^2}{2\pi} \left\{ \frac{\pi g}{b} + \frac{1}{2} \right\} \left( 1 - \left\{ \frac{\pi g}{b} + \frac{1}{2} \right\} \right).$$

A plot of $F^{\text{spin}}(b=1, g)$ is displayed in Figure 1. The constant $\frac{\pi}{2}$ in the $\frac{\pi}{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

$$(4.5) \qquad \frac{\pi}{2}g^2 - \frac{b^2}{8\pi} \le F^{\text{spin}}(b,g) \le \frac{\pi}{2}g^2.$$

<figure><img src="image_2.png" /><figcaption>FIGURE 1. (Left). The map $g \mapsto F(b=1, g)$ (red) together with its lower/upper bounds in (3.7). (Right) The map $g \mapsto F^{\text{spin}}(b=1, g)$ (red) and its bounds (4.5).</figcaption></figure>
REFERENCES

[1] X. Blanc and C. Le Bris. Thomas-Fermi type theories for polymers and thin films. *Advances in Differential Equations*, 5(7-9):977 – 1032, 2000.

[2] O. Bratteli and D.W. Robinson. *Operator Algebras and Quantum Statistical Mechanics: C\*- and W\*-algebras, symmetry groups, decomposition of states*. Texts and monographs in physics. Springer-Verlag, 1987.

[3] Vivek Chaudhary, Petr Neugebauer, Omar Mounkachi, Salma Lahbabi, and A. EL Fatimy. Phosphorene - an emerging two-dimensional material: recent advances in synthesis, functionalization, and applications. *2D Materials*, 9, 05 2022.

[4] J.H. Condon. Nonlinear de Haas-van Alphen effect and magnetic domains in Beryllium. *Phys. Rev.*, 145(2):526, 1966.

[5] W.J. De Haas and P.M. Van Alphen. The dependence of the susceptibility of diamagnetic metals upon the field. *Proc. Netherlands Roy. Acad. Sci*, 33(1106):170, 1930.

[6] I. Fushiki, E.H. Gudmundsson, C.J. Pethick, and J. Yngvason. Matter in a magnetic field in the Thomas-Fermi and related theories. *Annals of Physics*, 216(1):29–72, 1992.

[7] D. Gontier. N-representability in noncollinear spin-polarized density-functional theory. *Phys. Rev. Lett.*, 111:153001, 2013.

[8] D. Gontier. Existence of minimizers for Kohn–Sham within the local spin density approximation. *Nonlinearity*, 28(1):57, 2014.

[9] D. Gontier, S. Lahbabi, and Maichine A. Density functional theory for homogeneous two dimensional materials. *Comm. Math. Phys.*, 388:1475–1505, 2021.

[10] D. Gupta, V. Chauhan, and R. Kumar. A comprehensive review on synthesis and applications of molybdenum disulfide (mos2) material: Past and recent developments. *Inorganic Chemistry Communications*, 121:108200, 2020.

[11] Y. Kaddar, G. Tiouitchi, A. El Kenz, A. Benyoussef, P. Neugebauer, S. Lahbabi, A. El Fatimy, and O. Mounkachi. Black phosphorus synthesized by the chemical vapor transport: Magnetic properties induced by oxygen atoms. Submitted for publication, 2022.

[12] L.L. Lev, I.O. Maiboroda, M.A. Husanu, and et al. k-space imaging of anisotropic 2D electron gas in GaN/GaNl high-electron-mobility transistor heterostructures. *Nature Communication*, 9(2653), 2018.

[13] E.H. Lieb and M. Loss. *Analysis*, volume 14. American Mathematical Soc., 2001.

[14] E.H. Lieb, J.P. Solovej, and J. Yngvason. Ground states of large quantum dots in magnetic fields. *Phys Rev B Condens Matter*, 51(16):10646–10665, 1995.

[15] R.S. Markiewicz, M. Meskoob, and C. Zahopoulos. Giant magnetic interaction (Condon domains) in two dimensions. *Phys. Rev. Lett.*, 54(13):1436, 1985.

[16] J. Saha and A. Dutta. A review of graphene: Material synthesis from biomass sources. *Waste and Biomass Valorization*, 13, 03 2022.

[17] U. Von Barth and L. Hedin. A local exchange-correlation potential for the spin polarized case. I. *Journal of Physics C: Solid State Physics*, 5(13):1629, 1972.

[18] M. W. Wong. *Wigner Transforms*. Universitext. Springer, New York, 1998.

[19] Y. Yngvason. Thomas-Fermi theory for matter in a magnetic field as a limit of quantum mechanics. *Letters in Mathematical Physics*, 22:107–117, 1991.
(David Gontier) CEREMADE, UNIVERSITY OF PARIS-DAUPHINE, PSL UNIVERSITY, 75016 PARIS, FRANCE & ENS/PSL UNIVERSITY, DÉPARTEMENT DE MATHÉMATIQUES ET APPLICATIONS, F-75005, PARIS, FRANCE.

Email address: gontier@ceremade.dauphine.fr

(Salma Lahbabi) EMAMI, LRI, ENSEM, UHIIC, 7 ROUTE D'EL JADIDA, B.P. 8118 OASIS, CASABLANCA; MSDA, MOHAMMED VI POLYTECHNIC UNIVERSITY, LOT 660, HAY MOULAY RACHID BEN GUERIR, 43150, MOROCCO

Email address: s.lahbabi@ensem.ac.ma

(Abdallah Maichine) DEPARTMENT OF MATHEMATICS, FACULTY OF SCIENCES, MOHAMMED V UNIVERSITY IN RABAT, 4 AVENUE IBN BATTOUTA P.B. 1014 RP, RABAT, MOROCCO & ECOLE CENTRALE CASABLANCA, BOUSKOURA, VILLE VERTE, P.B. 27182, MOROCCO

Email address: abdallahmaichine@gmail.com