Reduced Density Matrix Functional Theory for Superconductors

TL;DR

This paper introduces a novel er initio reduced density matrix functional theory for superconductors, establishing a universal functional and deriving Kohn-Sham equations to improve upon existing superconductivity models.

Contribution

It develops a new formalism linking the density operator to reduced densities, deriving a universal functional and Kohn-Sham equations for superconductors at finite temperature.

Findings

Existence of a universal functional _eta[,] for superconductor ground states.

Derivation of Bogoliubov-de Gennes-like equations within the new formalism.

Extension of the Sham-Schlfcter connection to develop an exchange-correlation functional.

Abstract

We present an \textit{ab initio} theory for superconductors, based on a unique mapping between the statistical density operator at equilibrium, on the one hand, and the corresponding one-body reduced density matrix $\gamma$ and the anomalous density $\chi$, on the other. This new formalism for superconductivity yields the existence of a universal functional $\mathfrak{F}_\beta[\gamma,\chi]$ for the superconductor ground state, whose unique properties we derive. We then prove the existence of a Kohn-Sham system at finite temperature and derive the corresponding Bogoliubov-de Gennes-like single particle equations. By adapting the decoupling approximation from density functional theory for superconductors we bring these equations into a computationally feasible form. Finally, we use the existence of the Kohn-Sham system to extend the Sham-Schl\"uter connection and derive a first exchange-correlation functional for our theory. This reduced density matrix functional theory for superconductors has the potential of overcoming some of the shortcomings and fundamental limitations of density functional theory of superconductivity.

Equations

$$(n,\chi )\stackrel{1-1}{⟷}{\hat{\varrho }}_{\mathrm{eq}}\equiv \frac{e^{-\beta ({\hat{H}}_{v,\Delta }-\mu \hat{N})}}{\text{Tr}[e^{-\beta ({\hat{H}}_{v,\Delta }-\mu \hat{N})}]},$$

$$|{\Psi }_{\mathrm{gs}}⟩\stackrel{1-1}{⟷}{\hat{\gamma }}_{\mathrm{gs}}\equiv N{\mathrm{Tr}}_{N-1}[|{\Psi }_{\mathrm{gs}}⟩⟨{\Psi }_{\mathrm{gs}}|],$$

$$\mathcal{E}[\gamma ]=\int {\mathrm{d}}^{4}x\int {\mathrm{d}}^4{x}^{\prime }h(\mathbf{r},{\mathbf{r}}^{\prime })\gamma (\mathbf{x};{\mathbf{x}}^{\prime })+\mathcal{F}[\gamma ],$$

$${\hat{\gamma }}_2\equiv \left(\genfrac{}{}{0ex}{}{N}{2}\right){\mathrm{Tr}}_{N-2}[|\Psi ⟩⟨\Psi |].$$

$$\displaystyle\gamma_{2}(\bm{x};\bm{x}^{\prime})\equiv\tfrac{1}{2}\gamma(\bm{x};\bm{x})\big{[}\gamma(\bm{x}^{\prime};\bm{x}^{\prime})-\rho^{\rm hole}_{\rm xc}(\bm{x},\bm{x}^{\prime})\big{]}.$$

$${\left|\frac{{\gamma }^{1/2}(\mathbf{x};{\mathbf{x}}^{\prime })}{\sqrt{\gamma (\mathbf{x};\mathbf{x})}}\right|}^2,$$

$${\Omega }_{\mu ,\beta }\left[\varrho \right]=\text{Tr}\hat{\varrho }\left(\hat{H}-\mu \hat{N}+\frac{1}{\beta }\ln \hat{\varrho }\right).$$

$${\hat{v}}_{\mathrm{ext}}-\mu \stackrel{1-1}{⟷}{\hat{\varrho }}_{\mathrm{eq}}\stackrel{1-1}{⟷}{\hat{\gamma }}_{\mathrm{eq}}\equiv N{\mathrm{Tr}}_{N-1}[{\hat{\varrho }}_{\mathrm{eq}}].$$

$${\Omega }_{\mu ,\beta }\left[{\hat{\gamma }}_{\mathrm{eq}}\right]=\text{Tr}{\hat{\varrho }}_{\mathrm{eq}}[{\hat{\gamma }}_{\mathrm{eq}}]\left(\hat{H}-\mu \hat{N}+\frac{1}{\beta }\ln {\hat{\varrho }}_{\mathrm{eq}}[{\hat{\gamma }}_{\mathrm{eq}}]\right).$$

$$\displaystyle\bar{\Psi}_{\sigma}(\bm{r})=\bigg{(}\begin{array}[]{c}\psi^{\dagger}_{\sigma}(\bm{r})\\ \psi_{-\sigma}(\bm{r})\end{array}\bigg{)}.$$

$${\Gamma }_{\sigma \sigma }(\mathbf{r},{\mathbf{r}}^{\prime })$$

$$=\left(\begin{array}{cc}{\gamma }_{\sigma \sigma }(\mathbf{r},{\mathbf{r}}^{\prime })& {\chi }^{†}(\mathbf{r},{\mathbf{r}}^{\prime })\\ \chi (\mathbf{r},{\mathbf{r}}^{\prime })& \delta (\mathbf{r}-{\mathbf{r}}^{\prime })-{\gamma }_{-\sigma -\sigma }(\mathbf{r},{\mathbf{r}}^{\prime })\end{array}\right).$$

$$\mathbf{\Gamma }(\mathbf{r},{\mathbf{r}}^{\prime })=\left(\begin{array}{cc}{\Gamma }_{\uparrow \uparrow }(\mathbf{r},{\mathbf{r}}^{\prime })& {\Gamma }_{\uparrow \downarrow }(\mathbf{r},{\mathbf{r}}^{\prime })\\ {\Gamma }_{\downarrow \uparrow }(\mathbf{r},{\mathbf{r}}^{\prime })& {\Gamma }_{\downarrow \downarrow }(\mathbf{r},{\mathbf{r}}^{\prime })\end{array}\right)$$

$${\hat{H}}_{v,\Delta }$$

$$+\frac{1}{2}\sum _{\sigma {\sigma }^{\prime }}\int {\mathrm{d}}^{3}r\int {\mathrm{d}}^3{r}^{\prime }{\hat{\psi }}_{\sigma }^{†}(\mathbf{r}){\hat{\psi }}_{{\sigma }^{\prime }}^{†}({\mathbf{r}}^{\prime })\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\hat{\psi }}_{{\sigma }^{\prime }}({\mathbf{r}}^{\prime }){\hat{\psi }}_{\sigma }(\mathbf{r})$$

$$-\int {\mathrm{d}}^3{r}_1\int {\mathrm{d}}^3{r}_1^{\prime }\int {\mathrm{d}}^3{r}_2\int {\mathrm{d}}^3{r}_2^{\prime }{\hat{\psi }}_{\downarrow }^{†}({\mathbf{r}}_1^{\prime }){\hat{\psi }}_{\uparrow }^{†}({\mathbf{r}}_1)u({\mathbf{r}}_1^{\prime },{\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_2^{\prime }){\hat{\psi }}_{\uparrow }({\mathbf{r}}_2){\hat{\psi }}_{\downarrow }({\mathbf{r}}_2^{\prime })$$

$$\equiv \hat{T}+\hat{v}-(\hat{\Delta }+\text{H.c.})+\hat{W}+\hat{U}.$$

$${\Omega }_{v,\Delta }\left[\varrho \right]=\text{Tr}\hat{\varrho }\left({\hat{H}}_{v,\Delta }-\mu \hat{N}+\frac{1}{\beta }\ln \hat{\varrho }\right).$$

$${\hat{\varrho }}_{\mathrm{eq}}=\frac{1}{\mathcal{Z}}{e}^{-\beta \left({\hat{H}}_{v,\Delta }-\mu \hat{N}\right)},$$

$${\varrho }_{\mathrm{eq}}\stackrel{1-1}{⟷}\left(v,\Delta \right).$$

$$\frac{e^{-\beta ({\hat{H}}_{v,\Delta }-\mu \hat{N})}}{\mathcal{Z}}=\frac{e^{-\beta ({\hat{H}}_{v^{\prime },{\Delta }^{\prime }}-\mu \hat{N})}}{{\mathcal{Z}}^{\prime }}.$$

$${\hat{H}}_{v,\Delta }-{\hat{H}}_{v^{\prime },{\Delta }^{\prime }}=\hat{v}-{\hat{v}}^{\prime }-(\hat{\Delta }-{\hat{\Delta }}^{\prime }+h.c.)=\frac{1}{\beta }\ln \frac{{\mathcal{Z}}^{\prime }}{\mathcal{Z}}.$$

$${\Omega }_{v^{\prime },{\Delta }^{\prime }}\left[{\varrho }^{\prime }\right]$$

$$<\text{Tr}\hat{\varrho }\left({\hat{H}}_{v^{\prime },{\Delta }^{\prime }}-\mu \hat{N}+\frac{1}{\beta }\ln \hat{\varrho }\right)$$

$$={\Omega }_{v,\Delta }\left[\varrho \right]+\text{Tr}[\hat{\varrho }({\hat{v}}^{\prime }-\hat{v})-\hat{\varrho }({\hat{\Delta }}^{\prime }-\hat{\Delta }+h.c.)]$$

$$={\Omega }_{v,\Delta }\left[\varrho \right]+v^{\prime }[\gamma ]-v[\gamma ]-{\Delta }^{\prime }\left[\chi \right]+\Delta \left[\chi \right],$$

$$v[\gamma ]$$

$$\Delta \left[\chi \right]$$

$${\Omega }_{v,\Delta }\left[\varrho \right]<{\Omega }_{v^{\prime },{\Delta }^{\prime }}\left[{\varrho }^{\prime }\right]-v^{\prime }[\gamma ]+v[\gamma ]+{\Delta }^{\prime }\left[\chi \right]-\Delta \left[\chi \right].$$

$${\Omega }_{v,\Delta }\left[\varrho \right]+{\Omega }_{v^{\prime },{\Delta }^{\prime }}\left[{\varrho }^{\prime }\right]<{\Omega }_{v^{\prime },{\Delta }^{\prime }}\left[{\varrho }^{\prime }\right]+{\Omega }_{v,\Delta }\left[\varrho \right],$$

$${\Omega }_{v,\Delta }^{\mathrm{eq}}$$

$$=\text{Tr}[(\hat{T}+\hat{v}-\mu ){\hat{\gamma }}_{\mathrm{eq}}]-\Delta [{\chi }_{\mathrm{eq}}]+{\mathcal{F}}_{\beta }\left[{\gamma }_{\mathrm{eq}},{\chi }_{\mathrm{eq}}\right].$$

$${\mathcal{F}}_{\beta }\left[{\gamma }_{\mathrm{eq}},{\chi }_{\mathrm{eq}}\right]$$

$$\times \left(\hat{W}+\hat{U}+\frac{1}{\beta }\ln {\hat{\varrho }}_{\mathrm{eq}}\left[{\gamma }_{\mathrm{eq}},{\chi }_{\mathrm{eq}}\right]\right),$$

Full text

I Introduction

Superconductivity was discovered more than one century ago, when Kamerlingh Onnes observed, and wrote in his lab notes, that “mercury’s resistance is practically zero” at very low temperatures van Delft and Kes (2010). While the first phenomenological theory of superconductivity was developed 24 years later in the form of the London equations London and London (1935); de Gennes (1999), it took nearly half a century until the first successful microscopic description emerged Bardeen et al. (1957). The so-called Bardeen-Cooper-Schrieffer (BCS) theory, based on the concept of Cooper pairs and a non-particle number conserving Ansatz for the wave function, predicted a pair-bound-state in the presence of an attractive force, regardless of the strength of the interaction. Although BCS theory is able to explain some universal features of superconductors (such as the ratio between the energy gap and the critical temperature), a proper treatment of the strong electron-phonon coupling regime had to wait until the theory developed by Eliashberg Eliashberg (1960).

It is by and large accepted that the electron-phonon interaction is well accounted for within BCS and Eliashberg theories. Yet the correlation effects due to the electron-electron Coulomb repulsion are extraordinarily difficult to handle, and they are usually compressed in an adjustable parameter (namely, $\mu^{*}$). Since the adjustability of the parameter diminishes the predictive capability of the theory, it comes as no surprise that the correct description of superconductors remains one of the great questions in condensed matter theory. We still have a poor understanding of the superconducting features of unconventional superconductors, like high-$T_{c}$ cuprates Bednorz and Müller (1986); Lee et al. (2006); Berthod et al. (2017) or layered organic materials Powell and McKenzie (2011). Moreover, recent discoveries of superconductivity in two layers of graphene Cao et al. (2018); Xu and Balents (2018), and in compressed hydrides at extreme pressure Drozdov et al. (2015); Errea et al. (2015); Somayazulu et al. (2019), continue to challenge our understanding of superconducting systems.

Proposed by Oliveira, Gross and Kohn in 1988, density functional theory for the superconducting state (SC-DFT) provided an ab initio unified treatment of correlation and inhomogeneity effects in superconductors Oliveira et al. (1988). Extending the famous Hohenberg-Kohn theorems of normal density functional theory (DFT) Hohenberg and Kohn (1964), SC-DFT states a one-to-one mapping between the equilibrium statistical density operator $\hat{\varrho}_{\rm eq}$ and the corresponding electronic $n(\bm{r})\equiv\sum_{\sigma}\langle\psi_{\sigma}^{\dagger}(\bm{r})\psi_{\sigma}(\bm{r})\rangle$ and anomalous $\chi(\bm{r},\bm{r}^{\prime})\equiv\langle\psi_{\uparrow}(\bm{r})\psi_{\downarrow}(\bm{r}^{\prime})\rangle$ densities:

[TABLE]

with $\hat{H}_{v,\Delta}$ being the Hamiltonian for a superconductor in an external potential $v(\bm{r})$ and a non-local pairing potential $\Delta(\bm{r},\bm{r}^{\prime})$. This latter term is included to break the $U(1)$ symmetry in the superconducting phase. If not otherwise indicated, “$\Tr$” always means the trace in Fock space. Since SC-DFT accounts for electron-electron interactions in the same way as normal DFT, there is no use of any semi-empirical parameter. Hence, SC-DFT is a truly ab initio theory for superconductors.

Formally, SC-DFT is able to describe superconductivity in all systems, for the many-electron problem is cast into a universal exchange-correlation functional whose existence is ensured by the Oliveira-Gross-Kohn theorem. Further work on this topic Lüders et al. (2005); Marques et al. (2005) extended SC-DFT to include the nuclear density as a third component. Recently, the theory was revisited to allow the inclusion of magnetic fields Linscheid et al. (2015a, b).

SC-DFT has been extremely successful in predicting superconductivity in a wide variety of materials Sanna et al. (2007); Flores-Livas et al. (2012); Floris et al. (2007), and in particular in high-pressure research Profeta et al. (2006); Akashi et al. (2015); Flores-Livas et al. (2016). Yet, the practical applicability of SC-DFT rests (as in the standard DFT) upon a Kohn-Sham scheme, which maps the real interacting system of interest to an auxiliary noninteracting one with the same equilibrium electronic and anomalous densities $(n,\chi)$. The existence of such a Kohn-Sham system is problematic, however. In fact, in a recent paper we proved that such a noninteracting system does not exist, at least at zero temperature Schmidt et al. (2019). Furthermore, there is an unpleasant asymmetry in SC-DFT that complicates considerably certain key points in the derivations. This asymmetry is related to the use of the local electron density [there is only one $\bm{r}$ in $n(\bm{r})$] and the non-local anomalous density [there are two different $\bm{r}$ and $\bm{r}^{\prime}$ in $\chi(\bm{r},\bm{r}^{\prime})$]. Although no problem arises at the level of the one-to-one correspondence of Eq. (1), we can expect complications when developing exchange-correlation functionals for SC-DFT. Indeed, for superconductivity only a few functionals of the electronic and anomalous densities have been developed. Even worse, these functionals require completely ad hoc, and rather arbitrary, modifications in order to yield good results Mazin et al. (2008).

It is a well-know fact that the most spectacular failures of normal DFT (with standard exchange-correlation functionals) are related to an incorrect description of strongly (static, in quantum-chemistry jargon) correlated systems Cohen et al. (2008a, b). Since Kohn-Sham DFT does not deal with fractional occupation numbers, the electronic structure of multireference systems is still an open problem within the theory. Reduced density matrix functional theory (RDMFT) is a natural extension of DFT, whose aim is to exploit the one-electron picture of the many-body wave function $\ket{\Psi}$ by seeking a functional of the one-electron reduced density matrix $\hat{\gamma}=N{\Tr}_{N-1}[\ket{\Psi}\bra{\Psi}]$, allowing therefore fractional occupation numbers Gilbert (1975); Schade et al. (2017). RDMFT describes closed-shell molecular systems with accuracies higher by one order of magnitude than DFT Lathiotakis and Marques (2008a); Piris (2017a). It has also succeeded in predicting more accurate gaps of conventional semiconductors Lathiotakis et al. (2009). Furthermore, RDMFT correctly captures the physics of the insulator-metal phase transition of transition metal oxides Sharma et al. (2013), something, as is well known, DFT cannot do.

In this article, we develop a new ab initio theory of superconductivity, namely, a RDMFT formalism for superconductivity (from now on, SC-RDMFT). This is done by using $\gamma$ instead of $n$ as one of the basic variables of the theory. In this way, we expect that most of the advantages of RDMFT will translate well to strongly correlated superconducting systems. Furthermore, the use of the density-matrix restores the symmetry of the equations, and as we will see, simplifies considerably some derivations.

We then derive the existence of a Kohn-Sham system which results in a set of coupled Bogoliubov-de-Gennes-like equations. These equations are still extremely difficult to solve as they require precision on the level of the superconductive coupling while encompassing the energy scale of the whole band structure Zhu (2016). In order to simplify the Kohn-Sham equations into a form that can feasibly be solved we appropriate the so-called decoupling approximation from SC-DFT. It allows us to separate the solution of the standard electronic problem, the influence of the lattice and the superconducting coupling. Analogously to DFT, the challenges of the many-body problem cannot be avoided through a Kohn-Sham system but only be transferred to the exchange-correlation functional. We suggest a first solution to this challenge by introducing the Sham-Schlüter connection Godby et al. (1986, 1987a, 1987b) to SC-RDMFT and deriving a first functional. Moreover, the derivation avoids several ad-hoc approximations that are required in SC-DFT.

The paper is organized as follows. The first section is this introduction. For the paper’s completeness, Sec. II summarizes normal RDMFT and its advantages in the light of normal DFT. In Sec. III we present the generalized one-body reduce density matrix in Nambu-Gorkov space. Such a matrix is the correct variable to deal with the spontaneous breaking of the $U(1)$ symmetry in the superconducting phase. Section IV presents the theoretical foundations of SC-RDMFT. We derive a Gilbert theorem and state the mathematical properties of the universal exchange-correlation functional for superconductors. In Sec. V we prove the existence of the Kohn-Sham system at finite temperature, followed by the decoupling approximation in Sec. VI and the Sham-Schlüter connection in Sec. VII. Section VIII is devoted to the conclusions of the paper.

II Reduced density functional theory in a nutshell

As a solution of the Schrödinger equation, the state $\ket{\Psi_{\rm gs}}$ describes the ground-state quantum system completely. The Hohenberg-Kohn theorems imply the existence of a universal functional of the external potential and the electronic density which reaches a minimum when evaluated on the ground-state electronic density. In 1975, Gilbert proved an extension of such theorems, showing that there is also a one-to-one correspondence between the ground-state wave function of a non-degenerate many-body system and the corresponding one-body reduced density matrix Gilbert (1975):

[TABLE]

where $N-1$ particles are traced out. The main advantage of this theorem is that any observable of the system (in its ground state) can be written as a functional of $\hat{\gamma}_{\rm gs}$. In practice, this latter condition is relaxed and the functional is evaluated on the set of fermionic one-body reduced density matrices $\{\hat{\gamma}\}$. More importantly, the functional of the kinetic energy, which is unknown in DFT and should be included in the exchange-correlation functional, is known exactly in terms of $\hat{\gamma}$. In RDMFT, the functional for the ground-state energy reads:

[TABLE]

where $\gamma(\bm{x};\bm{x}^{\prime})=\bra{\bm{x}^{\prime}}\hat{\gamma}\ket{\bm{x}}$. The one-particle Hamiltonian $h(\bm{r},\bm{r}^{\prime})=-\tfrac{1}{2}\delta(\bm{r}-\bm{r}^{\prime})\nabla^{2}_{\bm{r}}+v_{\rm ext}(\bm{r},\bm{r}^{\prime})$ contains the kinetic and external-potential operators. $\mathcal{F}[\gamma]$ is an unknown universal functional of the ground-state one-body reduced density matrix. We used the customary compact notation for spin and position coordinates $\bm{x}\equiv(\bm{r},\sigma)$ (with $\sigma\in\{\uparrow,\downarrow\})$.

Another important advantage of Gilbert’s theorem is that one can consider a broader set of nonlocal external potentials $v_{\rm ext}(\bm{r},\bm{r}^{\prime})$, which arise when the quantum problem is formulated in terms of many valence electrons, subject to an external potential of fixed nuclei and core electrons Gilbert (1975). Yet, the exact form of the exchange-correlation functional $\mathcal{F}[\gamma]$ is, by and large, not available and therefore the predicted RDMFT energy, with approximate functionals, can be either lower or higher than the exact ground-state energy Müller (1984); Frank et al. (2007); Benavides-Riveros and Várilly (2012); Mitxelena et al. (2017). Functionals in RDMFT are often engineered as approximate expressions of the two-body reduced density matrix

[TABLE]

This is normally accomplished by writing $\hat{\gamma}_{2}$ in terms of the exchange-correlation hole $\rho^{\rm hole}_{\rm xc}(\bm{x},\bm{x}^{\prime})$, defined by the following relation:

[TABLE]

The electronic density, the main object in DFT, is of course the trace on spin space of the diagonal $\gamma(\bm{x},\bm{x})$.

For example, the famous Müller (for historical reasons also called Buijse-Baerends) functional describes the exchange-correlation hole as the square of a hole amplitude Buijse and Baerends (2002), reading

[TABLE]

where $\hat{\gamma}^{1/2}\equiv\sum_{i}n^{1/2}_{i}\ket{\varphi_{i}}\bra{\varphi_{i}}$ is written in terms of the so-called natural occupation numbers $\{n_{i}\}$ and natural orbitals $\{\varphi_{i}\}$, namely, the eigenvalues and eigenvectors of $\hat{\gamma}$. Further developments in RDMFT are inspired by this seminal functional Goedecker and Umrigar (1998); Cioslowski and Pernal (1999); Gritsenko et al. (2005); Csányi et al. (2002); López-Sandoval and Pastor (2002). There is another perspective by studying the cumulant part of $\gamma_{2}$ (i.e., $\gamma_{2}-\tfrac{1}{2}\gamma\wedge\gamma$) under some of its known representability conditions Piris (2007, 2017b). Remarkably, almost all RDMFT functionals fare quite well in benchmarking tests, yielding errors for the correlation energy an order of magnitude smaller than B3LYP, perhaps the most popular DFT functional in quantum chemistry, and a precision comparable to Møller-Plesset second-order perturbation theory Lathiotakis and Marques (2008b). RDMFT has also succeeded in predicting more accurate gaps of conventional semiconductors than semi-local DFT does. Furthermore it has demonstrated insulating behavior for Mott-type insulators Sharma et al. (2013); Shinohara et al. (2015); Pernal (2015). Recently, the role of the generalized Pauli exclusion principle has been stressed within RDMFT for pure states Benavides-Riveros (2018); Benavides-Riveros et al. (2017); Tennie et al. (2017); Schilling (2018); Schilling et al. (2017); Benavides-Riveros and Marques (2018); Schilling and Schilling (2019); Benavides-Riveros and Marques (2019).

When researching superconductors we naturally want to consider their behaviour at finite temperature as their most important property is their transition temperature. Fortunately, a finite-temperature RDMFT was already developed by Baldsiefen et al. Baldsiefen et al. (2015, 2017). Just as in finite-temperature DFT, the considered systems are grand canonical ensembles and instead of the energy, for a Hamiltonian $\hat{H}$, we consider the grand canonical potential $\Omega_{\mu,\beta}[\varrho]$ for fixed temperature $\beta^{-1}$ and chemical potential $\mu$:

[TABLE]

The appropriate space for the statistical density operators is the Fock space. Instead of the ground state, we are interested in the thermodynamic equilibrium. For general quantum systems, $\Omega_{\mu,\beta}\left[\varrho\right]$ is bounded from below by $\Omega_{\mu,\beta}\left[\varrho_{\rm eq}\right]$, where $\hat{\varrho}_{\rm eq}$ is the equilibrium statistical density operator. In contrast to zero-temperature RDMFT, Gilbert’s theorem in finite-temperature RDMFT is completely invertible: there is a one-to-one mapping between the equilibrium statistical density operator, the external potential (minus the chemical potential) and the corresponding equilibrium one-body reduced density operator Baldsiefen et al. (2015); Eschrig (2010):

[TABLE]

The proof is analogue to Mermin’s proof for finite-temperature DFT Mermin (1965). This result implies that $\hat{\varrho}_{\rm eq}$ can be written as a functional of $\hat{\gamma}_{\rm eq}$, and therefore, in the thermodynamic equilibrium, the functional (7) can be written as a functional of the equilibrium one-body reduced density matrix $\hat{\gamma}_{\rm eq}$, namely:

[TABLE]

It is worth saying that Eq. (8) is an advantage over the standard zero-temperature RDMFT. Indeed, there is no problem with degeneracy since degenerate states get the same equilibrium statistical density operator Eschrig (2010).

III Nambu-Gorkov one-body reduced density matrix

As presented in Section II, RDMFT is unfortunately not sufficient for the description of superconducting systems. The reason for this lies in the spontaneous breaking of the $U(1)$ symmetry, which implies that the particle number is not conserved. The breaking of such a symmetry also results in a finite expectation value of the anomalous density $\chi(\bm{r},\bm{r}^{\prime})=\langle\psi_{\uparrow}(\bm{r})\psi_{\downarrow}(\bm{r}^{\prime})\rangle$. In Nambu-Gorkov space, field operators are usually written as the following spinors:

[TABLE]

We write the generalized Nambu-Gorkov one-body reduced density matrix as the expected value of a product of Nambu-Gorkov field-operators, namely:

[TABLE]

Notice that $\Gamma_{\sigma\sigma}(\bm{r},\bm{r}^{\prime})$ can be denoted in terms of the reduced density matrix $\gamma_{\sigma\sigma}(\bm{r},\bm{r}^{\prime})$ and the anomalous density $\chi(\bm{r},\bm{r}^{\prime})$ as:

[TABLE]

Note also that the Nambu-Gorkov one-body reduced density matrix is not properly normalized. It integrates to the volume, instead of integrating to the particle number. However, one can still obtain the average particle number from the ordinary one-body reduced density matrix $\gamma$ and the standard deviation of the particle number from the anomalous density $\chi$.

Obviously this is only one of several possible choices. For instance, in order to describe a system that is not symmetric with respect to spin direction one would have to consider the full Nambu-Gorkov one-body reduced density matrix in spin space:

[TABLE]

In this notation, the anomalous densities with parallel spin describe triplet Cooper pairs (experimentally verified eight years ago by Sprungmann et al. Sprungmann et al. (2010)). For this paper we will limit ourselves to systems with spin-rotational symmetry. Therefore, as a matter of simplification, we will omit the index $\sigma\sigma$ from $\Gamma_{\sigma\sigma}$ from now on.

IV SC-RDMFT: Theoretical foundations

This section is devoted to prove the fundamental theorems of our new SC-RDMFT, namely, a Gilbert theorem for systems with spontaneous breaking of the $U(1)$ symmetry. We formulate the variational principle in order to include ensemble $N$-representable $\Gamma$ in a formalism à la Levy-Lieb. We then prove several properties of the universal functional that turn out to be crucial for minimization procedures and the development of SC-RDMFT functionals.

IV.1 Gilbert theorem for superconducting systems

We will develop RDMFT for superconductors analogously to SC-DFT and accordingly we start with the Hamiltonian for superconductors extended to non-local potentials:

[TABLE]

The terms $\hat{\Delta}$ and $\hat{\Delta}^{\dagger}$ have no analog in the non-superconducting electron gas. As in BCS theory, a non-zero expectation value for $\langle\hat{\Delta}\rangle$ must be understood as a consequence of an external source field which is included by hand to break the symmetry. It is switched off after the thermodynamic limit is taken. In Eq. (16), $\hat{T}$ denotes the kinetic energy operator, $\hat{v}$ the external potential, $\hat{W}$ the Coulomb repulsion, while the operator $\hat{U}$ is an effective electron-electron phonon-mediated attraction.

We consider systems described by such a Hamiltonian in a grand-canonical ensemble at finite temperature $T=1/(k_{\text{B}}\beta)$, where $k_{\text{B}}$ is the Boltzmann constant. The system is coupled to a particle and a heat bath through the Lagrange multipliers $\mu$ and $\beta$ (which we keep constant across this article). This results in the well known formula for the grand potential:

[TABLE]

The reader should keep in mind that $\Omega_{v,\Delta}\left[\varrho\right]$ is dressed in the symbols $\beta$ and $\mu$. We refrain from writing them in order to alleviate the notation. The minimizing equilibrium statistical density operator can be denoted as:

[TABLE]

where {\mathcal{Z}}={\Tr}\big{[}e^{-\beta\left(\hat{H}_{v,\Delta}-\mu\hat{N}\right)}\big{]}. Our objective is to prove (for fixed $\beta$ and $\mu$) a one-to-one relationship between the equilibrium statistical density operator $\varrho_{\rm eq}$, the external potentials $(v,\Delta)$, and the equilibrium Nambu-Gorkov one-body reduced density matrix $\Gamma_{\rm eq}$ or, equivalently, the pair $(\gamma_{\rm eq},\chi_{\rm eq})$. This mapping allows us to write all equilibrium observables as a functional of the equilibrium Nambu-Gorkov one-body reduced density matrix $\Gamma_{\rm eq}$.

Before beginning the proof, we would like to comment on the problem of the phonons. In order to introduce consistently the electron-phonon interaction in SC-RDMFT we can follow the same recipe as in Ref. Lüders et al. (2005); Marques et al. (2005). This involves treating the electronic and nuclear degrees of freedom at the same level, by introducing the diagonal of the $N$-particle nuclear density-matrix as an extra variable, and by developing a multi-component RDMFT Kreibich and Gross (2001). The derivation is straightforward, and follows closely the one for SC-DFT, however complicates considerably the formulas. We therefore decided to overlook the nuclei, and only to present, in Sect. VI.1, the electron-phonon term that stems from that formulation.

As a first step we prove (for fixed $\beta$ and $\mu$) the one-to-one mapping between the equilibrium statistical density operator and the pair of external potentials, namely:

[TABLE]

Thus, the density $\varrho_{\rm eq}$ is uniquely determined by the potentials $v$ and $\Delta$. The proof goes by contradiction and is adapted from the original proof for temperature-dependent RDMFT Baldsiefen et al. (2015).

We consider two Hamiltonians $\hat{H}_{v,\Delta}$ and $\hat{H}_{v^{\prime},\Delta^{\prime}}$ with the potentials $(v,\Delta)$ and $(v^{\prime},\Delta^{\prime})$ differing for more than a constant and assume that they lead to the same equilibrium statistical density operator. The following condition then holds:

[TABLE]

Here $\mathcal{Z}$ and $\mathcal{Z}^{\prime}$ are the partition functions corresponding to the Hamiltonians $\hat{H}_{v,\Delta}$ and $\hat{H}_{v^{\prime},\Delta^{\prime}}$ respectively. Solving Eq. (20) for the difference of the Hamiltonians leads to:

[TABLE]

Since the right hand of this equation is particle conserving, we obtain trivially that $\hat{\Delta}=\hat{\Delta}^{\prime}$. On the other hand, $\hat{v}$ and $\hat{v}^{\prime}$ differ by a constant, which contradicts our initial assumption. Therefore, we obtain the desired one-to-one relationship (19).

The second part of the theorem will follow the logic of the classic Mermin’s proof for finite-temperature DFT Mermin (1965). It will yield (again) the one-to-one correspondence between the equilibrium density operator $\varrho_{\rm eq}$ and $(\gamma_{\rm eq},\chi_{\rm eq})$. The proof goes (again) by reductio ad absurdum as follows. Suppose that there exist two equilibrium density operators $\varrho$ and $\varrho^{\prime}$ that amount to the same $\gamma$ and $\chi$. The Hamiltonian, statistical density operator, and grand potential associated with $v^{\prime}$ and $\Delta^{\prime}$ are labeled with $H_{v^{\prime},\Delta^{\prime}}$, $\varrho^{\prime}$ and $\Omega_{v^{\prime},\Delta^{\prime}}$. It follows that:

[TABLE]

where

[TABLE]

Interchanging the primed and unprimed variables yields the following equation:

[TABLE]

In (22) we have used the hypothesis that both $\varrho$ and $\varrho^{\prime}$ lead to the same reduced density matrices $\gamma$ and $\chi$. Adding Eq. (21) and (22) one obtains a contradiction à la Hogenberg-Kohn, namely:

[TABLE]

which is impossible. Therefore, the equilibrium reduced densities $(\gamma_{\rm eq},\chi_{\rm eq})$ are unique for every equilibrium statistical density operator $\varrho_{\rm eq}$ (for fixed $\beta$ and $\mu$). Hence, there exists $\varrho_{\rm eq}[\gamma_{\rm eq},\chi_{\rm eq}]$, a universal functional of the reduced densities.

Analogously to finite-temperature RDMFT we can write the potential $\Omega_{v,\Delta}\left[\varrho_{\rm eq}\right]$ as a functional of $(\gamma_{\rm eq},\chi_{\rm eq})$, namely,

[TABLE]

Here we define the universal functional $\mathcal{F}_{\beta}\left[\gamma_{\rm eq},\chi_{\rm eq}\right]$ whose existence is guaranteed by our two theorems. It reads

[TABLE]

where $\hat{\varrho}_{\rm eq}[\gamma_{\rm eq},\chi_{\rm eq}]$ is also a universal functional of the reduced densities. We should keep in mind that the functional $\mathcal{F}_{\beta}$ depends on $\beta$. Notice that, due to Eq. (14), the functional $\Omega_{v,\Delta}\left[\gamma_{\rm eq},\chi_{\rm eq}\right]$ can be written alternatively as a functional of the generalized one-body reduced density matrix (at equilibrium): $\Omega_{v,\Delta}\left[\Gamma_{\rm eq}\right]$.

Notice that up to this point our results are only defined on $\mathfrak{G}^{v,\Delta}$, the domain of equilibrium $(v,\Delta)$-representable Nambu-Gorkov one-body reduced density matrices, defined as:

[TABLE]

with $\hat{\Gamma}=\bar{\Psi}_{\sigma}(\bm{r})\otimes\bar{\Psi}^{\dagger}_{\sigma}(\bm{r}^{\prime}).$ While this restriction is unproblematic for now, it will be a cause of concern once we start any actual numerical minimization. The issue lies in the fact that we are not able to characterize in general the elements of $\mathfrak{G}^{v,\Delta}$. In order to overcome this challenge we will expand the minimization process to all ensemble-$N$-representable one-body reduced density matrices $\mathfrak{G}^{N}$, namely, the set of generalized one-body reduced density matrices which are produced by a physical statistical density operator $\varrho\in\mathfrak{D}^{N}$ where

[TABLE]

In a fashion similar to Levy-Lieb’s formulation of DFT Levy (1979); Lieb (1983), we can restate the functional as a variational principle. Indeed, since by definition

[TABLE]

we can divide the minimization as follows:

[TABLE]

Now the functional is minimized on the set $\mathfrak{G}^{N}$, which can easily be characterized by the condition $0\leq\Gamma\leq 1$ Bach et al. (1994); Hainzl and Seiringer (2016).

IV.2 Properties of the universal functional

The infimum in Eq. (29) is justified by the fact that for a Nambu-Gorkov one-body reduced density matrix $\Gamma\in\mathfrak{G}^{N}$ that is non-$(v,\Delta)$-representable (i.e., it comes from a non-equilibrium statistical density operator), it is not clear if there exists a minimizing $\varrho\rightarrow\Gamma$. We therefore prove the important theorem:

Theorem 1**.**

For all $(\gamma,\chi)\in\mathfrak{G}^{N}$, there exists $\varrho\in\mathfrak{D}^{N}$, such that $\varrho\rightarrow(\gamma,\chi)$ by contraction, and

[TABLE]

For the proof we require a sequence $\left\{\varrho_{k}\right\}$ with $\varrho_{k}\in\mathfrak{D}^{N}$ and $\varrho_{k}\rightarrow(\gamma,\chi)$. Just as in Ref. Baldsiefen et al. (2015) the proof is divided in three parts:

I. $\exists\,\varrho:\varrho_{k}\rightarrow\varrho$* on the weak-* topology.*

The proof of the first step remains the same as in Ref. Baldsiefen et al. (2015) and therefore we refrain from repeating it. It has to be noted that it yields the equation:

[TABLE]

for every compact operator $\hat{A}$.

II. $\varrho$* yields the pair $(\gamma,\chi)$.*

For clarity’s sake the reduced density matrices corresponding to $\varrho_{k}$ are labeled $\gamma_{k}$ and $\chi_{k}$ and the ones corresponding to $\varrho$ are denoted as $\bar{\gamma}$ and $\bar{\chi}$ even though they are all equal. We will use the weak-* convergence of $\varrho_{k}$ to prove the weak-* convergence of $(\gamma_{k},\chi_{k})$ which yields the strong convergence. $\left\{\gamma_{k}\right\}$/$\left\{\chi_{k}\right\}$ is weakly-* convergent if and only if all linear functionals $\mathcal{L}\left[\gamma_{k}\right]$/$\mathcal{L}\left[\chi_{k}\right]$ on $\mathfrak{G}^{N}$ converge to $\mathcal{L}\left[\bar{\gamma}\right]$/ $\mathcal{L}\left[\bar{\chi}\right]$. In other words if for all bounded functions $f$:

[TABLE]

and

[TABLE]

In order to connect the anomalous densities with their statistical density operator we use the following operator:

[TABLE]

Now, we can write Eq. (32) as:

[TABLE]

Unfortunately, we cannot apply Eq. (30) as $\hat{O}_{f}$ is not a compact operator. Nonetheless, as the eigenvalues of $\hat{O}_{f}$ are finite, the product of $\hat{O}_{f}$ with a finite dimensional projection operator $\hat{P}$ is compact. Using a projection $\hat{P}$ we transform the left hand of Eq. (34) to:

[TABLE]

We can show separately that each part goes to zero. By applying Eq. (30) to the left part of Eq. (35) we can see that it goes to zero. For the right part we prove separately that the terms with $\varrho$ and $\varrho_{k}$ both go to zero. To start with, we can choose $\hat{P}$ such that:

[TABLE]

where M depends on the dimension of $\hat{P}$ and $w_{ki}$ are the eigenvalues of $\varrho_{k}$. This means we can choose a $K$ such that

[TABLE]

for all $k\geq K$ and $\epsilon\geq 0$ and we can choose an $M$ such that the same is true for

[TABLE]

The same process can be repeated for $\gamma$. Indeed, by choosing the larger $K$ and $M$ from the two proofs we fulfill the conditions for the weak convergence of the pair $(\gamma_{k},\chi_{k})$. In other words the weak-* convergence of $\varrho_{k}$ implies the weak convergence of $(\gamma_{k},\chi_{k})$ and because we choose all $(\gamma_{k},\chi_{k})=(\gamma,\chi)$ it implies strong convergence and $\varrho\rightarrow(\gamma,\chi)$.

III. * $\mathcal{F}_{\beta}[\gamma,\chi]={\Tr}\hat{\varrho}\left(\hat{H}+\frac{1}{\beta}\ln\hat{\varrho}\right)$ for this $\varrho$.*

The third and last step is once again completely equivalent to Baldsiefen (2012).

Therefore, there is in fact a minimizing $\varrho$ in Eq. (29) for all $\Gamma\in\mathfrak{G}^{N}$, and we can replace the infimum by a minimum, and the functional can be written as

[TABLE]

This extension allows us to define a Levy-Lieb functional for SC-RDMFT, namely:

[TABLE]

This functional justifies functional differentiability within SC-RDMFT, which is required both to develop a Kohn-Sham framework and to engineer exchange-correlation functionals. Since $\Gamma$ and $(\gamma,\chi)$ are uniquely determined we will not distinguish between $\mathfrak{F}_{\beta}[\gamma,\chi]$ and $\mathfrak{F}_{\beta}[\Gamma]$ (or from $\mathcal{F}_{\beta}[\gamma,\chi]$ and $\mathcal{F}_{\beta}[\Gamma]$), from now on. A second important property of the functional (38) is its convexity, which follows directly from the convexity of the grand potential $\Omega_{v,\Delta}[\gamma,\chi]$.

A crucial problem in reduced density theories (DFT or RDMFT) is the problem of $v$-representability. We already defined the set $\mathfrak{G}^{v,\Delta}\subset\mathfrak{G}^{N}$ in Eq. (26) as the set of $(v,\Delta)$-representable Nambu-Gorkov reduced density matrices. For the following we will use a different but equivalent definition: $(\gamma_{0},\chi_{0})\in\mathfrak{G}^{v,\Delta}$ if there exists $(v,\Delta)$ such that

[TABLE]

Remarkably, $\mathfrak{G}^{v,\Delta}$ is dense in $\mathfrak{G}^{N}$. This result is crucial when deriving functionals through many-body perturbation theory, as we will do in the following. For instance, in $\mathfrak{G}^{v,\Delta}$ we can relatively safely assume that a grand potential approximation using a functional derived through many-body perturbation theory, specifically the Sham-Schlüter connection Godby et al. (1986, 1987a, 1987b), is minimized by the correct equilibrium Nambu-Gorkov one-body density matrix. We already discussed that, in practice, a minimization has to be carried out in the set of ensemble-$N$-representable Nambu-Gorkov one-body reduced density matrices $\mathfrak{G}^{N}$.

For the proof of the density of $\mathfrak{G}^{v,\Delta}$ on $\mathfrak{G}^{N}$ we require the concept of a continuous tangent functional (CTF): Consider a real functional $\mathcal{F}$ on a subset $\mathcal{S}$ of a Banach space $\mathfrak{H}$. Let $x_{0}\in\mathcal{S}$. A continuous linear functional $\mathcal{L}$ is said to exhibit a CTF at $x_{0}$ if Engel and Dreizler (2011):

[TABLE]

The crucial observation is the following theorem:

Theorem 2**.**

The universal functional $\mathfrak{F}_{\beta}$ has a unique CTF at every pair $(\gamma,\chi)\in\mathfrak{G}^{v,\Delta}$ and only at these pairs.

The proof of this latter statement follows by reductio ad absurdum. Consider a pair $(\gamma^{\prime},\chi^{\prime})$ that is not $(v,\Delta)$-representable and suppose that $\mathfrak{F}_{\beta}\left[\gamma^{\prime},\chi^{\prime}\right]$ has a CTF denoted as $(\tilde{v},\tilde{\Delta})$. Therefore,

[TABLE]

Notice the change of sign in front of the anomalous term in the Hamiltonian (16). Thus,

[TABLE]

The infimum is only assumed for $(\chi,\gamma)\in\mathfrak{G}^{v,\Delta}$ and only then is the equality fulfilled. This contradicts our assumption that $(\gamma^{\prime},\chi^{\prime})$ is not equilibrium $(v,\Delta)$-representable. In this sense, the existence of a CTF at $(\gamma^{\prime},\chi^{\prime})$ implies the equilibrium $(v,\Delta)$-representability of $(\gamma^{\prime},\chi^{\prime})$. The opposite direction follows trivially from the definition of equilibrium $(v,\Delta)$-representability.

One last theorem, a less known variant of the famous Bishop-Phelps theorem Bishop and Phelps (1986), is required to conclude the proof. Let $\mathcal{F}$ be a lower semi-continuous convex functional on a real Banach space $\mathfrak{H}$. Suppose $x_{0}=(\gamma_{0},\chi_{0})\in\mathfrak{H}$ and $V_{0}=(v_{0},\Delta_{0})\in\mathfrak{H}^{*}$. Then, for every $\epsilon>0$, there exists $x_{\epsilon}\in\mathfrak{H}$ and $V_{\epsilon}\in\mathfrak{H}^{*}$ such that Bishop and Phelps (1986):

$||V_{\epsilon}-V_{0}||\leq\epsilon$. 2. 2.

$V_{\epsilon}$ is CTF to $\mathcal{F}$ at $x_{\epsilon}$. 3. 3.

$\epsilon||x_{\epsilon}-x_{0}||\leq\mathcal{F}\left[x_{0}\right]+v_{0}[\gamma_{0}]+\Delta_{0}[\chi_{0}]-\underset{x\in\mathfrak{H}}{\inf}\left(\mathcal{F}\left[x\right]+v_{0}[\gamma]+\Delta_{0}[\chi]\right)$.

This theorem can be used to prove the following result, which is equivalent to our hypothesis that $\mathfrak{G}^{v,\Delta}$ is dense in $\mathfrak{G}^{N}$:

Theorem 3**.**

For all $x\in\mathfrak{G}^{N}$ there exists a sequence $\left\{x_{k}\right\}\subset\mathfrak{G}^{N}$, such that $x_{k}\rightarrow x$ and $x_{k}\in\mathfrak{G}^{v,\Delta},\,\forall\,k$.

To prove this, we start by denoting the right side of the third condition in the Bishop-Phelps Theorem as $C_{0}$, which is independent of $\epsilon$ and finite. Let $\epsilon_{k}=kC_{0}$. We can find an $x_{k}$ for each $\epsilon_{k}$ at which $\mathcal{F}$ has a CTF such that

[TABLE]

We know from Theorem 2 that the existence of a CTF to $\mathcal{F}$ at $x_{k}$ is equivalent to the equilibrium $(v,\Delta)$-representability of $x_{k}$. This concludes not only the proof of the theorem but also our original task of proving that the set $\mathfrak{G}^{v,\Delta}$ is dense in the set $\mathfrak{G}^{N}$.

V Kohn-Sham system for SC-RDMFT

After deriving a Hohenberg-Kohn theorem, the next logical step is to discuss a Kohn-Sham system for SC-RDMFT. We develop two different approaches to accomplish such a goal. The first is analogous to SC-DFT (or normal DFT as well) and works through functional derivatives. However, just like in SC-DFT it does not prove the existence of the Kohn-Sham system. The second approach, on the other hand, is analogous to finite-temperature RDMFT and yields a proof for the existence of the Kohn-Sham system.

V.1 Approach 1: analogous to SC-DFT

The following discussion is only valid at finite temperature, since no Kohn-Sham system exists at zero temperature Schmidt et al. (2019). In SC-DFT a Kohn-Sham system is defined as a non-interacting system that reproduces the electronic density $n_{\rm eq}$ and the anomalous density $\chi_{\rm eq}$ of the interacting system. In the same vein, we define the Kohn-Sham system in SC-RDMFT as the non-interacting system with the same (equilibrium) reduced density matrices $(\gamma,\chi)$. For systems with spin-rotational symmetry, the Nambu-Gorkov one-body reduced density matrix

[TABLE]

is equivalent to such densities. As a result of the one-to-one relation (19) the grand potential can be written as:

[TABLE]

where

[TABLE]

We divide the universal functional in meaningful parts:

[TABLE]

where $U_{\text{Hartree}}$ is the Hartree term of the total electron-electron interaction energy and $S_{\text{KS}}$ is the entropy of the Kohn-Sham system. A natural advantage of RMDFT, unlike DFT, is of course our explicit knowledge of the kinetic energy functional.

We denote the grand-canonical potential of the Kohn-Sham system as:

[TABLE]

Our aim is now to determine the Kohn-Sham potentials $v_{\text{KS}}$ and $\Delta_{\text{KS}}$ such that the Kohn-Sham grand canonical potential $\Omega_{\text{KS}}\left[\gamma,\chi\right]$ is minimized by the same $(\gamma,\chi)$ of the interacting system. As in normal DFT, at this point one needs the assumption that at finite temperatures the two relevant interacting densities $\gamma$ and $\chi$ are smooth (i.e., normalizable and twice differentiable) and therefore noninteracting $(v,\Delta$)-representable Engel and Dreizler (2011). We then obtain:

[TABLE]

and

[TABLE]

and therefore

[TABLE]

with $v_{\text{xc}}(\bm{r},\bm{r}^{\prime})=\delta\mathcal{F}_{\text{xc}}/\delta\gamma(\bm{r},\bm{r}^{\prime})$. The Kohn-Sham potential consists of an external potential, a Hartree term, and the exchange-correlation term. After analogous calculations for $\chi(\bm{r},\bm{r}^{\prime})$ we arrive at:

[TABLE]

Notice that the universal exchange-correlation functional energy $\mathcal{F}_{\text{xc}}$ is comprised of first the exchange-correlation energy of the electrons and second the difference between the interacting and non-interacting entropy. The existence of the Kohn-Sham potentials hinges on the existence of the functional derivatives of $\mathcal{F}_{\rm xc}\left[\gamma,\chi\right]$. Although we cannot answer this question in general, in practice we will use differentiable approximations for $\mathcal{F}_{\text{xc}}$.

Diagonalizing the Kohn-Sham Hamiltonian by a Bogoliubov transformation of the field operators $\hat{\psi}_{\sigma}$:

[TABLE]

where $\hat{\gamma}^{\dagger}_{k,\sigma}$ and $\hat{\gamma}_{k,\sigma}$ are creation and annihilation operators of fermionic quasiparticles, yields the following Bogoliubov-de-Gennes equations which are completely analogous to SC-DFT and only differ in the non-locality of the Kohn-Sham potential:

[TABLE]

These are the self-consistent SC-RDMFT equations for the Kohn-Sham orbitals. In Appendix A we present the full Nambu-Gorkov Green’s functions for SC-RDMFT.

V.2 Approach 2: analogous to finite-temperature RDMFT

In finite temperature RDMFT the existence of a Kohn-Sham system is proven by introducing a common system of eigenfunctions of the Kohn-Sham Hamiltonian and the one-body reduced density matrix. Unfortunately, the Hamiltonian of the superconducting Kohn-Sham system contains two anomalous terms ($\Delta$ and $\Delta^{*}$) which do not commute with the particle-number operator $\hat{N}$. However, we can derive the Kohn-Sham system through a second more original approach we shortly describe in this subsection, by directly searching for the Bogoliubov transform that diagonalizes $\Gamma$ in the following form:

[TABLE]

Such a Bogoliubov transform exists for every ensemble $N$-representable $\Gamma$ Bach et al. (1994). In this new auxiliary Bogoliubov-de-Gennes system, represented by only $\gamma$ (as the anomalous terms are zero) we can basically return to the normal finite-temperature RDMFT and the corresponding proof for the existence of a Kohn-Sham system Baldsiefen et al. (2015). We define new field operators belonging to the creation and annihilation operators $\hat{\gamma}_{i}$ and $\hat{\gamma}^{\dagger}_{i}$ of the Bogoliubov quasiparticles:

[TABLE]

where $\tilde{\phi}$ are the natural orbitals (eigenfunctions) of $\Gamma_{\rm diag}$. The density matrix $\gamma(\bm{r},\bm{r}^{\prime})$ is defined as:

[TABLE]

These equations highlight that the occupation numbers of the Bogoliubov quasiparticles $\lambda_{i}$ belong to eigenenergies $\epsilon_{i}$ and are connected through a Fermi-Dirac distribution: $\lambda_{i}=1/(1+e^{\beta(\epsilon_{i}-\mu)})$. Thus we arrive at an explicit equation for a Kohn-Sham potential in the Bogoliubov-de-Gennes system:

[TABLE]

where $t_{ij}$ are the entries of the kinetic-enegy operator. Obviously, the Hamiltonian can be written as $\hat{H}=\sum_{i}\epsilon_{i}\hat{\gamma}_{i}^{\dagger}\hat{\gamma}_{i}$. We now recognize this form as the result of the Bogoliubov transform of the superconducting Kohn-Sham Hamiltonian. If we reverse the Bogoliubov transform that diagonalized $\Gamma$ we arrive at Kohn-Sham potentials that fulfill our familiar Bogoliubov-de-Gennes equations (52). Consequently, unlike SC-DFT, we can solve the problem of non-interacting $(v,\Delta)$-representability at finite temperature in SC-RDMFT as we can prove that a Kohn-Sham system exists for every $\Gamma\in\mathfrak{G}^{N}$ at finite temperature.

Furthermore, due to the one-to-one mapping between $(v,\Delta)$, $(\gamma,\chi)$, and $\varrho_{\rm eq}$ we know that both ways of reaching a Kohn-Sham system result in exactly the same Kohn-Sham potentials and equilibrium state. It is worth noticing that at zero temperature the situation is completely different. Indeed, since at zero temperature the ground state of a normal non-interacting system is a Slater determinant, the corresponding one-body reduced density matrix is idempotent. The same holds for a superconducting system: it turns out that the Nambu-Gorkov one-body reduced density matrix is also idempotent, which implies:

[TABLE]

Therefore, $\chi$ and $\gamma$ are not independent. However, in perturbation theory, for an interacting system $\Gamma$ is idempotent in first order and only ceases to be so in second order Schmidt et al. (2019). The idempotence of $\Gamma$ in first-order perturbation theory hints at the fact that even though a Kohn-Sham system does not exist at zero temperature, a practical implementation might nevertheless not fail completely.

VI Decoupling Approximation

In this section we develop a useful computational framework of SC-RDMFT. As a starting point, we solve the electronic problem with the Kohn-Sham potential $v_{\text{KS}}\left[\gamma,\chi\right]$. Assuming a crystal lattice we can write the natural orbitals as Bloch states. Furthermore, we denote the functions $u_{i}(\bm{r})$, $v_{i}(\bm{r})$, and $\Delta(\bm{r},\bm{r}^{\prime})$ in the basis of the natural orbitals of the non-superconducting solution, namely:

[TABLE]

In what follows, we will use this eigenbasis of natural Bloch orbitals and the solution of the normal electronic system to solve the phononic and superconducting problem.

VI.1 Electron-phonon coupling

Several approximations are needed in order to treat phonons in SC-RDMFT. First of all, it is reasonable to assume that the atoms only move from their equilibrium lattice positions through small oscillations and therefore the usual harmonic approximation applies. Second, the lattice dynamics will be approximated by the ones of the corresponding non-superconducting system. In this way, we can define the electron-phonon scattering matrix elements in a completely analogous way to SC-DFT Sanna (2017):

[TABLE]

The main difference lies in the fact that $\phi_{{n\bm{k}}}$ are not DFT-Kohn-Sham orbitals but natural orbitals. In addition, the Kohn-Sham potential is non-local. In Eq. (62), $\bm{q}$ and $\bm{k}$ are the phonon and electron momenta, $n$ and $m$ are the natural-orbital band indices, $\omega_{\nu}$ the phonon frequency, $\nu$ the phonon branch and $\Delta V_{\rm scf}^{\bm{q}\nu}$ the variation in the Kohn-Sham potential due to ionic displacement. This results in the Hamiltonian:

[TABLE]

where $\hat{b}_{\nu\bm{q}}$ and $\hat{b}^{\dagger}_{\nu\bm{q}}$ are the annihilation and creation operators of the phonons. The form of the electron-phonon coupling in SC-RDMFT and SC-DFT seems very similar. However, we should notice that the electron-phonon coupling constants that enter our theory should, in principle, be calculated from RDMFT. While in DFT the matrix elements $g^{\nu}_{m\bm{k}+\bm{q},{n\bm{k}}}$ are easily obtained from density functional perturbation theory Giustino (2017), such a framework still does not exist for RDMFT for solids. Yet, as discussed in SC-DFT Lüders et al. (2005), to solve the corresponding gap equation are required the electron-phonon coupling constants, $g^{\nu}_{m\bm{k}+\bm{q},{n\bm{k}}}$, as well as the normal-state Kohn-Sham eigenenergies $\xi_{{n\bm{k}}}$. In practice, however, the phononic contributions to the functionals are often averaged on the Fermi surface Marques et al. (2005):

[TABLE]

VI.2 Band-decoupling approximation

In SC-DFT the decoupling approximation is based on the assumption that the superconducting transition does not introduce any structural transition. It is also assumed that it does not cause any hybridization between the bands. We already assumed that the first approximation is true when we treated the phonons. However, the first approximation also includes the fact that $\Delta_{\text{KS}}(\bm{r},\bm{r}^{\prime})$ has the periodicity of the lattice and it is therefore unnecessary to sum over the indices $k$. The second assumption reduces the equations for $u_{i}(\bm{r})$ and $v_{i}(\bm{r})$ to:

[TABLE]

This approximation is motivated by the difference in energy scale between the electronic bonding and superconducting pairing. If any of the bands are degenerate on the scale of $\Delta_{\text{KS}}$ one can question the validity of the approximation. Yet, even then, we expect the hybridization to be negligible.

Inserting these approximations into the Bogoliubov-de-Gennes equations (52), and using the orthogonality of the basis set we arrive at:

[TABLE]

The energies $\tilde{\epsilon}_{{n\bm{k}}}$ are defined as the energies from the non-superconducting RDMFT calculation minus the chemical potential. This yields equations of the same form as in SC-DFT. Indeed, for the eigenenergies we obtain:

[TABLE]

and for the amplitudes:

[TABLE]

where $e^{\Phi_{{n\bm{k}}}}$ is the phase of $\Delta_{{\rm{KS}};{n\bm{k}}}$. We are already familiar with the first equation from BCS-theory where $\Delta_{{\rm{KS}},{n\bm{k}}}$ is the energy gap. We will interpret $\Delta_{{\rm{KS}},{n\bm{k}}}$ in the same fashion in SC-RDMFT. Of course the solution of the Bogoliubov-de-Gennes equations also results in new formulas for $\gamma$ and $\chi$. Here we use the fact that the Bogoliubov-quasiparticle states are occupied according to a Fermi-Dirac distribution. The densities $\gamma$ and $\chi$ take on the form:

[TABLE]

and

[TABLE]

We defined the anomalous exchange-correlation potential $\Delta_{\text{xc}}$ as the functional derivative of the universal functional with respect to $\chi$. This, in combination with Eq. (69), yields a self consistent equation for the anomalous exchange-correlation potential:

[TABLE]

Here we used the fact that $\chi$ is a functional of $\Delta$ and $\gamma$ and consequently we can calculate the functional derivative in Eq. (71) through the chain rule as a derivative with respect to $\Delta_{\text{KS}}$ and $\gamma$. As we consider $\Delta(\bm{r},\bm{r}^{\prime})$ as the superconducting energy gap we arrive at a self-consistent gap equation which we can solve in order to determine the transition temperature.

VII The Sham-Schlüter connection

The last missing piece in the SC-RDMFT puzzle is an expression for the universal exchange-correlation functional. The Sham-Schlüter connection Godby et al. (1986, 1987a, 1987b) is a many-body perturbation approach for the derivation of exchange-correlation potentials. The idea was already introduced to SC-DFT in the original paper by Oliveira, Gross and Kohn Oliveira et al. (1988). We will use the existence of the Kohn-Sham system to show that the Sham-Schlüter connection can be easily extended to SC-RDMFT.

The starting point is the Dyson equation in Nambu-Gorkov space:

[TABLE]

where $\bar{G}^{\text{KS}}$ is the non-interacting Green’s function which corresponds to the SC-RDMFT Kohn-Sham Hamiltonian:

[TABLE]

$\bar{\Sigma}^{\text{KS}}(\bm{r},\bm{r}^{\prime},\omega_{i})$ is the self energy of the system:

[TABLE]

with a double counting correction $\bar{\Sigma}^{DC}(\bm{r},\bm{r}^{\prime},\omega_{i})$ that is comprised of the SC-RDMFT exchange-correlation potentials:

[TABLE]

We now arrive at the reason for using the Kohn-Sham system as the non-interacting system of the Dyson equation. Indeed, the densities $(\gamma,\chi)$ can be written as the equal time limit of the Nambu-Gorkov Green’s function. Here, $\gamma$ corresponds to the normal Green’s function and $\chi$ to the anomalous propagator:

[TABLE]

Obviously, this is true for the Nambu-Gorkov Kohn-Sham Green’s function as well as for the Nambu-Gorkov Green’s function of the interacting system. Since we defined the Kohn-Sham system as the non-interacting system reproducing $\chi$ and $\gamma$ of the interacting system, the equal time limit of $\bar{G}^{\text{KS}}(\bm{r}\tau,\bm{r}^{\prime}\tau^{+})=\bar{G}(\bm{r}\tau,\bm{r}^{\prime}\tau^{+})$. Inserting the equality into the corresponding Dyson equation yields the following two equations:

[TABLE]

The analogue of this relation for the density is known for DFT Godby et al. (1986, 1987a, 1987b) as the Sham-Schlüter connection. After decomposing the self-energy we arrive at the following system of integral equations for the exchange-correlation potentials. As a matter of notation we use $\bar{G}^{\text{KS}}(\bm{r},\bm{r}_{1},\omega_{i})=\bar{G}^{\text{KS}}_{\downarrow\downarrow}(\bm{r},\bm{r}_{1},\omega_{i})=\bar{G}_{\uparrow\uparrow}^{\text{KS}}(\bm{r},\bm{r}_{1},\omega_{i})$ and analogously for $F(\bm{r},\bm{r}^{\prime},\omega_{i})$ and $\Sigma^{\text{xc}}(\bm{r},\bm{r}^{\prime},\omega_{i})$.

[TABLE]

and

[TABLE]

One can slightly simplify these equations by transforming to Fourier-space (the convergence factor $\lim_{\eta\to 0^{+}}$ is left out for simplicity’s sake). Indeed, by defining:

[TABLE]

the decoupling approximation reduces the anomalous potential to only one index ${n\bm{k}}$, namely:

[TABLE]

At this point the first major advantage in comparison to SC-DFT emerges. In SC-DFT the equation corresponding to the density is of course local while the anomalous equation is non-local. The two equations arising from the Sham-Schlüter connection in SC-RDMFT are symmetric in the sense that they are both purely non-local. Consequently, we do not sum over any index ${n\bm{k}}$ in Eq. (85) and we are not plagued with the SC-DFT-problem of reconciling a local and a non-local equation. This greatly simplifies the calculations: the off-diagonal elements of $v_{\text{xc}}({n\bm{k}},{n\bm{k}}^{\prime})$ do not depend on the anomalous potential and can be calculated separately. Moreover, if $v_{xc}$ satisfies the lattice periodicity (namely, $v_{xc}(\bm{r},\bm{r}^{\prime})=v_{xc}(\bm{r}+\bm{R},\bm{r}^{\prime}+\bm{R})$ with $\bm{R}$ being any lattice vector), the nondiagonal elements are non-zero only when they belong to different bands $v_{\text{xc}}({n\bm{k}},\bm{n}^{\prime}\bm{k})$. Since the off-diagonal component $(\bm{k},\bm{k}^{\prime})$ becomes zero under the lattice periodicity of the exchange-correlation potential, Eqs. (84) and (85) can be further simplified.

In order to continue analogously to SC-DFT we simplify the equation by limiting ourselves to terms linear in $\Delta_{{n\bm{k}}}$:

[TABLE]

Naturally, this approximation is only valid close to the transition temperature where the energy gap is small. As a last approximation, we replace the interacting Green’s function with the Kohn-Sham Green’s function. From the perspective of Migdal’s theorem Migdal (1958), this approximation is not completely sound. For this reason, a promising avenue for further research will be to dress the Green’s function with the phononic self energy Marques et al. (2005). The final form of the functional is derived in detail in Appendix B.

We finally arrive at the linearized gap equation

[TABLE]

that we divide into one Coulomb term,

[TABLE]

and two phonon terms

[TABLE]

and

[TABLE]

The Coulomb term (88) and phonon terms (89) and (90) are completely analogous to the corresponding terms in SC-DFT. The only difference in the $A$ term lies in the fact that the terms $\xi_{{n\bm{k}}}$, $\Delta_{{n\bm{k}}}$ and $g^{{n\bm{k}}{n\bm{k}}^{\prime}}_{\lambda q}$ correspond to the results of RDMFT and not DFT calculations. In SCD-DFT there is an additional repulsive $B$ term, namely,

[TABLE]

which turns out to be rather problematic as it suddenly jumps to zero close to the Fermi surface and is consequently neglected Lüders et al. (2005); Marques et al. (2005). Therefore, it speaks rather favorably of SC-RDMFT that the term is zero without any ad hoc approximations. Once again the reason why the term is simpler in SC-RDMFT lies in the symmetry of the equations. The sums over ${n\bm{k}}_{1}$ in the last term arise due the use of the density in SC-DFT and consequently do not appear in SC-RDMFT.

The C-term (89) is exactly analogous to superconducting DFT term. The D-term (90) is different from SC-DFT. The difference lies in the first summand which was neglected in SC-DFT as it diverges Lüders et al. (2005). Although the second summand also has a divergent contribution in SC-RDMFT, the divergences of the first and second summands asymptotically cancel Akashi and Arita (2013). Consequently, the final form of the D-term is equal in both theories. Comparing this functional with its analogue from SC-DFT we only expect differences due to the different Kohn-Sham energies and electron-phonon coupling in SC-RDMFT. All in all, the process of developing the first anomalous exchange-correlation functional for SC-RDMFT turned out to be simpler and relied on fewer approximations to reach an equivalent result to SC-DFT.

VIII Conclusions

In this paper we presented the theoretical foundations of a new ab initio theory of superconductivity. SC-RDMFT is based on a unique relationship between the statistical density operator at equilibrium $\varrho_{\rm eq}$ and the corresponding one-body reduced density matrix $\gamma$ and the anomalous density $\chi$. This reduced density matrix formalism for superconductivity yields the existence of a universal functional $\mathfrak{F}_{\beta}[\gamma,\chi]$ of these two quantities, whose universal properties we derived. The derived theorems prove the possibility of a completely ab initio reduced density matrix formalism for superconductors and demonstrate the elegance and advantages of using the Nambu-Gorkov one-particle reduced density matrix or its substitutes to describe superconducting systems.

We used a Bogoliubov transform to prove the existence of the Kohn-Sham system at finite temperature and obtain a system of coupled Bogoliubov-de Gennes-like Kohn-Sham equations. As even such coupled Kohn-Sham equations in their original form are extremely challenging to solve we introduced the decoupling approximation to SC-RDMFT. By decoupling the changes of the superconducting phase transition from the electronic and phononic calculations we were able to arrive at a linear system of equations. Solving the system leads to a set of self-consistent equations for the gap and our primary variables $\gamma$ and $\chi$.

We then used the existence of the Kohn-Sham system to derive an exchange-correlation functional based on the Sham-Schlüter connection in Nambu-Gorkov space. Formally the derivation was extremely similar to SC-DFT but a few essential differences existed. Because SC-RDMFT is concerned with two non-local variables, the resulting equations exhibited a symmetry in their form that is missing in SC-DFT. Due to this missing symmetry an extra strong approximation is necessary in SC-DFT, namely that the system is nearly homogeneous. In the end we simplified our functional to the linear regime in order to arrive at a BCS-like gap equation. The resulting approximate functional has three terms that have the same analytical form as the analogous functional in SC-DFT. However, in SC-DFT there is a fourth term that has a rather pathological behavior, that is set to zero by hand. This problematic term turns out to be absent in SC-RDMFT.

Our functional was developed in close analogy to SC-DFT. It would be certainly interesting to follow a different path, and to derive it by generalizing existing (non-superconducting) RDMFT functionals (such as the Müller functional). In view of the success of RDMFT for correlated systems, this would maybe allow us to tackle, for the first time, the problem of superconductivity in correlated systems in a completely ab initio fashion.

Acknowledgements.

We thank Hardy Gross for helpful discussions.

Appendix A Nambu-Gorkov Green’s function

Just like the Nambu-Gorkov one-body reduced density matrix, the Green’s function in Nambu-Gorkov space is a $2\times 2$ matrix:

[TABLE]

Here $G_{\sigma,\sigma^{\prime}}(\bm{r}\tau,\bm{r}^{\prime}\tau^{\prime})$ corresponds to the normal Green’s function and $F^{\dagger}_{-\sigma,\sigma^{\prime}}(\bm{r}\tau,\bm{r}^{\prime}\tau^{\prime})$ is known as the anomalous propagator.

We can denote the constituents of the Kohn-Sham Green’s function as:

[TABLE]

in position space. A Fourier transformation assuming the decoupling approximation yields:

[TABLE]

Appendix B Linear Regime and Frequency Sums

In the linear regime we arrive at:

[TABLE]

Fortunately, Marques Marques (2000) already calculated all the required frequency sums analytically 111$\tilde{\epsilon}_{{n\bm{k}}}$ is replaced with $\xi_{{n\bm{k}}}$ to avoid any confusion with the other variables in the linear regime:

[TABLE]

[TABLE]

[TABLE]

For the electronic self-energy expressions we arrive at:

[TABLE]

[TABLE]

[TABLE]

And the phononic terms yield:

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

and $D_{\lambda q}(\omega_{1}-\omega_{2})$ being the phonon propagator.