RETRACTED ARTICLE: Quasi-dimensional models applied to kinetic and exchange energy density functionals

Abstract

We investigate the behavior of three-dimensional 3D exchange energy functional of density-functional theory in anisotropic systems with two-dimensional 2D character and 1D character. The local density approximation (LDA), the generalized gradient approximation (GGA), and the meta-GGA behave as functions of quantum well width. We use the infinite-barrier model (IBM) for the quantum well. In the first section, we describe the problem of three-dimensional exchange functional, in the second section we introduce the quasi-2D IBM system, in the third section we introduce the quasi-1D IBM system. Using that an exact-exchange functional provides the correct approach to the true two-dimensional limit, we want to show that the 2D limit can be considered as a constraint on approximate functionals. For the 1D limit case we also propose a new functional obtained with methods completely similar to those of 2D limit.

1 Introduction

The Kohn–Sham (KS) Density Functional Theory (DFT) [1] is the most used method for electronic structure calculations in quantum chemistry and material science [2,3,4]. In the KSDFT self-consistent scheme, the noninteracting kinetic energy functional $T_s[n]=\int d\mathbf {r}\;\tau (\mathbf {r})=\int d\mathbf {r}\; \sum _{j}^{occ}|\nabla \phi _j(\mathbf {r})|^2/2$ is treated exactly using the one-particle occupied KS orbitals $\{ \phi _j(\mathbf {r})\}$, and only the exchange-correlation (XC) energy functional $E_{xc}[n]=\int d\mathbf {r}\;n(\mathbf {r})\epsilon _{xc}(\mathbf {r})$ must be approximated [5]. Here $n(\mathbf {r})$ is the electronic density, $\epsilon _{xc}$ is the XC energy per particle, and $\tau$ is the positive defined kinetic energy density.

The XC functional $E_{xc}[n]$ contains all many-body quantum effects beyond the Hartree method, being intensively investigated for several decades. Nowadays, there are known many exact properties of $E_{xc}[n]$, such as the Görling–Levy (GL) perturbation theory [6,7,8], scaling relations due to coordinate transformations [9,10,11], gradient expansions [12,13,14,15,16,17], asymptotic behaviors [18,19,20,21,22,23,24], XC hole sum rules and on-top hole behaviors [25,26,27,28,29,30,31,32,33,34,35]. Many of these exact XC properties, have been used in the construction of XC functional approximations, that are classified on the so called Jacob’s ladder [36, 37].

Using two simple models, the quasi-2D electron gas and the quasi-1D electron gas we show a fundamental limitation of the 3 nonempirical rungs of the Jacob’s ladder, namely the local density approximation LDA and its semilocal extensions, generalized gradient approximation GGA and meta-GGA MGGA, the most widely used forms of which are worse than the LDA in the strong 2D limit. An exact-exchange functional provides the correct approach to the true two-dimensional limit. We investigate the performance of three-dimensional density functional $E_{xc}[n]$ in the quasi-two-dimensional electron gas, showing how all three semi-local approximations behave as functions of quantum well width.

2 Quasi-2D IBM system

The quasi-2D IBM system is defined by the KS potential

$$\begin{aligned} v^{KS}(x,y,z)=\left\{ \begin{array}{ll} 0, &{}\quad x\in [0,L]\\ \infty , &{}\quad {\text {otherwise}},\\ \end{array} \right. \end{aligned}$$

(1)

Then the KS orbitals are plane waves in the (yz)-plane, having the following expression

$$\begin{aligned} \Psi _{l,\mathbf {k}_\parallel }(x)=\frac{\sqrt{2}}{\sqrt{L}}\sin \Big (\frac{l \pi x}{L}\Big )\frac{\sqrt{2}}{\sqrt{A}}e^{i \mathbf {k}_\parallel \mathbf {r}_\parallel }, \end{aligned}$$

(2)

with A and $\mathbf {k}_\parallel$ are the normalization area and the wave vector in the (yz)-plane, and $l=1,2,3,..$ is the sub-band index. For the quasi-2D regime, we consider only the lowest level is occupied. Then the density is

$$\begin{aligned} n(x)=\frac{2}{L\pi (r_s^{2D})^2}\sin ^2\Big (\frac{\pi x}{L}\Big ), \end{aligned}$$

(3)

where $r_s^{2D}$ is the bulk parameter of the 2D uniform electron gas (UEG), which will be kept fixed when L is changed from the maximum value $L_{\max }$ [38] to zero. Solving the Schrödinger equation, we find the energy levels

$$\begin{aligned} E_{l,k}=\frac{1}{2}\left[ \left( \frac{l\pi }{L}\right) ^2+k_y^2+k_z^2\right] \end{aligned}$$

(4)

When only the lowest level is occupied $(E_{1,k^{2D}_F} < E_{2,0})$, the density of states of this system begins to resemble the density of states of a 2D electron gas [39]. Here $k^{2D}_F = (2n^{2D})^{1/2}$ is the two-dimensional Fermi wavevector, so

$$\begin{aligned} L< \sqrt{\frac{3}{2}}\pi r_s^{2D}=3.85r_s^{2D}=L_{\text {max}} \end{aligned}$$

(5)

Note that $r_s^{2D}=\sqrt{2}/k_F^{2D}$, where $k_F^{2D}$ is the 2D bulk Fermi wave vector.

In the limit $L\rightarrow 0$, the system approaches the 2D UEG. Varying L is equivalent to performing a non-uniform density scaling in one dimension,

$$\begin{aligned} n_\lambda (x,y,z)=\lambda n(\lambda x,y,z) \end{aligned}$$

(6)

with $\lambda \sim 1/L \rightarrow \infty$. Under this scaling, the reduced gradient $s=|\nabla n|/[2(3\pi ^2)^{1/3}n^{4/3}]$ behaves as $s_\lambda \sim \lambda ^{2/3}$.

2.1 Kinetic energy of the quasi-2DEG

The kinetic energy density (defined by $T_s=\int d\mathbf {r}\tau$) is

$$\begin{aligned}&\tau =\tau ^W+\tau ^P,\nonumber \\&\tau ^W=\frac{\pi (k_F^{2D})^2}{2L^3}\cos ^2(\frac{\pi x}{L}),\nonumber \\&\tau ^P=\frac{(k_F^{2D})^4}{4\pi L}\sin ^2(\frac{\pi x}{L}), \end{aligned}$$

(7)

where $\tau ^W$ and $\tau ^P$ are the von Weizsäcker and Pauli kinetic energy densities [40, 41]. The first of eq. (7) follows from the following equation

$$\begin{aligned} T_s[n] = T^{W}_s + \int \tau ^{TF}F_s^p(s,q)d\mathbf{r} \end{aligned}$$

(8)

where $\tau ^{TF} = (3/10)(k^{2D}_F)^2 n$ is the Thomas-Fermi kinetic energy (KE) density, with $k^{2D}_F = (3\pi ^2n)^ {1/3}$ being the Fermi wave vector, $s = |\nabla n|/(2k_F n)$ and $q = \nabla ^2 n/(4k^2_F n)$ are the reduced gradient and Laplacian, $T^W_s =\int \tau ^{T F} (5s^2/3)d \mathbf{r}$ is the von Weizsäcker kinetic energy and $F^p_s=F^{KS}_s-F^W_s$ is the Pauli KE enhancement factor.

The averaged kinetic energies per particle (defined by $T_s/N$) are

$$\begin{aligned}&T_s/N=T_s^W/N+T_s^P/N,\nonumber \\&T_s^W/N=\pi ^2/[2 L^2],\nonumber \\&T_s^P/N =(k_F^{2D})^2/4=t_s^{2D}, \end{aligned}$$

(9)

where $t_s^{2D}=\tau ^{2D}/n^{2D}$ is the 2D UEG kinetic energy per particle. Then, the Pauli KE per particle fully recovers the 2D UEG, while the von Weizsäcker part diverges as $\sim L^{-2}$, representing the short-wavelengths oscillations in the x-direction. Noting that

$$\begin{aligned} \frac{\delta T_s^W}{\delta n}=\frac{\pi ^2}{2L^2}, \end{aligned}$$

(10)

and using the Euler equation [42]

$$\begin{aligned} \frac{\delta T_s[n]}{\delta n(\mathbf{r})}+\nu _{ext}(\mathbf{r})+\nu _J(\mathbf{r};[n])+\frac{\delta E_{xc}[n]}{\delta n(\mathbf{r})}= \mu \end{aligned}$$

(11)

and Eq. (1), we conclude that

$$\begin{aligned} \frac{\delta T_s^P}{\delta n}=\pi n^{2D}=(k_F^{2D})^2/2. \end{aligned}$$

(12)

Fig. 1

Quasi-2D IBM Pauli kinetic energies per particle ($T_s^P/N$) from various KE functionals, versus $L/L_{\text {max}}$, for $r_s^{2D}=2$ (upper panel) and $r_s^{2D}=5$ (lower panel)

In Fig. 1, we show $T_s^P/N$ computed from several KE functionals, versus $L/L_{\text {max}}$. We consider the recently proposed PG1 GGA [41], PGint GGA [43], LKT GGA [44], as well as the popular Thomas–Fermi–Weizsäcker (TFW), second-order gradient expansion (GE2) [45], E00 GGA [46], and the Perdew–Constantin (PC) Laplacian-level meta-GGA [47]. We recall that PG1 and LKT are very accurate for the orbital-free DFT (OFDFT) calculations of bulk solids. On the other hand, PGint is based on the second-order gradient singularity expansion which can mimic the singularity of the jellium linear response at the wave vector $k=2k_F$. Finally, the PC meta-GGA is a very good model for the kinetic energy density $\tau$, but its functional derivative shows strong unphysical oscillations [48] A reparametrization of the PC KE functional has been proposed in Ref. [49].

The best performance in the quasi-2D regime, is found from PGint GGA, closely followed by PG1 and LKT functionals. The worst performances are given by GE2, E00 GGA, and PC meta-GGA, all of them failing badly in the strong quasi-2D regions (when $L/L_{\text {max}}\rightarrow 0$). Moreover, GE2 and E00 give wrongly negative $T_s^P/N$ when $L/L_{\text {max}}\le 0.3$.

Fig. 2

Upper panel: Pauli potential $v_s^P=\delta T_s^P/\delta n$ versus the scaled distance x/L, for the quasi-2D IBM quantum well with $L=L_{\text {max}}/2$ and $r_s^{2D}=2$. Lower panel: The averaged Pauli potential $\bar{v}_s^P=\int dx\;n v_s^P /N$ versus $L/L_{\text {max}}$, for the 2D bulk parameter $r_s^{2D}=2$

In the upper panel of Fig. 2, we show the Pauli potential $v_s^P=\delta T_s^P/\delta n$ computed from several KE functionals, using the exact density of the quasi-2D IBM with $r_s^{2D}=2$ and $L=L_{\text {max}}/2$. Due to the Euler equation, the exact curve must be a constant representing the kinetic potential of the 2D UEG. Any tested KE functional cannot give a constant Pauli potential. However, LKT and PG1 Pauli potentials have less structure in the region $0.2 \le x/L \le 0.8$. Noting that this is a very hard test for any functional, we consider that LKT and PG1 performances are quite remarkable.

In the lower panel of Fig. 2, we show the averaged Pauli potential $\bar{v}_s^P= \int dx\;n v_s^P /N$ versus $L/L_{\text {max}}$ for the case $r_s^{2D}=2$. The trend is similar with the one of Fig. 1, with PGint providing the best performance.

2.2 Exchange energy of the quasi-2DEG

The first-order density matrix of the quasi-2D IBM is

$$\begin{aligned} n_1(\mathbf {r},\mathbf {r}')=\frac{2}{\pi L}\sin (\frac{\pi x}{L})\sin (\frac{\pi x'}{L})\frac{k_F^{2D}}{\rho '}J_1(k_F^{2D}\rho '), \end{aligned}$$

(13)

where, without any loss of generality, we chose $\mathbf {r}=(x,0,0)$, and $\mathbf {r}'=(x',\rho ',\alpha )$ in cylindrical coordinates [50]. Note that $n(x)=n_1(\mathbf {r},\mathbf {r})$. We also recall that the density matrix of the 2D UEG is $n_1^{2D-UEG}(|\mathbf {r}_\parallel -\mathbf {r}'_\parallel |)= \frac{k_F^{2D}}{\pi }\frac{J_1(k_F^{2D}|\mathbf {r}_\parallel -\mathbf {r}'_\parallel |)}{|\mathbf {r}_\parallel -\mathbf {r}'_\parallel |}$, where $\mathbf {r}_\parallel$ is the 2D vector in the (yz)-plane.

We calculate the exact exchange energy (EXX) per particle $E_x/N$, where $E_x=\frac{1}{2}\int d\mathbf {r}\int d\mathbf {r}'\;n(\mathbf {r})\frac{n_x(\mathbf {r},\mathbf {r}')}{|\mathbf {r}-\mathbf {r}'|}$, with $n_x(\mathbf {r},\mathbf {r}')=-|n_1(\mathbf {r},\mathbf {r}')|^2/[2 n(\mathbf {r})]$ being the exchange hole.

Fig. 3

Exchange energy per particle ($\epsilon _x=E_x/N$) versus $L/L_{\text {max}}$ for the quasi-2D IBM with 2D bulk parameter $r_s^{2D}=2$ (upper panel), and $r_s^{2D}=5$ (lower panel)

In Fig. 3, we show a comparison between EXX, MVS meta-GGA [51], SCAN meta-GGA [52], and Q2D GGA [53] exchange energies per particle ($\epsilon _x=E_x/N$), in the whole quasi-2D regime ($0\le L/L_{\text {max}}\le 1$). We recall that all tested semilocal functionals (MVS, SCAN, and Q2D) have been constructed from the quasi-2D condition $F_x\propto s^{-1/2}$ when $s\rightarrow \infty$, that gives a finite value of $\epsilon _x=E_x/N$ in the 2D limit [54]. Here $F_x$ is the exchange enhancement factor, defined by $E_x=\int d\mathbf {r}\; e_x^{LDA} F_x$, with $e_x^{LDA}=-3(3\pi ^2)^{1/3}n^{4/3}/[4\pi ]$.

The best accuracy is obtained with Q2D GGA that, by construction, recovers the 2D LDA exchange energy per particle, when $L\rightarrow 0$. We note that EXX behaves as 2D LDA exchange in the 2D limit, similar with the Pauli kinetic energy that also recovers the 2D LDA kinetic energy. This fact shows that the short-wavelengths oscillations in the x-direction which are essential for the divergence of the von Weizsäcker KE per particle, do not contribute to the exchange energy.

3 Quasi-1D IBM system

The quasi-1D IBM system is defined by the KS potential

$$\begin{aligned} v^{KS}(x,y,z)=\left\{ \begin{array}{ll} 0, &{}\quad \rho \in [0,L]\\ \infty , &{}\quad {\text {otherwise}},\\ \end{array} \right. \end{aligned}$$

(14)

where $\rho =\sqrt{x^2+y^2}$ is the radial distance from the z-axis. The Kohn–Sham orbitals have the form

$$\begin{aligned} \Psi _{l,k_z}(\rho )=\frac{1}{\sqrt{L_z}}e^{ik_z z}\phi _l(\rho ), \end{aligned}$$

(15)

where $L_z$ is a normalization length, and $\phi _l(\rho )$ satisfies the equation

$$\begin{aligned} -\frac{1}{2}\Big [\frac{d^2 \phi _l}{d\rho ^2}+\frac{1}{\rho }\frac{d\phi _l}{d\rho }\Big ]=\epsilon _l \phi _l(\rho ), \end{aligned}$$

(16)

whose solutions are of the form $\phi _l(\rho )=A J_0(\frac{x_{0l}}{L}\rho )$, where $x_{0l}$ is the lth zero of the Bessel function $J_0$. The total energy is

$$\begin{aligned} E_{l,k_z}=\frac{x_{ol}^2}{2 L^2}+\frac{k_z^2}{2}. \end{aligned}$$

(17)

The quasi-1D regime is defined by the condition $E_{1,k_F^{1D}}\le E_{2,0}$, which defines the maximum length

$$\begin{aligned} L\le L_{\text {max}}=\sqrt{x_{02}^2-x_{01}^2}/k_F^{1D}, \end{aligned}$$

(18)

where $k_F^{1D}$ is the 1D Fermi wave vector, which we will keep fixed when we shrink $L\rightarrow 0$. Then the quasi-1D IBM KS orbitals are

$$\begin{aligned} \Psi _{1,k_z}(\rho )= & {} \frac{1}{\sqrt{L_z}}e^{ik_z z}\phi _1(\rho ), \nonumber \\ \phi _1(\rho )= & {} \frac{1}{L\sqrt{\pi }J_1(x_{01})}J_0(\frac{x_{01}}{L}\rho ). \end{aligned}$$

(19)

Using the rule $\sum _{k_z} \rightarrow \frac{L_z}{2\pi }\int ^{k_F^{1D}}_{-k_F^{1D}} d k_z$, we obtain the quasi-1D IBM density

$$\begin{aligned} n(\rho )=2\frac{k_F^{1D}}{\pi }\phi _1^2(\rho )=\frac{2 k_F^{1D}}{L^2\pi ^2 J_1^2(x_{01})}J_0^2(\frac{x_{01}}{L}\rho ). \end{aligned}$$

(20)

In the limit $L\rightarrow 0$, the system approaches the 1D UEG. Varying L is equivalent to performing a non-uniform density scaling in two dimensions,

$$\begin{aligned} n_\lambda ^{xy}(x,y,z)=\lambda ^2 n(\lambda x,\lambda y,z). \end{aligned}$$

(21)

with $\lambda \sim 1/L \rightarrow \infty$. Under this scaling, the reduced gradient $s=|\nabla n|/[2(3\pi ^2)^{1/3}n^{4/3}]$ behaves as $s_\lambda \sim \lambda ^{1/3}$.

Fig. 4

Upper panel: Quasi-1D IBM density versus the radial distance $\rho$ (in unit of L), for several values of L ($L=L_{\text {max}}/2$,$L=L_{\text {max}}/5$, and $L=L_{\text {max}}/10$), and for $k_F^{1D}=0.5$. Lower panel: The reduced gradient (s) and reduced Laplacian ($q=\nabla ^2 n /[4 (3\pi ^2)^{2/3}n^{5/3}]$) versus the radial distance $\rho$ for the quasi-1D IBM systems of the upper panel, s (solid line) and q (dashed line)

In the upper panel of Fig. 4, we show how the quasi-1D IBM density (with fixed $k_F^{1D}=0.5$) changes when L decreases. We consider the cases $L=L_{\text {max}}/2$,$L=L_{\text {max}}/5$, and $L=L_{\text {max}}/10$. Note that $N=\int d\rho \; 2\pi \rho n(\rho )=2 k_F^{1D}/\pi$ is dependent only on the 1D Fermi wave vector. In the lower panel of Fig. 4, we show the reduced gradient s and Laplacian q for these quasi-1D IBM systems. When $\rho \rightarrow 0$ then $s\sim - L^{-4/3}(k_F^{1D})^{-1/3}\rho + \mathcal {O}(\rho ^2)$ and $q\sim - 0.365/[k_F^{1D} L]^{2/3} + \mathcal {O}(\rho ^2)$. On the other hand, when $\rho \rightarrow L$, both s and q are diverging. Finally, we observe that q is almost constant for a large part of the quasi-1D region, and after that increases very sharply.

3.1 Kinetic energy of the quasi-1DEG

The quasi-1D IBM kinetic energy density is

$$\begin{aligned} \tau= & {} \tau ^W+\tau ^P,\nonumber \\ \tau ^W= & {} \frac{k_F^{1D}}{\pi }\Big ( \frac{d \phi _1(\rho )}{d \rho } \Big )^2,\nonumber \\ \tau ^P= & {} \frac{(k_F^{1D})^3}{3\pi }\phi _1^2(\rho ). \end{aligned}$$

(22)

The averaged kinetic energies per particle are

$$\begin{aligned} T_s/N= & {} T_s^W/N+T_s^P/N,\nonumber \\ T_s^W/N= & {} x_{01}^2/[2 L^2],\nonumber \\ T_s^P/N= & {} (k_F^{1D})^2/6=t_s^{1D}, \end{aligned}$$

(23)

where $t_s^{1D}=\tau ^{1D}/n^{1D}$ is the 1D UEG kinetic energy per particle. Then the Pauli KE per particle fully recovers the 1D UEG, while the von Weizsäcker part diverges as $\sim L^{-2}$, representing the short-wavelengths oscillations in the circular $\rho$-direction. We observe strong similarities with the quasi-2D case. The von Weizsäcker functional derivative is a constant

$$\begin{aligned} \frac{\delta T_s^W}{\delta n}=\frac{x_{01}^2}{2L^2}, \end{aligned}$$

(24)

and using the Euler equation and Eq. (14), we conclude that

$$\begin{aligned} \frac{\delta T_s^P}{\delta n}=\frac{d\tau ^{1D-UEG}}{d n^{1D-UEG}}=\pi ^2 (n^{1D})^2 /8. \end{aligned}$$

(25)

Fig. 5

Quasi-1D IBM Pauli kinetic energies per particle ($T_s^P/N$) from various KE functionals, versus $L/L_{\text {max}}$, for $k_F^{1D}=0.5$. The lower panel shows the strong quasi-1D regime ($L/L_{\text {max}}\le 0.1$)

In Fig. 5, we show the quasi-1D $T_s^P/N$ versus $L/L_{\text {max}}$, for the same KE functionals used in Fig. 1, when $k_F^{1D}=0.5$. We observe that all functionals fail badly, diverging in the 1D limit. The best functional in the strong quasi-1D regime is the PGint GGA, while PC meta-GGA and E00 GGA are relatively accurate in the moderate quasi-1D regime. However, we mention that the quasi-1D IBM is one of the most difficult tests for KE functionals.

Fig. 6

Upper panel: Pauli potential $v_s^P=\delta T_s^P/\delta n$ versus the scaled distance $\rho /L$, for the quasi-1D IBM quantum well with $L=L_{\text {max}}/2$ and $k_F^{1D}=0.5$. Lower panel: The averaged Pauli potential $\bar{v}_s^P=\int d\rho \;2\pi \rho \;n v_s^P /N$ versus $L/L_{\text {max}}$, for the 1D Fermi wave vector $k_F^{1D}=0.5$

In the upper panel of Fig. 6, we show the Pauli potential $v_s^P=\delta T_s^P/\delta n$ computed from several KE functionals, using the exact density of the quasi-1D IBM with $k_F^{1D}=0.5$ and $L=L_{\text {max}}/2$. Due to the Euler equation, the exact curve must be a constant representing the kinetic potential of the 1D UEG. Any tested KE functional can not give a constant Pauli potential. However, LKT and PG1 Pauli potentials have less structure in the region $0 \le \rho /L \le 0.6$. All semilocal KE functionals tested in Fig. 6, have $v^P\rightarrow 0$ when $\rho /L\rightarrow 1$. This feature is a consequence of their behavior in the tail of the density. Finally, we observe close similarities of this figure with the upper panel of Fig. 2.

In the lower panel of Fig. 6, we show the averaged Pauli potential $\bar{v}_s^P= \int dx\;n v_s^P /N$ versus $L/L_{\text {max}}$ for the case $k_F^{1D}=0.5$. All functionals perform similarly, in the weak and moderate quasi-1D regime (a.i. when $L/L_{\text {max}}\ge 0.4$), they give $\bar{v}_s^P$ quite close to the exact 1D UEG potential, but in the strong quasi-1D regime they fail badly, with $\bar{v}_s^P$ diverging when 1D limit is approaching. The trend is similar with the one of Fig. 5, with PGint providing the best performance.

Fig. 7

Exchange enhancement factor $F_x$ versus the reduced gradient s, for several exchange functionals. Here $\alpha =(\tau -\tau ^W)/\tau ^{\text {unif}}$

Fig. 8

Exchange energy per particle ($\epsilon _x=E_x/N$) versus $L/L_{\text {max}}$ for the quasi-1D IBM with 1D Fermi wave vector $k_F^{1D}=2$ (upper panel), and $k_F^{1D}=0.5$ (lower panel)

3.2 Exchange energy of the quasi-1DEG

The first-order density matrix of the quasi-1D IBM is

$$\begin{aligned} n_1(\mathbf {r},\mathbf {r}')=\frac{1}{\pi }\phi _1(\rho )\phi _1(\rho ') \frac{2\sin (k_F^{1D}(z'-z))}{z'-z}, \end{aligned}$$

(26)

such that $n(\rho )=n_1(\mathbf {r},\mathbf {r})$. On the other hand, the density matrix of the 1D UEG is

$$\begin{aligned} n_1^{2D-UEG}(z,z')=\frac{2}{\pi }\frac{\sin (k_F^{1D}(z'-z))}{z'-z}. \end{aligned}$$

(27)

The exchange energy is

$$\begin{aligned} E_x= & {} -\frac{2}{\pi ^3 L^2 J_1(x_{01})^4}\int _0^1 dt \int _0^1 dt' \int _0^{2\pi }d\theta \int _{-\infty }^\infty dy \; t t' \nonumber \\&\frac{\sin ^2(k_F^{1D}Ly)}{y^2 \sqrt{t^2+t'^2-2tt'\cos (\theta )+y^2}}J_0(x_{01}t)^2 J_0(x_{01}t')^2,\nonumber \\ \end{aligned}$$

(28)

where, without any loss of generality, we chose $\mathbf {r}=(\rho ,0,0)$, and we consider the changes of variables $t=\rho /L$, $t'=\rho '/L$, and $y=z'/L$. We note that the 1D UEG exchange energy per particle diverges, because of the known Coulomb divergence in 1D.

Using the non-uniform density scaling in two dimensions, see Eq. (21), we find that a GGA exchange enhancement factor (defined by $E_x=\int d\mathbf {r}\;n\epsilon _x^{LDA}F_x$) must behave in the quasi-1D regime as

$$\begin{aligned} F_x={\text {constant}}/s^2. \end{aligned}$$

(29)

This is very different from the quasi-2D case, where $F_x\propto s^{-1/2}$, see Ref. [54]. Then we propose the exchange functional, named Q1D GGA, with the following exchange enhancement factor:

$$\begin{aligned} F_x^{Q1D}=F_x^{PBEsol}+\frac{s^4+s^6}{1+s^4+s^6} \Big (-F_x^{PBEsol}+\frac{a}{s^2}\Big ),\nonumber \\ \end{aligned}$$

(30)

where $a=0.06525$ has been fitted to the quasi-1D IBM. By construction, Q1D GGA is accurate in the quasi-1D density regime, and recovers PBEsol GGA exchange functional at small reduced gradients.

In Fig. 7, we show a comparison of several exchange enhancement factors. The Q1D GGA recovers PBEsol only at very small reduced gradients ($s\le 0.5$), and after that $F_x^{Q1D}$ sharply decays as $\frac{a}{s^2}$.

In Fig. 8, we report a comparison between EXX, MVS meta-GGA [51], SCAN meta-GGA [52], Q2D GGA [53], and Q1D GGA exchange energies per particle ($\epsilon _x=E_x/N$), in the whole quasi-1D regime ($0\le L/L_{\text {max}}\le 1$). MVS and SCAN perform almost identical, and Q2D GGA is just a little better, all of them failing badly when $L\rightarrow 0$. On the other hand, Q1D GGA is remarkably accurate for the quasi-1D IBM, in both $k_F^{1D}=2$ and $k_F^{1D}=0.5$ cases. In fact, the same accuracy is obtained for any value of 1D Fermi wave vector.

4 Conclusions

The purpose of this work was to show the fundamental limitation of the 3D local/semilocal exchange-correlation energy functional approximations of DFT by considering systems with 2D characteristics. We have shown that the dimensional crossover from 3D to 2D of the exact xc energy can be significantly improved at a meta-GGA level, and we derive different exact constraints using an IBM quasi-2D. The 2D limit can be considered as a constraint on approximate functionals. For the 1D limit case we have obtained the $F_x \propto 1/s^2$ constraint with the IBM quasi-1D model and we have proposed a new functional that works well in this limit: the Q1D GGA functional.

Appendices

Appendix A: Overview of KE functionals (PGint, PG1, LKT, E00, TFW, PC)

In this work, we have considered the following kinetic energy functionals:

• (1) PGint[MGGA], whose enhancement factor has the form

$$\begin{aligned} F_s^{PGint}(s)=e^{-\mu (s)s^2}+\frac{5}{3}s^2, \quad \mu (s)=\mu _1+(\mu _2-\mu _1)\frac{\alpha s^2}{1+\alpha s^2} \end{aligned}$$

(A1)

with $\mu _1=40/27$, $\mu _2=20/9$, and the parameter $\alpha =10$ has been chosen such that PGint will be close to PG$\frac{20}{9}$ for $s\ge 0.2$.

• (2) PG1[MGGA], whose enhancement factor has the form

$$\begin{aligned} F_s^{PG1}=e^{-s^2}+\frac{5}{3}s^2. \end{aligned}$$

(A2)

• (3) LKT[MGGA], whose enhancement factor has the form

$$\begin{aligned} F_s^{LKT}=\frac{1}{\text {cosh}(as)}+ \frac{5}{3}s^2, \quad \text {with} \quad a=1.3. \end{aligned}$$

(A3)

• (4) E00[GGA], whose enhancement factor has the form

$$\begin{aligned} F_s^{E00}=\frac{135+28s^2+5s^4}{135+3s^2}. \end{aligned}$$

(A4)

• (5) TFW[GGA], whose enhancement factor has the form

$$\begin{aligned} F_s^{TFW}=1+\frac{5}{3}s^2. \end{aligned}$$

(A5)

Appendix B: Overview of XC functionals (MVS, SCAN, Q2D)

In this work, we have considered the following exchange energy functionals:

• (1) MVS[MGGA],[55] whose enhancement factor has the form

$$\begin{aligned} F_x^{MVS}(s)=\frac{1+f_x(\alpha ^{KS})k_0}{(1+b s^4)^{1/8}}, \end{aligned}$$

(B1)

where $k_0=0.174$, $b=0.0233$ and

$$\begin{aligned} f_x(\alpha )=\frac{1-\alpha }{[(1+e_1\alpha ^2)^2+c_1\alpha ^4]^{1/4}} \end{aligned}$$

(B2)

with $e_1=-1.6665$ and $c_1=0.7438$.

• (2) SCAN[MGGA],[55] whose enhancement factor has the form

$$\begin{aligned} F_x^{SCAN}=\{h_x^1+f_x(\alpha ^{KS})[h_x^0-h_x^1]\}g_x(s) \end{aligned}$$

(B3)

with

$$\begin{aligned} h_x^1=1+k_1-\frac{k_1}{1+x/k_1}, \end{aligned}$$

(B4)

$$\begin{aligned} x=\mu s^2[1+\frac{b_4s^2}{\mu }e^{-|b_4|s^2/\mu }]+ [b_1s^2+b_2(1-\alpha ^{KS})e^{-b_3(1-\alpha ^{KS})^2}]^2, \end{aligned}$$

(B5)

$$\begin{aligned} f_x(\alpha )=e^{-c_{1x}\alpha /(1-\alpha )}H(1-\alpha ) -d_xe^{c_{2x}/(1-\alpha )}h(\alpha -1), \end{aligned}$$

(B6)

$$\begin{aligned} g_x(s)=1-e^{-a_1s^{-1/2}}, \end{aligned}$$

(B7)

where H(x) is the Heaviside step function. Several of the parameters have been fixed considering exact constraints: $h_x^0=1.174$, enforces the strongly tightened bound $F_x<1.174$; $\mu =\mu _{x}^{GE2}=10/81$, $b_1=(511/13500)/(2b_2)=0.1566$, $b_2=(5913/405000)^{1/2}=0.1208$, $b_3=0.5$, and $b_4=\mu ^2/k_1-1606/18225-b_1^2=0.1218$ are fixed from the GE4 behavior; $a_1=4.9479$, is a norm related to the exact-exchange energy of the hydrogen atom. The other parameters are $c_{1x}=0.667, c_{2x}=0.8, d_{x}=1.24$, and $k_1=0.065$.

• (3) Q2D[GGA], [56] whose enhancement factor has the form

$$\begin{aligned} F_x^{Q2D}(s)=\frac{F_x^{PBEsol}(s)(c-s^4)+0.5217s^{7/2}(1+s^2)}{c+s^6} \end{aligned}$$

(B8)

with $c=10^2$ that are fixed by fitting to the IBM reference system.