Calculating Bardeen-Cooper-Schrieffer and Magnetic Superstructure Electronic States with $\Theta\Phi$

TL;DR

The paper introduces the $ heta ext{ extPhi}$ package that enables first-principles calculations of superconducting and magnetic states with complex orderings in multi-band systems, bridging mean-field theory and quantum solid state physics.

Contribution

It presents a novel computational tool that incorporates BCS superconductivity and magnetic superstructures into first-principles calculations for multi-band systems.

Findings

Demonstrates the implementation with example calculations.

Allows for first-principles estimates of model Hamiltonian matrix elements.

Facilitates study of interplay between superconductivity and magnetism.

Abstract

We propose the $\Theta\Phi$ (Theta-Phi) package which addresses two of the most important extensions of the essentially single-particle mean-field paradigm of the computational solid state physics: the admission of the Bardeen-Cooper-Schrieffer electronic ground state and allowance of the magnetically ordered states with an arbitrary superstructure (pitch) wave vector. Both features are implemented in the context of multi-band systems which paves the way to an interplay with the solid state quantum physics packages eventually providing access to the first-principles estimates of the relevant matrix elements of the model Hamiltonians derived from the standard DFT calculations. Several examples showing the workability of the proposed code are given.

Tables

Table 1. Table 1: k 𝑘 k -dependence of the gap function in case of extended- s 𝑠 s , d x 2 − y 2 subscript 𝑑 superscript 𝑥 2 superscript 𝑦 2 d_{x^{2}-y^{2}} , and d x ​ y subscript 𝑑 𝑥 𝑦 d_{xy} -wave superconductivity on 2D honeycomb lattice. The vectors { δ i } subscript 𝛿 𝑖 \{\delta_{i}\} , connecting a graphene site to its three nearest neighbors are defined in the text.

	extended-$s$	$d_{x^2-y^2}$	$d_{xy}$
$\Delta \left(𝐤\right)\sim$	$e^{i𝐤{\delta }_1}+e^{i𝐤{\delta }_2}+e^{i𝐤{\delta }_3}$	$e^{i𝐤{\delta }_2}-e^{i𝐤{\delta }_3}$	$2e^{i𝐤{\delta }_1}-e^{i𝐤{\delta }_2}-e^{i𝐤{\delta }_3}$

Table 2. Table 2: Total energy gain E t ​ o ​ t subscript 𝐸 𝑡 𝑜 𝑡 E_{tot} and Fe magnetic moment comparison between Θ ​ Φ Θ Φ \Theta\Phi and VASP for ferromagnetic calculations in α 𝛼 \alpha -Fe.

	NM	FM ($\Theta \Phi$)	FM(VASP)
$\Delta E_{tot}$ (eV)	$0$	$-6.4$	$-5.0$
$\mu$ (${\mu }_B$)	$0$	$3.1$	$3.2$

Table 3. Table 3: Total energy gain Δ ​ E t ​ o ​ t Δ subscript 𝐸 𝑡 𝑜 𝑡 \Delta E_{tot} and Fe d 𝑑 d -shell magnetic moment comparison between Θ ​ Φ Θ Φ \Theta\Phi and VASP for ferromagnetic and anti-ferromagnetic [ 𝐐 = π ​ ( 0 , 1 , 1 ) ] delimited-[] 𝐐 𝜋 0 1 1 \left[\mathbf{Q}=\pi(0,1,1)\right] calculations in γ 𝛾 \gamma -Fe.

	NM	FM ($\Theta \Phi$)	AFM ($\Theta \Phi$)	FM(VASP)	AFM (VASP)

Equations

$$H=H_K+H_U+H_J+H_V,$$

$$H_K=\sum _{\begin{array}{c}l,l^{\prime },R\\ \tau ,s,s^{\prime }\end{array}}{t}_{l,l^{\prime }}^{s,s^{\prime }}(\tau )c_{l,R,s}^{†}{c}_{l^{\prime },R+\tau ,s^{\prime }}^{},$$

$$H_U=\sum _{\begin{array}{c}l,l^{\prime },R\\ s,s^{\prime }\end{array}}{U}_{l,l^{\prime }}^{s,s^{\prime }}{n}_{l,R,s}{n}_{l^{\prime },R,s^{\prime }},$$

$$H_J=\sum _{\begin{array}{c}l,l^{\prime },R\\ \tau ,\alpha ,\beta \end{array}}{J}_{l,l^{\prime }}^{\alpha \beta }(\tau )S_{l,R}^{\alpha }{S}_{l^{\prime },R+\tau }^{\beta }.$$

$$H_V=\sum _{l,l^{\prime },R,\tau }{V}_{l,l^{\prime }}(\tau )n_{l,R}{n}_{l^{\prime },R+\tau }.$$

$$n_{l,R}$$

$$S_{l,R}^{\alpha }$$

$${\Psi }_R^{†}=\left(c_{1R\downarrow }^{†},c_{1R\uparrow }^{†},...,c_{LR\downarrow }^{†},c_{LR\uparrow }^{†},c_{1R\downarrow }^{},c_{1R\uparrow }^{},...,c_{LR\downarrow }^{},c_{LR\uparrow }^{}\right)$$

$$|{\Psi }_g⟩=\prod _{\begin{array}{c}j,k,\\ j<M\end{array}}{\gamma }_{j,k}^{†}|\tilde{0}⟩.$$

$${\Gamma }_k^{†}=\left({\gamma }_{1k}^{†},{\gamma }_{2k}^{†},\dots ,{\gamma }_{4Lk}^{†}\right).$$

$${\Gamma }_k={\Xi }_k{\Psi }_k.$$

$$|\tilde{0}⟩=\prod _{j,k}{c}_{j,k,\downarrow }^{†}|0⟩.$$

$$H_{var}=\frac{⟨{\Psi }_g|H-\mu N|{\Psi }_g⟩}{⟨{\Psi }_g|{\Psi }_g⟩}\to min.$$

$$A=H_{var}-TS\to \min .$$

$$N=\sum _{j,R}{n}_{j,R}.$$

$$\rho =\Phi \left[\rho \right],$$

$$H^{MF}$$

$$H^{MF}=\sum _{R,\tau }{\Psi }^{†}(R)\mathcal{H}(\tau )\Psi (R+\tau )+2L,$$

$$H^{MF}=\sum _k{\Psi }_k^{†}{\mathcal{H}}_k^{}{\Psi }_k^{},$$

$${\mathcal{H}}_k=\sum _{\tau }{e}^{ik\tau }\mathcal{H}(\tau ).$$

$$\rho (\tau )=1{\delta }_{\tau ,0}-⟨{\Psi }_R^{}{\Psi }_{R+\tau }^{†}⟩.$$

$${\rho }_k={\left(1+e^{\beta {\mathcal{H}}_k}\right)}^{-1},$$

$${\rho }_k=\sum _{\tau }{e}^{ik\tau }\rho (\tau ).$$

$${\rho }_{new}=\Phi [{\rho }_{old}],$$

$$\Vert {\rho }_{new}-{\rho }_{old}\Vert$$

$$\left|N-N_0\right|$$

$${\Omega }_{\mathbf{n},\vartheta }={\sigma }_0\cos \frac{\vartheta }{2}-i\left(\mathbf{n},\sigma \right)\sin \frac{\vartheta }{2},$$

$$\vartheta =\left(\tau ,\mathbf{Q}\right)$$

$$\Omega (\mathbf{n},\tau ,\mathbf{Q})$$

$${\Omega }^{†}(l)\Omega (\mathbf{n},\tau ,\mathbf{Q})\Omega (l^{\prime })$$

$$U=⟨H^{MF}⟩=\frac{{\sum }_k\mathrm{Tr}\left({\mathcal{H}}_k{\rho }_k\right)}{{\sum }_k\mathrm{Tr}{\rho }_k}.$$

Full text

††journal: Computer Physics Communications

1 INTRODUCTION

The solid state quantum physics packages available to the students in the field are all based on the Hartree-Fock approximation for the electronic wave function [1]. Extensions to it are restricted to the so called “post-Hartree-Fock” methods and largely reduce to the perturbation (Möller-Plesset order n

• MPn) corrections to the Hartree-Fock approximate ground state. This significantly restricts the repertory of the types of the ground states accessible to the available software. Practically, some of the important types of the electronic states of solids cannot be reproduced since they simply have not been programmed in. The most striking (and scandalous) examples of unaccessible states are the Bardeen-Cooper-Schrieffer (BCS) state [2, 3] necessary for description of superconductors. Even the intuitively more transparent states of solids - the magnetically ordered ones with an arbitrary pitch (superstructure) vector are not directly accessible by the available numerical tools. The magnetically ordered states can be obtained by extending the chemical unit cells to (magnetic) super-cells and setting the primeval magnetic moments with broken symmetry (BS) in the input file. Within this technology only simplest ordered magnetic states - ferromagnetic or antiferromagnetic with very simple magnetic super-lattices can be accessed. Thus, the presence of either BCS or complex magnetic states of solids is incurred indirectly. Namely, the presence of the occupied non-bonding one-electron states in the vicinity of the Fermi level indicates a possibility of that or another kind of instabilities of the primary symmetric structure leading to some BS solution not allowing to establish its specificity. Clearly, the BS solutions of non-programmed types never come to surface. Another feature so far missing in the available software is the temperature dependence of the solutions of the electronic problems. This feature is, however, important due to characteristic physical effects: transitions among the high-temperature symmetric and various low-temperature BS phases occurring in the experiment. Thus, we undertake the present development with a goal to heel the outlined deficiencies in the existing software. The paper is organized as follows: in Section 2 we present necessary theoretical concepts; in Section 3 we describe implementation details; in Section 4 we give test examples of applying developed software. In Section 5 we discuss the results of test calculations and give some perspectives.

2 THEORY ACCOUNT

2.1 Hamiltonian

The most general form of the Hamiltonian considered in $\Theta\Phi$ package is:

[TABLE]

where $H_{K}$ is the one-electron part, containing the kinetic energy and local terms:

[TABLE]

$H_{U}$ is the on-site multi-orbital Coulomb repulsion:

[TABLE]

and $H_{J}$ is the multi-site, multi-orbital Heisenberg term, which, in general, can be anisotropic:

[TABLE]

Finally, $H_{V}$ is the inter-site, inter-orbital Coulomb repulsion:

[TABLE]

Here, $n_{l,R,s}=c_{l,R,s}^{\dagger}c_{l,R,s}^{\phantom{\dagger}}$ is the electron density operator on $l$-th orbital of site $R$, while $n_{l,R}=n_{l,R,\uparrow}+n_{l,R,\downarrow}$ is the total density operator and $S_{l,R}^{\alpha}$ is the operator for the spin component $\alpha(=x,y,z)$:

[TABLE]

expressed as usual through Pauli matrices $\sigma_{s,s^{\prime}}^{\alpha}$.

In order to appropriately treat magnetic superstructure with arbitrary orbital-dependent pitch vectors, we implement individual spin quantization axis rotation for each orbital. Thus, there is no need to increase the size of unit cell at the expense of modifying the Hamiltonian terms. The orbital spin quantization axis rotation parameters together with the superstructure pitch vector become additional variational parameters and determine the optimal spin configurations. The details of this technique are given in subsection2.5.

2.2 Variational wave-function

The Hamiltonian (1) is a true many-body operator, impossible to solve exactly. We obtain an approximate variational solution by optimizing its expectation value with respect to the variational parameters of a trial wave-function within the mean-field solution.

The mean-field solution of eq.(1) is most easily obtained with the use of the Nambu formalism: that is replacing the Fermi operators $c_{l,R,s}^{\dagger}(c_{l,R,s}^{\phantom{\dagger}})$ creating an electron (hole) with the spin projection $\sigma$ in the state single state pertaining to the site (unit cell of the single band model) by a Nambu vector

[TABLE]

composed of the operators creating an electron (hole) with the spin projection $s=\downarrow,\uparrow$ in one of the states $l=1,\dots,L$ in the unit cell $R$. Obviously, since $\Psi_{R}$ is a column, then $\Psi_{R}^{\dagger}$ is a row consisting of corresponding hermitian conjugated operators. Here we include in the basis both spin projections as well as particle and hole creation operators in order to allow for superconducting terms in the Hamiltonian as well as off-diagonal spin exchange.

At this point we specify the class of the variational wave function used in $\Theta\Phi$ which is a generalization of Anderson’s RVB wave function:

[TABLE]

Here $\gamma_{j,k}$ are the so-called canonical quasi-particles:

[TABLE]

The number of operators in the Nambu vector $\Gamma_{k}$ is equal to the length of $\Psi_{k}$: $4L$. There exists a unitary transformation which relates the two:

[TABLE]

In eq.(5) we have introduced a new vacuum $\left|\tilde{0}\right\rangle$ which is defined in a way, similar to that of the canonical quasiparticles’ vacuum in BCS theory[4]:

[TABLE]

It is the coefficients of the matrices $\Xi_{k}$ which play the role of the variational wave function parameters, although they are not independent. In the present manuscript, we denote as $\left\{\zeta\right\}$ the set of all independent variables determining the $4L\times 4L$ elements of matrices $\Xi_{k}$ in eq. (6) at all $k$. The term “canonical” for quasiparticles in the context of $\left|\Psi_{g}\right\rangle$ means that this wave function is build up by filling the vacuum $\left|\tilde{0}\right\rangle$ with canonical quasi-particles. In other words, eq.(5) is a Slater determinant of canonical quasiparticles. In eq.(5) the filling occurs up to some band index $M$, which, in principle, has to be considered as another variational parameter. In the case of a particle-hole symmetric Hamiltonian, the optimal value of $M$ can be shown[5] to be $M=2L$ - a choice currently realized in $\Theta\Phi$. The form of the wave function (5) as a set of independent quasiparticles uniquely defines the thermodynamics of such an state: thermal excitations are quasiparticle-quasihole pairs. This corresponds to the partition function being the product of Fermi functions of quasiparticles energies determined in the self-consistent procedure described below.

2.3 Mean-field and self-consistency

The mean-field solution of eq.(1) is accomplished by minimizing the expectation value of either the variational energy (at zero-temperature):

[TABLE]

or of the Helmholtz free energy (at finite temperature) with respect to a set of variational parameters $\left\{\zeta\right\}$ (defined below).

[TABLE]

Here, $\mu$ is the chemical potential (taking care about the correct number of electrons - see below), $T$ is the temperature in the energy units (meaning $k_{B}=1$ in our notations), while $S$ is the (information/Shannon) entropy (defined below). The number of particles operator is defined as:

[TABLE]

The necessary (but not sufficient) condition for the minimum of the(Helmholtz) energy is the stationarity of the former with respect to the variational parameters $\left\{\zeta\right\}$ characterizing the trial wave function eq.(5). The stationarity conditions can be always expressed as a set of non-linear equations in the space of variational parameters. This can be symbolically written as an equation:

[TABLE]

where $\rho$ is the generalized one-particle density matrix of the system corresponding to the wave function eq.(5) and $\varPhi$ symbolically represents the mean-field self-consistency procedure allowing to derive the new density matrix given the old one. The number of particles constraint narrows down even more the search domain. In this restricted domain, the point giving the minimal (Helmholtz) energy represents the solution to the problem eq.(7). In addition, the equation $\left\langle N\right\rangle=N_{0}$ determines the chemical potential $\mu$ and should be considered as a complement to the system eq.(9).

Within the mean-field approach, the bare Hamiltonian $H$ is replaced with a one-electron Hamiltonian $H^{MF}$:

[TABLE]

The terms with the superscript $MF$ are obtained from the respective many-body operators by a “linearization” procedure (applying the Wick’s theorem[6]). The detailed forms of the linearized operators are given in the A. The fulfillment of the condition eq.(8) amounts to solve a non-linear optimization problem, which is a well-known challenge in both quantum chemistry and solid state physics. It is outline in the following subsections.

In terms of Nambu vectors the mean-field Hamiltonian reads:

[TABLE]

where the last term accounts for the anti-commutation relations of the operators comprising $\Psi_{R}$ ($L$ is the number of orbitals of the problem). It brings a constant energy shift (arising from fermionic anti-commutation rules) and will be omitted in what follows. The details of the contributions to the Hamiltonian used are also given in the A. Here $\mathscr{H}(\tau)$ is the matrix representation of $H^{MF}$ in the basis of operators, composing the Nambu vector $\Psi(R)$.

Finally, we take advantage of periodicity of the crystal and transform $H^{MF}$ to the reciprocal space:

[TABLE]

with

[TABLE]

2.4 Density matrix and fixing $\mu$

The central object in any mean-field theory is the density matrix, which is defined through the Nambu notations eq.(4) as:

[TABLE]

Here $\rho(\tau)$ are $4L\times 4L$ matrices, independent on $R$ due to translational invariance and $\mathbf{1}$ is a $4L\times 4L$ identity matrix. By using the hermiticity of$\mathscr{H}_{k}$ at each $k$ we can define an arbitrary analytic function of matrix argument in a usual way. As a consequence of the condition eq.(8), the density matrix can be compactly written as:

[TABLE]

where $\beta=1/T$ - inverse temperature. eq.(11) is a generalization of Fermi-Dirac distribution. Here, we have also introduced the Fourier transform of the density matrix:

[TABLE]

As described in Section 3, the density matrix contains all the parameters necessary to define the trial wave function eq.(5) and the electronic phases of the system. The optimization procedure is then equivalent to finding a self-consistent solution for the density matrix obtained from eq.(9)

[TABLE]

such that the newly generated $\rho_{new}$ does not deviate from $\rho_{old}$ at the previous iteration within an accuracy $\eta_{\rho}$, while the chemical potential is kept so that the total number of particles is within an accuracy $\eta_{\mu}$ equal to $N_{0}$:

[TABLE]

2.5 Rotations of the local quantization axes

In order to deal with arbitrary direction of spin-$\frac{1}{2}$ quantization axes for each local orbital we implement the standard quantization axis rotation formulae:

[TABLE]

which describes a rotation by angle $\vartheta$ around a rotation axis with the unit vector $\mathbf{n}=(n_{x},n_{y},n_{z})$. Here $\mathbf{\mathbf{\bm{\mathbf{\sigma}}}}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is the vector composed of Pauli matrices, while $\sigma_{0}$ is the $2\times 2$ identity matrix.

In paper [7] it is shown that the Bloch states for the electrons in a lattice can be built from the local states defined with respect to the local (unit cell related) axes of the spin quantization. These authors profit from the fact (explicitly not formulated) that rotations in the spin space around the $y$-axis form a commutative (sub)group (in fact the $SO(2)$ group) in the general group $SU(2)$ of the spin rotations and thus its elements can serve as images for a representation of the elements of another commutative group: that of the lattice translations. The same is, however, true for an arbitrary vector $\mathbf{n}$ fixed for all unit cells of the lattice not necessarily for $\mathbf{n}=(0,1,0)$. We use this option and allow the direction of $\mathbf{n}$ to be an optimized variable. Anyhow the relative rotation angle between two unit cells shifted by a lattice vector (one with integer components) $\mathbf{\mathbf{\tau}}$ is given by

[TABLE]

where $\mathbf{Q}$ pitch or superstructure (wave)vector (belongs to the reciprocal space). This defines the matrix

[TABLE]

which must appear inside each product of the vector creation/destruction operators representing the electron hopping between the unit cells separated by the shift vector $\mathbf{\tau}$. For $\mathbf{\tau}=0$, $\Omega(\mathbf{n,\tau,Q})=\sigma_{0}$.

It is not necessary, however, to assume that the spins within a unit cell are all quantized along the same axis. That is to say that one can assign to each orbital in the unit cell its own rotation matrix determining the direction of its quantization axis in the global “laboratory” frame. Let $SU(2)$-matrices $\Omega(l)$ and $\Omega(l^{\prime})$ be those which rotate the spin quantization axes of the $l$-th and $l^{\prime}$-th orbitals in the unit cell. Then for the pair of such orbitals when the $l^{\prime}$-th orbital is located in the unit cell shifted by the lattice vector $\mathbf{\tau}$ the matrix multiplier

[TABLE]

must be inserted between the fermion-vector multipliers representing the electron being destroyed in the $l^{\prime}$-th orbital and one being created in the $l$-th orbital. This matrix is as well of the form given by eq.(13). The complete formulae describing the rotation of various Hamiltonian terms are reported in the B.

2.6 Physical properties

Once the self-consistency is reached, the observables: correlation functions as well as band structure and electronic density of states (DOS) can be calculated. Below we summarize the formulas for some of them, but many others can be defined in addition. At the self-consistency, the internal energy $U$ is just an average of the mean field Hamiltonian:

[TABLE]

The (Helmholtz) free energy $A$ and the entropy $\varsigma$ are defined through the spectrum of the canonical quasi-particles:

[TABLE]

[TABLE]

The expression eq.(17) is used in practical calculations as anticipated above. The specific heat can also be easily obtained:

[TABLE]

The uniform Pauli magnetic susceptibility is given by:

[TABLE]

where $M_{tot}^{z}$ is the total magnetic moment of the site (unit cell):

[TABLE]

expressed through the density matrix elements between the local states.

The band structure along the path connecting the high symmetry points of the Brillouin zone as well as the DOS are both obtained from the spectral density $A_{l}(k,\epsilon)$:

[TABLE]

so that the total DOS: $\sum_{k,l}A_{l}(k,\epsilon)$, the partial one: $\sum_{k}A_{l}(k,\epsilon)$ and the band structure (as a set of points in $\epsilon,k$ plane where $\sum_{l}A_{l}(k,\epsilon)$ diverges) can be obtained.

Closing the discussion of the observables, we note that the solutions of self-consistency equations have their domains of existence as functions of Hamiltonian parameters, filling, temperature and any other external condition. These domains form the system’s phase diagrams, which can be directly compared with the experiment and/or other theories etc. As well, there might exist (and this is the standard situation) several solutions or phases for the same values of the Hamiltonian parameters, and/or temperature etc. In this case, owing to the variational principle, the phase with the lowest free energy must be chosen, which paves the way to the phase transitions at finite temperature.

3 IMPLEMENTATION DETAILS

The general energy optimization procedure as applied to the electronic phases of a crystal can be in principle solved in two ways:

as a self-consistency procedure: starting as described above from a trial $\rho_{old}$ and iteratively calculating new $\rho_{new}$ at each step finding the chemical potential. In order to stabilize such a procedure usually some forms of mixing are used. Here we use the simplest linear mixing, such that $\rho_{new}=\alpha\rho_{new}+(1-\alpha)\rho_{old}$ with $0\leqslant\alpha\leqslant 1$; 2. 2.

the self-consistency condition eq.(9) supplemented with the constraint fixing the number of particles can be also solved by using the globally convergent Newton-Raphson method[8] for root finding. It can be shown, that the self-consistency condition is equivalent to considering the evanescence conditions for the derivatives of the energy with respect to the matrix elements of the density. Finally, the energy is compared at each step to ensure the the chosen trial step actually leads to a minimum.

Both approaches are implemented in $\Theta\Phi$. In the second case it is crucial to use all possible symmetries of $\rho$ in order to reduce the number of unknowns. This will be highlighted in the following subsection. As will become evident from the following subsection, the density matrix contains all the (order) parameters necessary to define the trial wave function eq.(5).

3.1 Broken symmetry and generalized density matrix

In this subsection we briefly discuss the properties of the generalized density matrix $\rho$. This is the central object of mean-field theory which determines the electronic phase of the system. Variational wave function (5) depends on a set of variational parameters $\left\{\zeta\right\}$. On the other hand, the density matrix $\rho$, defined as a matrix of pairwise averages of all fermionic operators present in the Nambu vector eq.(4), also depends on $\left\{\zeta\right\}$. The number of independent parameters $\left\{\zeta\right\}$ is in general much smaller than the number of elements in $\rho$. However, analyzing $\rho$ helps to give a physical meaning to the variational parameters $\left\{\zeta\right\}$ as well as to determine the electronic phase of the system. Moreover, by triggering some of $\left\{\zeta\right\}$ it is possible to induce a phase change in the system. The wave function eq.(5) allows for two solutions breaking fundamental symmetries (and their combination): i) particle number conservation and ii) time-reversal invariance. The former corresponds to the superconducting (BCS) states, while the later - to all kinds of magnetic states (when supplied with the spin quantization axis rotation). In this context, we call the inter-relations among the $\rho$’s matrix elements (or constraints on their values)

• a particular density matrix symmetry. These symmetries should not be confused with the spatial or point group ones. Typically, the so-called anomalous terms in $\rho$ (those of the form $\left\langle c_{\uparrow}^{\dagger}c_{\downarrow}^{\dagger}\right\rangle$ or $\left\langle c_{\downarrow}^{\phantom{\dagger}}c_{\uparrow}^{\phantom{\dagger}}\right\rangle$) if non-vanishing lead to BCS states, while the imbalance between “up” and “down” spin averages ($\left\langle c_{\uparrow}^{\dagger}c_{\uparrow}^{\phantom{\dagger}}\right\rangle\neq\left\langle c_{\downarrow}^{\dagger}c_{\downarrow}^{\phantom{\dagger}}\right\rangle$) leads to magnetically ordered phases. Finally, we note that $\rho(\tau=0)$ is in general Hermitian, while at a finite lattice shifts this is not true any more and the following relation holds:

[TABLE]

i.e. Hermitian conjugate of a density matrix at the lattice shift $\tau$ is the density matrix at the lattice shift $-\tau$ assuring the hermiticity of the density matrix as a whole.

To be more specific, we adapt the general definition eq.(10) to the case of a two-orbital system where we allow for all possible types of pairing. The density matrix at zero distance in that case can be written as:

[TABLE]

Here the parameters $n_{1(2)\sigma},\Delta_{1(2)},\Delta_{12},r_{12\sigma}$ have the meaning of local occupation, local and inter-orbital pairing and inter-orbital bonding respectively. The density matrix at $\tau\neq 0$ has the same structure, except for the fact that on the diagonal there would be the inter-site intra-orbital hopping: $n_{1(2)\sigma}\to h_{1(2)\sigma}$ and all the other parameters will be sampled at the corresponding distance $r$. The zero matrix elements of $\rho$ stem from the definition of the trial wave-function (5) and from the commutation relations for fermionic operators. Already at this step it is clear that a $16\times 16$ matrix actually has only $9$ independent matrix elements. This symmetry is checked and enforced during $\Theta\Phi$ run, using iterative self-consistent scheme, while the reduction from $\rho(\tau=0)$ to parameters $n_{1(2)\sigma},\Delta_{1(2)},\Delta_{12},r_{12\sigma}$ and back is used in Newton-Raphson method in order to reduce the number of self-consistency equations and unknowns.

The switching between different phases in eq.(5) is accomplished by zeroing some of the parameters. E.g. setting all $\Delta=0$ brings the system into a Fermi-liquid phase with no superconducting fluctuations. Setting $\Delta_{1(2)}=0$ excludes local pairing, setting $\Delta_{12}$ on the contrary forbids the inter-orbital pairing. The well-known $d_{x^{2}-y^{2}}$-wave pairing can be achieved on a square lattice by imposing that the superconducting part of the density matrix along e.g. $y$ direction be equal to minus that along $x$ direction, which in case of a single orbital model can be symbolically written as:

[TABLE]

Of course, in case of multi-orbital systems and complicated lattices even more involved inter-relations between $\rho$ matrix elements at various distances (see e.g. Ref.[9]), may exist which are advisable to be analyzed before using $\Theta\Phi$.

3.2 Matrix elements and porting with existing programs

Thus, $\Theta\Phi$ requires the hoping matrix elements $t_{l,l^{\prime}}^{s,s^{\prime}}(\tau)$ as well as $U_{l,l^{\prime}}^{s,s^{\prime}}$,$V_{l,l^{\prime}}(\tau)$, and $J_{l,l^{\prime}}^{\alpha\beta}(\tau)$, as input. While the user can specify all the Hamiltonian parameters manually, assuming an abstract theoretical model, we additionally introduced a functionality allowing the import of the hopping matrix elements $t_{l,l^{\prime}}^{s,s^{\prime}}(\tau)$ from existing first-principles codes. This is achieved by reading the tight-binding Hamiltonians, generated by either of the two well-known post-processing codes: wannier90 [10] and lobster [11, 12, 13]. By doing this, the precise, material-specific predictions can be made. $\Theta\Phi$ will read the files called wannier90_hr.dat (for wannier90) or RealSpaceHamiltonians.lobster (for lobster) and populate the hopping matrix, superseding any previously stored. Both wannier90 and lobster have in turn interfaces to all major ab-initio codes (like vasp, abinit, castep, QuantumEspresso etc.), which allows the $\Theta\Phi$ user to profit from the results obtained by all these codes. Additionally, we can recommend using the MagAîxTic package [14] to obtain the required estimates of the$J_{l,l^{\prime}}^{\alpha\beta}(\tau)$ parameters.

4 EXAMPLES

The present examples Section pursues a dual purpose: to illustrate the capacities of the proposed package and by this to contribute to bridging certain gap between solid state quantum chemistry and solid state theoretical physics communities. For the former, the variety of the electronic states/phases considered and studied by the latter remains to a significant extent inaccessible together with the fascinating properties of these phases since they are simply not present in the tools (solid state quantum chemistry codes) the former use. For the latter, incidentally, remains inaccessible the possibility to independently estimate the parameters (matrix elements) of the model Hamiltonians featuring the most interesting phases/hottest topics of experimental and theoretical interest. Thus in this Section we start with simple example of the 1D Hubbard model.

4.1 Superconducting state in 1D Hubbard model

We start by considering the simplest model for strongly correlated electrons on a 1D ring of sites, each holding one orbital - single-orbital Hubbard model, which has the Hamiltonian:

[TABLE]

Here, we denote by $r$ the coordinate of the current cell (site), and we treat the interaction term according to the prescriptions outlined in Appendix A. If $U<0$, it is possible to stabilize a simple BCS superconducting solution, thanks to the electron-electron attraction term in the BCS Hamiltonian: $U\langle c_{r\uparrow}^{\dagger}c_{r\downarrow}^{\dagger}\rangle\langle c_{r\downarrow}^{\phantom{\dagger}}c_{r\uparrow}^{\phantom{\dagger}}\rangle$. The BCS state corresponds to the density matrix of the form:

[TABLE]

Here the components of the density matrix are: $n=\left\langle c_{rs}^{\dagger}c_{rs}^{\phantom{\dagger}}\right\rangle$, $\Delta=\langle c_{r\uparrow}^{\dagger}c_{r\downarrow}^{\dagger}\rangle$ being respectively the site population and the anomalous average. In the case of an on-site Hubbard attraction, only $\rho$ at zero lattice shift enters into the self-consistency equations, and the above averages do not depend on $r$ due to the translational invariance.

The $\Theta\Phi$ results for this simple model eq.(18) are shown in Fig.(1). As expected from the BCS theory, the temperature dependence of the a superconducting order parameter $\Delta$is described with a very good accuracy by:

[TABLE]

in the vicinity of the transition temperature $T_{c}$. This is very remarkable feature provided by the \mathit{\Theta}$$\mathit{\Phi} program - the temperature dependence of the solution of the electronic problem. Incidentally, the $U$-dependence of $\Delta_{0}$ at lowest temperature is very well described by a BCS-like formula:

[TABLE]

for small absolute values of *negative i.e. *attractive $U$.111We note that this simple variational ground state does not represent adequately the true ground state of the 1D Hubbard model. Indeed, by using the bosonization technique[15], it was demonstrated that in a wide range of 1D models with short-ranged interactions there occurs the separation of spin and charge degrees of freedom, while the correlation functions decay with a power-law as a function of distance. Thus the more accurate treatment reveals no long range correlations in 1D systems, although a sizable superconducting correlations are also predicted by bosonization. By this we demonstrate the workability of two important features of $\Theta\Phi$ - the accessibility of the BCS state of the system and the temperature dependence of the solution of an electronic problem.

4.2 Superconducting states of Graphene

Superconductivity in pristine graphene does not occur naturally, although recently the superconductivity has been reported in graphene bilayers rotated at a series of “magic” angles[16, 17]. On the other hand, combination of Dirac-like linear spectrum with superconductivity may lead to a bunch of novel phenomena, such as specular Andreev reflection [18]. Moreover, recently an artificial honeycomb 2D super-lattice - effectively equivalent to graphene layer - has been realized by using the self-assembly of PbSe, PbS or CdSe nano-crystals matching their (100) facets [19]. Such super-lattices, possess well defined Dirac spectrum, although with several new features arising from stronger spin-orbit coupling [20, 21]. Several approaches predicted superconductivity with various order parameter symmetries on 2D honeycomb lattice, among which auxiliary-boson theory (ABF) [22, 9] (also called the Resonating Valence Bond - RVB - state in this context), renormalization-group[23], or even self-consistent Bogoliubov-de Gennes approach to spatially inhomogeneous superconductivity [24, 23].

Here, we benchmark $\Theta\Phi$ by reproducing several ABF superconducting solutions of various symmetries in graphene and by comparing them with the reference results [22, 9]. We remind that on 2D honeycomb lattice each lattice site connects to three nearest neighbor (NN) sites by vectors $\{\delta_{i}\},\;i=1,2,3$. These vectors can be defined in several ways, and here we adopt the following definition (see Fig. 2):

[TABLE]

The unit cell of graphene contains two atoms and is defined by the following unit cell vectors: $a_{1}=\delta_{1}-\delta_{3}$ and $a_{2}=\delta_{2}-\delta_{3}$. If the coordinates of the first site in the unit cell are chosen to be zero, then the coordinates of the second one will be $u=\delta_{1}+\delta_{2}-\delta_{3}$. With only NN hopping active, each unit cell is connected to its four nearest neighbors at lattice shifts $\tau=\pm a_{1},\pm a_{2}$.

Superconducting phases on honeycomb lattice can come out with several symmetries: extended $s$-wave, $d_{x^{2}-y^{2}}$-wave, $d_{xy}$-wave and several triplet pairing symmetries (like $p$-wave); corresponding to different fixed relations between the superconducting order parameters $\langle c_{r\uparrow}^{\dagger}c_{r+\delta_{i}\downarrow}^{\dagger}\rangle$. In the present work, we consider only the extended $s$-wave, $d_{x^{2}-y^{2}}$-wave, $d_{xy}$-wave and report the $k$-space dependence of their gap functions in Tab. 1.

The Hamiltonian of free electrons, hopping at NN distance on honeycomb lattice, reads as follows:

[TABLE]

where $t$ is the hopping matrix element. It can be readily diagonalized, yielding the well-known graphene spectrum:

[TABLE]

It manifests the Dirac cones at $K=\pi(\frac{1}{3},\frac{2}{3})$ and $K^{\prime}=\pi(\frac{2}{3},\frac{1}{3})$ points of the Brillouin zone.

In Ref. [22], it was argued that a sizable exchange correlation $J$ in graphene could arise from on-site Hubbard repulsion $U$: $J=\frac{2t^{2}}{U}.$ We remind the value of $t$ known for long time in graphene: $t=2.8$ eV, (see [25, 26] and reference therein) while $U$ estimated in Ref. [22, 9] is $U=3.3t$. These estimates give for the ratio $J/t\approx 1.7$.

With the exchange term the full interacting Hamiltonian becomes:

[TABLE]

where the spin operators are defined through the fermionic operators according to eq.(3).

We treat the Hamiltonian eq.(21) within the mean-field approach in $\Theta\Phi$ according to the formalism outlined in the previous Section2 and compare our results with those of Refs. [22, 9]. In the case of graphene, which has two atoms in the unit cell (site $A$ and site $B$), we set $L=2$. Therefore, the density matrix in real space is a $8\times 8$ matrix for each distance. Here we consider the $t-J$ model on honeycomb lattice in which all the terms (both kinetic and interaction parts) have the range extending to at most NN inter-cell lattice shift $\tau=a_{1(2)}$, as shown in Fig.2. In addition, there is always present the local density matrix at $\tau=0$. For the sake of illustration, we show below local part of $\rho$:

[TABLE]

where various parameters appearing inside eq.(22) are $\tau=0$ limit of the following correlation functions:

[TABLE]

Here, $i=1,2$ is the orbital index, $\tau$ is the given unit cell separation, $s$ is the spin index. In addition, we consider only the non-local superconducting averages.

In the case of the extended-$s$ phase, all $\Delta_{12}$ at lattice shifts $\tau=0$, $\tau=a_{1}$ and $\tau=a_{2}$ are equal, while in the $d_{x^{2}-y^{2}}$ case: $\Delta_{12}=0$ at $\tau=0$ and $\Delta_{12}(\tau=a_{1})=-\Delta_{12}(\tau=a_{2})$, finally in the $d_{xy}$ case: $\Delta_{12}(\tau=0)=-2\Delta_{12}(\tau=a_{1})=-2\Delta_{12}(\tau=a_{2})$. $\sum_{i,\sigma}n_{i\sigma}(\tau=0)=\mathscr{\mathcal{N}}$ fixes the chemical potential by imposing the average occupation per unit cell to be $\mathcal{N}=2$. The angular dependence of the superconducting gap in $k$-space for the three symmetry cases[9] is summarized in Table1.

A few words should be said about the treatment of the exchange term in $\Theta\Phi$ and in the reference articles. Namely, in the former we retain all non-superconducting terms in the mean-field decoupling of the Hamiltonian term proportional to $J$ (see A), while in the latter the hopping renormalization, arising from the linearized exchange term was not taken into account. That is to say, discarding the last term in eq.(28). Moreover, when making the mean-field average, the coefficient of $\frac{3}{2}$ in front of $J$ was omitted too. In order to compare our results with those of Ref. [9], we rescale $J$ down to a factor of $\frac{2}{3}$ and discard the hopping renormalization in this subsection. With this assumptions, the agreement with the reference data of Ref. [9] is excellent, as can be seen from the right panel of Fig.2, where we compare the critical superconducting temperature $T_{c}$ as a function of doping $\eta$ (per C atom) at $J=1$ in the notations of Ref. [9]. This agreement holds for both extended-$s$ and $d_{x^{2}-y^{2}}$-wave phases (we did not found the data for $T_{c}$ in $d_{xy}$-wave phase). We would like to emphasize that each point in the right panel of Fig.2 is a result of a graph similar to the left panel of Fig. 1 with many points, each being a solution of self-consistency equations (LABEL:eq:self-cons). This proves the solidity of our benchmark.

To have an insight into the electronic structure of these two solutions, we compare the density of states (DOS) and the band structure in the extended-$s$ and $d_{x^{2}-y^{2}}$-wave phases, as shown in Fig.3. One can see that in the former case, there is a full gap in the DOS, while in the latter case DOS goes to zero linearly in the vicinity of $\omega=0$ as expected for $d_{x^{2}-y^{2}}$-symmetry gap. In the $k$-space, the true gap opens at $K$ (and $K^{\prime}$) points for the $s$-wave phase, while for the $d_{x^{2}-y^{2}}$-wave case there are only two nodal points on the line $M\to K\to\Gamma^{\prime}$ whose exact position depends on doping $\eta$ and $J$. Notice that in the latter case, there appears a difference between $K$ and $K^{\prime}$ as well as between $M=\pi(\frac{1}{2},\frac{1}{2})$ and $M^{\prime}=\pi(\frac{1}{2},0)$. Peculiarly, in the latter case the linear spectrum, intrinsic to $H_{0}$ is restored in the $d_{x^{2}-y^{2}}$-wave solution, although bearing a different meaning of Cooper pair nodal quasi-particles. In addition, the overall appearance of the band structure is very different in the two cases: in the extended-$s$ case $k$-dependence of the bands is that of the non-interacting dispersion (although renormalized by a coefficient proportional to $J$), which is also the gap function in this case, while in the $d_{x^{2}-y^{2}}$-wave case the overall $k$-dependence is a result of both non-interacting dispersion and the corresponding gap function. Such a big difference in the band structure leads to a sizable change in the total free energy among the phases. In Fig.4, we show the comparison of the free energy as a function of $J$ the three superconducting phases. The extended-$s$ phase is the lowest one at $J>0.75t$ (in the notations of Ref.[22]), while for $J<0.75t$ the $d$-wave phases (almost indistinguishable between each other in free energy) are most favorable.

We emphasize that the symmetry classification, derived in Ref.[22], is only rigorous in the proximity of the transition temperature, since the linearized self-consistency equations were used. Another validity condition for this classification is neglecting the band renormalization as explained above. Treatment of the full interaction term would introduce additional non-superconducting order parameters and would complicate significantly the symmetry analysis. Under the conditions of Ref. [22], we have found that for the extended-$s$ and $d_{x^{2}-y^{2}}$ phases the respective symmetry relations hold exactly in the whole range of their existence, while for the $d_{xy}$ phase the symmetry relations are only asymptotic in the limit $T\to T_{c}$. By this we show the workability of the $\Theta\Phi$ for a multi-band model and accessibility of the BCS states of different symmetry in it together with their temperature dependence.

4.3 Magnetic phases of iron

The magnetism of iron is of vital fundamental and technological importance and was extensively studied (see e.g. Ref. [28] for review and references therein as well as Ref.[29, 30]). Among various iron phases we chose $\alpha$- and $\gamma$-Fe in order to benchmark $\Theta\Phi$. In $\alpha$-Fe we compared the FM phase as obtained in $\Theta\Phi$ and DFT+U calculations in VASP, while in $\gamma$-Fe after comparing the NM and FM phase with VASP DFT+U, with use of $\Theta\Phi$ we addressed the AFM phase as well and showed that the system in this phase gains the energy with respect to the FM one.

We have performed DFT simulations using vasp [31] in the PBE/GGA [32] version. The energy cutoff was set at 500 eV. In the DFT+U calculations, $U=9$ eV, $J=1$ eV within Dudarev’s formulation were used[33].

The effective tight-binding Hamiltonian was obtained by wannier90 [10], with the frozen window method for disentanglement with inner window of $[-2.61:21]$ eV and the outer one of $[-2.61:70]$ eV. Local orbitals of $s,p$ and $d$ character were included and the final Hamiltonian, expressed in the maximally localized Wannier orbitals basis was used as input for $\Theta\Phi$. The Hubbard-$U$ correction terms where applied to Fe $d$ orbitals both in DFT and $\Theta\Phi$ settings. We have used the same values of $U$ and $J$ for comparison with vasp. In order to tackle the problem of double counting, typical in post-DFT approaches, we implemented in $\Theta\Phi$ the modified Coulomb term of the form [34]:

[TABLE]

which has both $U_{ll^{\prime}}$ and $J_{ll^{\prime}}$ terms orbital dependent while $\left\langle n_{l,r,\sigma}^{0}\right\rangle$ are average occupations of the “correlated” orbitals as coming from DFT. Rewritten this way, the Hubbard $U$ term describes the electronic fluctuations around the DFT values, without shifting the center of mass of a “correlated” band as a whole. From the practical side, the interaction term eq.(23) can be easily “assembled” from the elemental bits defined in A. For the sake of comparison with VASP we used isotropic $U_{jj^{\prime}}=U$ and $J_{jj^{\prime}}=J$, although $\Theta\Phi$ can handle arbitrary interaction matrices.

A few words should be said about which terms contribute to the Wick’s theorem decoupling of $H_{U}$ in $\Theta\Phi$ as compared to DFT+U. In the former, we retain all the terms in eq.(27), while in DFT+U only the first, “decoupled” term is retained [35, 36]. Even if the ground state without “anomalous” averages is considered, there are always sizable terms, contributing to the inter-orbital hopping, missing in traditional DFT+U. Our treatment, therefore, is more complete and consistent.

For $\alpha$-Fe we used the bcc lattice with $a=2.858$Å  and a $8\times 8\times 8$ Monkhorst-Pack $k$-point mesh [37] (both in $\Theta\Phi$ and VASP). In this case we considered a ferromagnetic solution (FM) and a non-magnetic solution (NM) and used the energy of the latter as a reference point for comparing the relative stability of the former. It can be seen from Table 2 that the agreement between $\Theta\Phi$ and VASP is in general good: FM state gains the energy with respect to the NM one and the $d$-shell magnetic moment is very similar. The energy gain in $\Theta\Phi$ is approximately $30\%$ greater, due to correct treatment of the intra-atomic orbital hopping terms as explained above. The comparison of the band structures is presented in Figure5.

For $\gamma$-Fe we considered the fcc lattice with $a=3.583$Å  and a $8\times 8\times 8$ Monkhorst-Pack $k$-point mesh [37]. In addition to the FM and NM states defined above, in this case we also consider an antiferromagnetic state (AFM) with the pitch vector $\mathbf{Q}=\pi(0,1,1)$ in terms of the vectors of the reciprocal lattice given relative to the primitive cell, which corresponds to an alternation of oppositely magnetized ferromagnetic planes along the $z$-axis of the conventional cubic unit cell.

We found that in $\gamma$-Fe the NM state has the highest energy among all the states considered. It can be seen from Table 4.3, that among the magnetic states, the AFM spiral state with the above value of $\mathbf{Q}$ is lower than the FM one by $1.6$ eV, which is larger than the corresponding VASP figure ($0.6$ eV). The FM state is $1.8$ eV higher than the NM one in $\Theta\Phi$, as compared to VASP, while for the AFM state this offset amounts to $0.8$ eV. This follows that the AFM state description in $\Theta\Phi$ is somehow closer to VASP as compared to the description of the FM state. Such a difference in energetics between $\Theta\Phi$ and VASP can be ascribed to the different treatment of the Hubbard term mean-field decoupling as explained above. Finally, the $d$-orbital magnetic moment in $\Theta\Phi$ is comparable to the VASP one. The very good comparison of the band structure results for FM $\gamma$-Fe is presented in Figure 6.