Development of $ab~{}initio$ method for exciton condensation and its application to $\bf TiSe_{2}$

The ab initio approach adopted in this work is based on the multicomponent DFT [61]. By generalizing the electronic density, this approach has opened the gate to computing molecule bonding length, magnetization, and superconductivity from first principles [61, 62, 63]. Starting from the many-electron interacting Hamiltonian

$$\hat{H}=\hat{T}+\hat{U}^{\rm ee}+\hat{V}_{{\rm ext}}$$

(1)

where the $\hat{T}$ is the electronic kinetic energy, $\hat{U}^{\rm ee}$ is the Coulomb interaction among electron gas, and $\hat{V}_{\rm ext}$ is the external field including the potential from the nuclei. The Kohn-Sham (KS) method [64] replaces the many-body interaction by a local XC potential and writes down an effective single-particle Hamiltonian:

$$\hat{H}^{{\rm KS}}=\hat{T}[n]+\hat{V}_{{\rm XC}}[n]$$

(2)

which can reproduce the same ground state energy as Eq. (1) with the same local electronic density according to the Hohenberg-Kohn theorem [65]. Thus, if we have the true XC potential, we can obtain the exact solution for the target system, but $\hat{V}_{\rm XC}$ can only be acquired approximately.
In correlated systems, like superconductor [63, 59], the local potential is insufficient to represent the many-body interaction, and NL potential and density must be adopted to characterize the electronic order. Here, we start from a DFT with the local charge density $n(\mathbf{r})$ and the NL charge densities $\chi(\mathbf{r},\mathbf{r}^{\prime})$ in the form:

	$\displaystyle n(\mathbf{r})=\sum_{\sigma}\langle\Psi_{\sigma}^{\dagger}(\mathbf{r})\Psi_{\sigma}(\mathbf{r})\rangle$
	$\displaystyle\chi(\mathbf{r},\mathbf{r}^{\prime})=\sum_{\sigma}\langle\Psi_{\sigma}^{\dagger}(\mathbf{r})\Psi_{\sigma}(\mathbf{r}^{\prime})\rangle$		(3)

where $\Psi_{\sigma}$ is the electron operator, $\sigma$ denotes the spin index, and $\langle A\rangle={\rm Tr}\left[e^{-\beta H}A\right]$ defines the thermal expectation value. Therefore, following the KS construction of DFT, we can include possible external potential acting on the NL charge density by considering the Hamiltonian:

$$\hat{H}=\hat{T}+\hat{U}^{\rm ee}+\hat{V}_{\rm ext}+\hat{\Delta}_{\rm ext}$$

(4)

where

	$\displaystyle\hat{V}_{\rm ext}=\sum_{\sigma}\int d^{3}{\mathbf{r}}\Psi^{\dagger}_{\sigma}(\mathbf{r})v_{\rm ext}(\mathbf{r})\Psi_{\sigma}(\mathbf{r})$
	$\displaystyle\hat{\Delta}_{\rm ext}=\sum_{\sigma}\int d^{3}{\mathbf{r}}\Psi^{\dagger}_{\sigma}(\mathbf{r})\Delta_{\rm ext}(\mathbf{r},\mathbf{r}^{\prime})\Psi_{\sigma}(\mathbf{r}^{\prime}).$		(5)

At finite temperature, this Hamiltonian reproduces the thermodynamic potential as a functional of local and NL charge density:

	$\displaystyle\Omega[n,\chi]=$	$\displaystyle F[n,\chi]+\int d^{3}{\mathbf{r}}~{}n(\mathbf{r})\left[v_{\rm ext}(r)-\mu\right]$
		$\displaystyle+\int d^{3}{\mathbf{r}}d^{3}{\mathbf{r}^{\prime}}\left[\chi(\mathbf{r},\mathbf{r}^{\prime})\Delta_{\rm ext}(\mathbf{r},\mathbf{r}^{\prime})+h.c.\right]$		(6)

where the free energy $F[n,\chi]$ is universal and independent from external potentials.
Compared to the interacting Hamiltonian, Eq. (4), we can write down the corresponding KS Hamiltonian of non-interacting orbitals [64]:

$$\hat{H}^{\rm KS}=\hat{T}+\hat{V}_{\rm XC}+\hat{\Delta}_{\rm XC}$$

(7)

and the thermodynamic potential becomes:

	$\displaystyle\Omega^{\rm KS}[n,\chi]$	$\displaystyle=F^{\rm KS}[n,\chi]+F^{\rm L}_{\rm XC}[n]+F^{\rm NL}_{\rm XC}[\chi]$
		$\displaystyle+\frac{1}{2}\int d^{3}{\mathbf{r}}d^{3}{\mathbf{r}^{\prime}}\frac{n(\mathbf{r})n(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}-\mu\int d^{3}{\mathbf{r}}n(\mathbf{r})$		(8)

where $F^{\rm KS}[n,\chi]$ is the free energy of non-interacting KS system, and $F^{\rm L}_{\rm XC}[n]$ and $F^{\rm NL}_{\rm XC}[\chi]$ are the XC free energy from the local and NL densities. To make the KS Hamiltonian Eq. (7) reproduce the same thermal ground state as the interacting system, both thermodynamic potentials must be minimized by the same ground state densities:

		$\displaystyle\frac{\partial\Omega}{\partial n(\mathbf{r})}=\frac{\partial F}{\partial n(\mathbf{r})}+v_{\rm ext}(\mathbf{r})-\mu=0$
		$\displaystyle\frac{\partial\Omega^{\rm KS}}{\partial n(\mathbf{r})}=\frac{\partial F^{\rm KS}}{\partial n(\mathbf{r})}+\frac{\partial F^{L}_{\rm XC}}{\partial n(\mathbf{r})}+\int d^{3}{\mathbf{r}^{\prime}}\frac{n(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}-\mu=0$		(9)

and

		$\displaystyle\frac{\partial\Omega}{\partial\chi(\mathbf{r},\mathbf{r}^{\prime})}=\frac{\partial F}{\partial\chi(\mathbf{r},\mathbf{r}^{\prime})}+\Delta_{\rm ext}(\mathbf{r},\mathbf{r}^{\prime})=0$
		$\displaystyle\frac{\partial\Omega^{\rm KS}}{\partial\chi(\mathbf{r},\mathbf{r}^{\prime})}=\frac{\partial F^{\rm KS}}{\partial\chi(\mathbf{r},\mathbf{r}^{\prime})}+\frac{\partial F^{\rm NL}_{\rm XC}}{\partial\chi(\mathbf{r},\mathbf{r}^{\prime})}=0.$		(10)

By defining the exchange-correlation potential:

	$\displaystyle v_{\rm xc}[n](\mathbf{r})\equiv\frac{\partial F^{L}_{\rm XC}}{\partial n(\mathbf{r})}$
	$\displaystyle\Delta_{\rm xc}[\chi](\mathbf{r},\mathbf{r}^{\prime})\equiv\frac{\partial F^{\rm NL}_{\rm XC}}{\partial\chi(\mathbf{r},\mathbf{r}^{\prime})},$		(11)

we can write down the generalized KS equation:

		$\displaystyle\left[-\frac{\nabla^{2}}{2}+\int d^{3}{\mathbf{r}^{\prime}}\frac{n(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}+v_{\rm xc}-\mu\right]\psi_{i}(\mathbf{r})$
		$\displaystyle~{}~{}~{}~{}~{}+\int d^{3}{\mathbf{r}^{\prime}}\Delta_{\rm xc}(\mathbf{r},\mathbf{r}^{\prime})\psi_{i}(\mathbf{r}^{\prime})=E_{i}\psi_{i}(\mathbf{r})$		(12)

The result shows that if we can obtain the exact form of the XC free energy, we can construct the exact ground state of the interacting system based on the KS theorem. This derivation is general for all NL-densities beyond the definition in Eq. (3). For example, by choosing $\chi(\mathbf{r},\mathbf{r}^{\prime})=\langle\Psi_{\uparrow}(\mathbf{r})\Psi_{\downarrow}(\mathbf{r}^{\prime})\rangle$, i.e. the Cooper pairing, we can obtain the DFT for superconductivity [63, 59].
In this work, we apply the KS perturbation theory [66] following the similar approach of SCDFT [58, 59] to determine the $F_{\rm XC}$’s. By identifying the KS Hamiltonian, Eq. (7), as the unperturbed part from the interacting Hamiltonian, Eq (4), we obtain a perturbation theory:

$$\hat{H}_{0}=\hat{H}^{\rm KS}~{}~{};~{}~{}\hat{H}_{I}=\hat{U}^{\rm ee}-\hat{V}_{\rm XC}-\hat{\Delta}_{\rm XC}.$$

(13)

Provided with all order contributions in the $\hat{\Delta}_{\rm XC}$ operator, the difference between the interacting and KS thermodynamic potential can be expanded by the linked-cluster theorem [67] as:

$\displaystyle\Omega=\Omega^{\rm KS}-\frac{1}{\beta}\sum_{l=1}^{\infty}U_{l}$

(14)

where $U_{l}$ are different connected Feynman diagrams. According to Eq. (II.1) and Eq. (II.1), both thermal potentials are minimized by the ground state density so that the derivative of Eq. (14) must vanish and provide the condition to solve for the XC potential. In Appendix A, we show that this method can also be applied to derive the GW correction, which is commonly obtained from the Green’s function method [68].

To discuss the structural transition, we must adopt a theory involving the atomic degree of freedom. In this work, we restrict the discussion to the Born-Oppenheimer approximation which uses the phonon vibration to represent the atomic motion. The standard approach expands the interatomic potential from the equilibrium structure in the Taylor series, where the lowest quadratic term defines the non-interacting phonon normal mode while higher-order terms take account of the anharmonic phonon-phonon (ph-ph) interactions [69]. As a result, the spontaneous lattice distortion can be characterized by the thermal average of the phonon operator, which is proportional to the static atomic displacement.
To accommodate the lattice deformation, we extend the KS Hamiltonian by introducing a phonon Hamiltonian:

$$\hat{H}_{\rm ph}=\hat{H}_{\rm 0,ph}+\hat{U}^{\rm ph}_{\rm anh.}+\hat{U}^{\rm e-ph}+\hat{\Delta}_{b}$$

(15)

where the free phonon Hamiltonian $\hat{H}_{\rm 0,ph}$, ph-ph anharmonic interaction $\hat{U}^{\rm ph}_{\rm anh.}$, and e-ph interaction $\hat{U}^{\rm e-ph}$ are defined as standard forms found in literature [60, 70]. The last term denotes a linear deformation potential of the form:

$$\hat{\Delta}_{b}=\sum_{\mathbf{q}\nu}\left(\Delta_{b,\mathbf{q}\nu}^{*}\hat{b}_{\mathbf{q}\nu}+\Delta_{b,\mathbf{q}\nu}\hat{b}^{\dagger}_{\mathbf{q}\nu}\right),$$

(16)

which can push atoms away from their equilibrium position without phonon softening. Therefore, the phonon operator acquires a finite thermal average:

$$\langle\hat{b}_{\mathbf{q}\nu}\rangle=\chi_{b,\mathbf{q}\nu}$$

(17)

such that if we treat $\chi_{b,\mathbf{q}\nu}$ on the same footing as the densities in Eq .(3), the $\Delta_{b}$’s can be determined by the perturbation expansion method introduced in the previous section where the Hamiltonian now includes the phonon degree of freedom:

	$\displaystyle\hat{H}=\hat{H}_{0}+\hat{H}_{\rm I}$		(18a)
	$\displaystyle\hat{H}_{0}=\hat{T}+\hat{V}_{\rm XC}+\hat{\Delta}_{\rm XC}+\hat{H}_{\rm 0,ph}+\hat{\Delta}_{b}$		(18b)
	$\displaystyle\hat{H}_{\rm I}=\hat{U}^{\rm ee}-\hat{V}_{\rm XC}-\hat{\Delta}_{\rm XC}+\hat{U}^{\rm ph}_{\rm anh.}+\hat{U}^{\rm e-ph}-\hat{\Delta}_{b}.$		(18c)

This is the general formalism to describe lattice deformation at finite temperature. The distorted pattern follows the phonon normal vector denoted by the mode index $\nu$, and the displacement amplitude is proportional to the absolute value of $\chi_{b,\mathbf{q}\nu}$. The momentum $\mathbf{q}$ will determine the periodicity of the new crystal structure. In Sec. III we will apply Eq. (18a-18c) to derive the gap equation for the exciton condensation.

In this section, we provide a discussion concerning the phonon properties obtained by the first-principles approaches. The ab initio phonon vibration can be calculated by diagonalizing the dynamical matrix, which encodes the total energy change when atoms move away from their equilibrium position. In practical approach, the dynamical matrix is computed by density functional perturbation theory (DFPT) [69] or using the frozen phonon method [71], which gives the phonon vibration at zero temperature. Phonon frequencies obtained by these methods contain the ”bare” part as well as the renormalization effect:

	$\displaystyle\omega^{2}_{\mathbf{q}\nu}={\omega^{\rm bare}_{\mathbf{q}\nu}}^{2}+2\omega^{\rm bare}_{\mathbf{q}\nu}\Pi^{\rm bare}_{\mathbf{q}\nu}(T)$		(19a)
	$\displaystyle\Pi^{\rm bare}_{\mathbf{q}\nu}(T)=\frac{2}{\mathcal{N}_{k}}\sum_{\mathbf{k}mn}|g^{\rm bare}_{mn\nu}(\mathbf{k},\mathbf{q})|^{2}\frac{f_{\mathbf{k}+\mathbf{q},m}-f_{\mathbf{k},n}}{\xi^{m}_{\mathbf{k}+\mathbf{q}}-\xi^{n}_{\mathbf{k}}}$		(19b)

where $g^{\rm bare}_{mn\nu}(\mathbf{k},\mathbf{q})$ is the bare e-ph matrix element:

$$g^{\rm bare}_{mn\nu}(\mathbf{k},\mathbf{q})=\left(\frac{\hbar}{2\omega_{\mathbf{q}\nu}^{\rm bare}}\right)^{1/2}\langle\psi_{m\mathbf{k}+\mathbf{q}}|\partial_{\mathbf{q}\nu}V|\psi_{n\mathbf{k}}\rangle,$$

(20)

and $\Pi^{\rm bare}_{\mathbf{q}\nu}$ is the self-energy due to the Coulomb screening, where the superscript ”bare” denotes that we used the bare coupling constant from Eq. (20). The factor of 2 in Eq. (19b) appears for the spin degeneracy. $\mathcal{N}_{k}$ is the number of k-points, $\xi^{n}_{\mathbf{k}}$ is the electron energy of the $n$-th band with momentum $\mathbf{k}$ measured from the chemical potential, $f_{\mathbf{k},n}=1/(e^{\beta\xi^{n}_{\mathbf{k}}}+1)$ is the Fermi-Dirac distribution function. Here we note that although the self-energy depends on the phonon frequency via the coupling constant, the renormalization only depends on the deformation potential $\partial_{\mathbf{q}\nu}V$ and the electronic occupations $f_{\mathbf{k},n}$. Therefore, for later use, we define a frequency-independent factor

$$X_{\mathbf{q}\nu}(T)=\frac{\hbar}{\mathcal{N}_{k}}\sum_{\mathbf{k}mn}|\langle\psi_{m\mathbf{k}+\mathbf{q}}|\partial_{\mathbf{q}\nu}V|\psi_{n\mathbf{k}}\rangle|^{2}\frac{f_{\mathbf{k}+\mathbf{q},m}-f_{\mathbf{k},n}}{\xi^{m}_{\mathbf{k}+\mathbf{q}}-\xi^{n}_{\mathbf{k}}},$$

(21)

such that the phonon frequency difference between two temperatures due to the screening effect can be written as

$$\omega^{2}_{\mathbf{q}\nu}(T_{1})-\omega^{2}_{\mathbf{q}\nu}(T_{2})=2X_{\mathbf{q}\nu}(T)\Big{|}^{T_{1}}_{T_{2}},$$

(22)

where we neglect the temperature dependence in the deformation potential.
In general, since experimental observables are the fully renormalized quantities, researchers focus on the explanation of the direct result $\omega$ as the experiment measurement, while the isolation between the bare frequency $\omega^{\rm bare}$ and self-energy is rarely studied [72]. However, when a single phonon momentum $\mathbf{q}$ connects the electron pocket and hole pocket in some metallic systems, $\Pi_{\mathbf{q}\nu}$ tends to diverge, which will force the corresponding renormalized phonon becomes soft with imaginary frequency since the self-energy is always negative. This so-called Fermi-surface nesting effect will destabilize the crystal structure, result in lattice distortion, and cause the transition into a CDW phase.
Phonons derived from the above method are limited to the harmonic approximation, which expands the energy variation to the second order of atomic displacement. Although the harmonic phonon approach has great success in studying transport properties, carrier relaxation, and polaronic system [73, 74, 75], it fails to describe many important properties related to the lattice anharmonicity, such as thermal expansion, lattice transition, the temperature dependence of phonon frequency, and prediction of superconductivity critical temperature [76, 77, 78, 79, 80, 81].
In this work, we adopt the self-consistent phonon (SCP) theory [82, 60] to compute the anharmonic phonon frequency. The SCP theory expands the energy to quartic order of the atomic displacement:

$$U_{n}=\frac{1}{n!}\left(\frac{\hbar}{2}\right)^{\frac{n}{2}}\sum_{\{q_{n}\}}\delta({\Sigma\mathbf{q}_{n},\mathbf{G}})\frac{\Phi(q_{1},\cdots q_{n})}{\sqrt{\omega_{q_{1}}\cdots\omega_{q_{n}}}}\hat{A}_{q_{1}}\cdots\hat{A}_{q_{n}}$$

(23)

where $n$ is the number of moved atom, $q=(\mathbf{q},\nu)$ is the collective index, $\delta({\Sigma\mathbf{q}_{n},\mathbf{G}})$ impose the momentum conservation, $\Phi(q_{1},\cdots q_{n})$ is the interatomic force constant, and $\hat{A}_{q}$ is the displacement operator. To the lowest order, the self-energy comes from the loop diagram of the 4-points vertex:

$$\Sigma_{\mathbf{q},\nu\nu^{\prime}}=-\frac{1}{2}\sum_{q_{1}}\frac{\hbar\Phi(\mathbf{q}\nu,-\mathbf{q}\nu^{\prime},q_{1},-q_{1})}{4\sqrt{\omega^{\rm scp}_{\mathbf{q}\nu}\omega^{\rm scp}_{-\mathbf{q}\nu^{\prime}}}\omega^{\rm scp}_{q_{1}}}\times(1+2n_{{\rm B}}(\omega^{\rm scp}_{q_{1}}))$$

(24)

where $n_{\rm B}$ is the Bose-Einstein distribution function. This anharmonic self-energy contribution contains the mixing between different phonon normal modes $\nu,\nu^{\prime}$ such that it modifies both the phonon frequency and the vibration pattern. When the off-diagonal term can be neglected, the SCP equation can be simplified as:

$\displaystyle{\omega^{\rm scp}_{\mathbf{q}\nu}}^{2}=\omega^{2}_{\mathbf{q}\nu}+2\omega^{\rm scp}_{\mathbf{q}\nu}\Sigma_{\mathbf{q},\nu\nu}$

(25)

It has been shown that the SCP equation requires a real solution for all $\Omega$’s such that the self-consistent solution breaks down when any phonon becomes soft.

We first focus on Fig. 1(a) and derive the basic structure of the gap equation. The mathematical representation for the diagram is straightforward:

	$\displaystyle U^{\rm(a)}=\sum_{cv\mathbf{k}}g_{cv\nu_{0}}(\mathbf{k},{\bf M})\chi^{cv}_{\mathbf{k}}\left(\chi_{b,{\bf M}\nu_{0}}+\chi^{*}_{b,{\bf M}\nu_{0}}\right)$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\Delta^{cv}_{\mathbf{k}}\chi^{cv}_{\mathbf{k}}-\Delta_{b}^{*}\chi_{b,{\bf M}\nu_{0}}+c.c.$		(36)

where we have used the property, $\hat{b}_{\bf M}=\hat{b}_{-{\bf M}}$. Applying the minimum condition in Eq. (14) and taking the derivative respective to the NL densities, by keeping $\chi_{b,{\bf M}}(\chi_{b})$ and $\chi_{b,{\bf M}}^{*}(\chi^{*}_{b})$ as independent variables, we can obtain:

	$\displaystyle\Delta_{b}=\sum_{cv,\mathbf{k}}g_{cv\nu_{0}}(\mathbf{k},{\bf M})\chi^{cv}_{\mathbf{k}}+g^{*}_{cv\nu_{0}}(\mathbf{k},{\bf M}){\chi^{cv}_{\mathbf{k}}}^{*}$		(37a)
	$\displaystyle\Delta^{cv}_{\mathbf{k}}=g_{cv\nu_{0}}(\mathbf{k},{\bf M})\left(\chi_{b,{\bf M}\nu_{0}}+\chi^{*}_{b,{\bf M}\nu_{0}}\right).$		(37b)

These two equations are the same gap equations derived from the mean-field approach [9], and we will refer them as ”MF gap equation” for the rest discussion. They can be solved consistently using the solution of Eq. (31a) and the definition in Eq. (17) and Eq. (27b). To present the detail, we carry out an example of an analytically solvable model in a one-dimensional system and left the discussion in Appendix C.

The mass renormalization effect can be separated into two parts corresponding to the gap function $\Delta_{b}$ and $\Delta_{\mathbf{k}}$, respectively, and here we focus on the $\Delta_{b}$. We isolate the $\bf q=M$ part (Fig. 1(c)) from the bubble diagram (Fig. 1(d)) and combine it with the contribution from the effective self-energy term. To compute Fig. 1(c), we note that due to the $\hat{\Delta}$ potential term, the phonon operator $\hat{b}^{(\dagger)}_{\bf M}$ in the unperturbed Hamiltonian, Eq. (31a), acquires a non-zero thermal average such that its thermal Green’s function becomes (see Appendix D):

$$D_{\nu}(\mathbf{q},i\omega_{n})=\frac{-2\Omega_{\mathbf{q}\nu}}{\omega^{2}_{n}+\Omega^{2}_{\mathbf{q}\nu}}-\delta_{\omega_{n},0}\delta_{\mathbf{q},{\bf M}}\delta_{\nu,\nu_{0}}\frac{\beta(\Delta_{b}^{*}+\Delta_{b})^{2}}{\Omega_{\mathbf{q}\nu}^{2}}.$$

(38)

The Green’s function now takes an additional static ($\omega_{n}=0$) term for the $\nu_{0}$-phonon with $\mathbf{q}={\bf M}$ which is quadratically proportional to the potential $\Delta_{b}$. Therefore, the $\Delta_{b}$-relevant part of Fig. 1(b,c) becomes:

	$\displaystyle U^{\rm(b+c)}$	$\displaystyle=$	$\displaystyle\frac{(\Delta_{b}^{*}+\Delta_{b})^{2}}{\Omega_{{\bf M}\nu_{0}}^{2}}\Pi_{{\bf M}\nu_{0}}(T)$		(39)
		$\displaystyle+$	$\displaystyle\frac{|\Delta_{b}|^{2}}{\Omega^{2}_{{\bf M}\nu_{0}}}\Bigg{(}\Sigma_{\mathbf{q},\nu\nu^{\prime}}(T)-\Pi_{\bf M\nu_{0}}(T=0~{}{\rm K})\Bigg{)}$

where $\Pi_{\bf M\nu_{0}}$ is defined as Eq. (19b) with the e-ph coupling replaced by Eq. (35). Using the relation Eq. (78) and taking the derivative respective to $\chi^{*}_{b}$, we can obtain:

	$\displaystyle\frac{\partial U^{\rm(b+c)}}{\partial\chi^{*}_{b}}$	$\displaystyle=\frac{-2\Delta_{b}}{\Omega^{2}_{{\bf M}\nu_{0}}}\bigg{(}2\Omega_{{\bf M}\nu_{0}}\Sigma_{{\bf M}\nu_{0}}(T)+2X_{{\bf M}\nu_{0}}(T)\Big{|}^{T}_{T=0}\bigg{)}$
		$\displaystyle\equiv-\mathcal{Z}_{b}\Delta_{b}$		(40)

The mass renormalization for the phonon displacement potential can thus be obtained by adding Eq. (III.2) to the right hand side of Eq. (37a), such that the gap equation is modified by a overall ratio:

$$\Delta_{b}\equiv\lambda_{b}\left[\sum_{cv,\mathbf{k}}g_{cv\nu_{0}}(\mathbf{k},{\bf M})\chi^{cv}_{\mathbf{k}}+g^{*}_{cv\nu_{0}}(\mathbf{k},{\bf M}){\chi^{cv}_{\mathbf{k}}}^{*}\right]$$

(41)

where $\lambda_{b}=1/(1+\mathcal{Z}_{b})$. Since the factor $\mathcal{Z}_{b}$ is always positive, the mass renormalization will reduce the displacement potential and thus decrease the critical temperature for exciton condensation.

The mass renormalization on the NL potential $\Delta_{\mathbf{k}}$ can be obtained from Fig. 1(d) which is the bubble diagram involving only normal propagators, and here we include all the phonon momentum $\mathbf{q}$ and normal modes $\nu$. To compute Fig. 1(d), we first note that for the state of Eq. (26), the Green’s function can be written as:

$$\hat{G}_{n(\mathbf{k})}(i\omega_{n})=\sum_{cv}\frac{|t^{n}_{v\mathbf{k}}|^{2}|{v\mathbf{k}}\rangle\langle{v\mathbf{k}}|+|t^{n}_{c\mathbf{k}}|^{2}|{c\mathbf{k}+{\bf M}}\rangle\langle{c\mathbf{k}+{\bf M}}|}{i\omega_{n}-E_{n(\mathbf{k})}}$$

(42)

which already mixes states of momentum $\mathbf{k}$ and $\mathbf{k}+M$ following the Hamiltonian Eq. (31a). Therefore, Fig. 1(d) can be expressed as:

	$\displaystyle U^{\rm(d)}$	$\displaystyle=$	$\displaystyle\sum_{\mathbf{k}\mathbf{q},nm,\nu}I(E_{m(\mathbf{k}+\mathbf{q}),E_{n(\mathbf{k})}},\Omega_{\mathbf{q}\nu})$
		$\displaystyle\times$	$\displaystyle\bigg{[}\sum_{v_{i}v_{j}}|t^{n}_{v_{i}\mathbf{k}}|^{2}|t^{m}_{v_{j}\mathbf{k}+\mathbf{q}}|^{2}|g_{v_{j}v_{i}\nu}(\mathbf{k},\mathbf{q})|^{2}$
			$\displaystyle+\sum_{c_{i}c_{j}}|t^{n}_{c_{i}\mathbf{k}}|^{2}|t^{m}_{c_{j}\mathbf{k}+\mathbf{q}}|^{2}|g_{c_{j}c_{i}\nu}(\mathbf{k}+{\bf M},\mathbf{q})|^{2}\bigg{]}$

where we define

$$I(E_{1},E_{2},\omega)=\sum_{\omega_{n}\omega_{m}}\frac{1}{i\omega_{n}-E_{1}}\frac{1}{i\omega_{m}-E_{2}}\frac{-2\omega}{(\omega_{n}-\omega_{m})^{2}+\omega^{2}}$$

(44)

which is the common expression appearing in the bubble diagram including two fermions and one boson propagator. However, in a general system with bands of more than two, there is no analytical expression of coefficient $t^{n}_{v\mathbf{k}}$’s in terms of coupling potential $\Delta_{\mathbf{k}}$. As a result, we take the small $\Delta_{\mathbf{k}}$ approximation and use the perturbation theory to expand the coefficients. In the linear order, since the eigenvalue is the same as the diagonal entries, we can use the original band index to denote the new states such that:

	$\displaystyle E_{n=v_{i}(\mathbf{k})}=\xi^{v_{i}}_{\mathbf{k}}~{};~{}t^{n=v_{i}}_{v_{j}\mathbf{k}}=\delta_{v_{i},v_{j}}~{};~{}t^{n=v_{i}}_{c_{j}\mathbf{k}}=\frac{\Delta^{v_{i}c_{j}}_{\mathbf{k}}}{\xi^{c_{j}}_{\mathbf{k}+{\bf M}}-\xi^{v_{i}}_{\mathbf{k}}}$
	$\displaystyle E_{n=c_{i}(\mathbf{k})}=\xi^{c_{i}}_{\mathbf{k}+{\bf M}}~{};~{}t^{n=c_{i}}_{c_{j}\mathbf{k}}=\delta_{c_{i},c_{j}}~{};~{}t^{n=c_{i}}_{v_{j}\mathbf{k}}=\frac{\Delta^{c_{i}v_{j}}_{\mathbf{k}}}{\xi^{v_{j}}_{\mathbf{k}}-\xi^{c_{i}}_{\mathbf{k}+{\bf M}}}$

However, the linear order is insufficient to express the eigenvector terms in $U^{\rm(d)}$ since the terms always appear as an absolute square, $|t^{n}_{c,v}|^{2}$. We compute the perturbation expansion to the second order for the explicit $\Delta^{cv}_{\mathbf{k}}$-dependence and take the derivative on Eq. (III.3). The algebraic detail is left to the Appendix E, and we quote the result here. The derivative of the bubble diagram respective to NL-densities can be written down as:

		$\displaystyle\frac{\partial U^{\rm(d)}}{\partial\chi_{\mathbf{k}}^{cv*}}=\frac{\Delta^{cv}_{\mathbf{k}}}{(\xi^{c}_{\mathbf{k}+{\bf M}}-\xi^{v}_{\mathbf{k}})\left(f_{\mathbf{k}+{\bf M},c}-f_{\mathbf{k},v}\right)}$
		$\displaystyle~{}~{}~{}~{}~{}\times\sum_{\mathbf{q},\nu}\bigg{[}\sum_{v_{2}}I(\xi^{v_{2}}_{\mathbf{k}+\mathbf{q}},\xi^{v}_{\mathbf{k}},\omega_{\mathbf{q}\nu})|g_{v_{2}v\nu}(\mathbf{k},\mathbf{q})|^{2}$
		$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{c_{2}}I(\xi^{c_{2}}_{\mathbf{k}+{\bf M}+\mathbf{q}},\xi^{c}_{\mathbf{k}+{\bf M}},\omega_{\mathbf{q}\nu})|g_{c_{2}c\nu}(\mathbf{k}+{\bf M},\mathbf{q})|^{2}\bigg{]}$
		$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\equiv-\mathcal{Z}^{cv}_{\mathbf{k}}\Delta^{cv}_{\mathbf{k}}.$		(46)

By adding Eq. (III.3) to the right-hand side of Eq. (37b), the gap equation is then modified by an overall ratio from the mass renormalization term:

$$\Delta^{cv}_{\mathbf{k}}\equiv\lambda^{cv}_{\mathbf{k}}g_{cv\nu_{0}}(\mathbf{k},{\bf M})\left(\chi_{b,{\bf M}\nu_{0}}+\chi^{*}_{b,{\bf M}\nu_{0}}\right)$$

(47)

where $\lambda^{cv}_{\mathbf{k}}=1/(1+\mathcal{Z}^{cv}_{\mathbf{k}})$. Due to the occupation difference in the denominator, $\mathcal{Z}_{\mathbf{k}}$ diverges and suppresses the potential $\Delta_{\mathbf{k}}$ for pairing between electron-electron($e$-$e$) and hole-hole($h$-$h$), such that the remaining configuration can only contain one electron and one hole ($e$-$h$ pairing).
In conclusion, Eq. (41) and Eq. (47) determine the gap which characterizes the excitonic insulator triggered by the e-ph coupling in a distorted lattice structure. The results contain the crucial mass renormalization factors for the electronic and atomic sectors, which are missed in the conventional MF approach. This derivation starts from first principles, requires no empirical parameters, and can be applied to study natural systems.

The first-principles electronic structure of monolayer $\rm TiSe_{2}$ can be found in the literature [44], and here we present our calculation as a starting point. We perform the DFT calculation on the relaxed hexagonal structure using the QUANTUM ESPRESSO package [84], with the Perdew–Burke–Ernzerhof exchange-correlation functional and fully relativistic norm-conserving pseudopotentials generated with Pseudo Dojo [85, 51, 86]. Going beyond DFT, we use the VASP code to include the GW correction in an additional calculation [87]. All the computational details are summarized in Appendix F.
In Fig. 2(a), we present the computed band structure where we match two results by their Fermi level at 0 eV and connect isolated data points by the Wannier interpolation method using the interface between VASP and wannier90 [88]. Both DFT and GW calculation shows that monolayer $\rm TiSe_{2}$ is a semi-metal with a hole pocket near the $\Gamma$ point and electron pocket in the M point while the Fermi surface is much smaller in the GW bands. By resolving the orbital in the band structure [44], we found that, with respect to the Fermi level, bands containing the Se p-orbitals tend to be suppressed to lower energy while the bands with the Ti d-orbitals can be promoted to higher energy in the GW correction. Near the $\Gamma$ point, two kinds of orbital are entangled, so the bands have Mexican hat-like dispersions. On the other hand, away from the $\Gamma$ point, the overall vertical stretch between occupied and unoccupied bands is about 1 eV.

In bulk $\rm TiSe_{2}$ a soft phonon with imaginary frequency has been observed and computed with wave vector $\mathbf{q}=L=(\frac{\pi}{2},\frac{\pi}{2},\frac{\pi}{2})$ which corresponds to the formation of $2\times 2\times 2$ CDW in three-dimensional crystal [31, 89, 90]. The soft phonon can be ”hardened” by increasing the environment temperature. However, in first-principles studies where occupation broadening was used for charge density to simulate the finite temperature effect, extremely large broadening is required for dynamical stability and unrealistic high transition temperature was obtained [89, 44] for both the bulk and monolayer structures. The failure of this approach indicates a strong anharmonicity in $\rm TiSe_{2}$ beyond the harmonic approximation, which is highlighted in a recent study [83].
Here we apply the method introduced in Sect. II.3 to compute the self-consistent anharmonic phonon frequency and study the phonon softening in monolayer $\rm TiSe_{2}$. In the practical procedure, we carry out the DFPT calculation using QUANTUM ESPRESSO to obtain the dynamical matrix and compute the e-ph matrix element with PERTURBO code [91] which is then used in Eq. (19b) to compute the screened self-energy. On the other hand, to compute the phonon anharmonicity, we first use the ab initio molecular dynamics method implemented in VASP to generate the interatomic force and then compute the anharmonic phonon frequency by the ALAMODE code [60].
We present the temperature-dependent phonon dispersion in Fig. 2(b), where we plot the imaginary frequency in the minus axis. The black line showing a robust soft mode at the M point is the direct result of DFPT, which is consistent with the formation of $2\times 2$ CDW of monolayer $\rm TiSe_{2}$. In the inset, we plot the corresponding normal mode for the soft phonon, which has atomic displacement perpendicular to the wave vector. In the hexagonal Brillouin zone (BZ), there are three inequivalent M points, and the condensation of these three normal modes can build up the lattice distortion observed in the CDW phase. Beyond the harmonic approximation, we can see that, by disregarding some temperature-independent correction, the soft phonon is the only mode sensitive to the anharmonic effect. Overall, the properties presented in the above results validate the single-mode assumption applied in Sect. III.
In Fig. 2(c), we analyze the soft-phonon at the M point. We present the phonon frequency as a function of temperature. We consider three different cases: that considers the correction only from the screening (Eq. (19b)), that considers the correction only from the anharmonic effect (Eq. (24)), and that considers both the two effects. First, we note that, at zero-temperature, the harmonic DFPT gives us a phonon frequency $\omega^{\rm DFPT}_{\bf M}=116i~{}{\rm cm^{-1}}$. When considering the temperature dependence in the Coulomb screened self-energy, the frequency varies by $30i~{}{\rm cm^{-1}}$ from 0 K to 600 K, which indicates that the phonon will always be soft in the harmonic approximation. On the other hand, the pure anharmonic effect can raise the phonon frequency dramatically, forcing the phonon mode to become ”hard” above $T=270$ K. This critical temperature can even be reduced to $T=190$ K when we combine the contribution from the two effects and solve the phonon frequency self-consistently. In Fig. 3, we present the self-energy change with respect to the zero-temperature value for phonon momentum $\mathbf{q}$ along the high symmetry line M-$\Gamma$ under different temperatures. We note that the most significant temperature dependence does not occur for the momentum $\mathbf{q}={\bf M}$ since the Fermi surface is not ideally ”nesting.” Instead, it is approximately a constant for $\mathbf{q}$ near the M point within a specific range and gradually decreases outside.

Based on the discussion in the last section, the EI/CDW lattice distortion is attributed to the combined effect of the three inequivalent M phonons. This so-called ”triple-$\mathbf{q}$” effect requires us to extend the ”single-$\mathbf{q}$” formalism presented in Sect. III. More explicitly, we write down the new wave function by introducing the pairing direction in the Eq. (26):

$$\phi^{n}_{(\mathbf{k})}=\sum_{v}t^{n}_{v\mathbf{k}}\psi_{v\mathbf{k}}+\sum_{c,I}t^{n}_{c\mathbf{k},I}\psi_{c\mathbf{k}+{\bf M}_{I}}$$

(48)

where $I=1,2,3$ is the index for the three M directions. In a simplified two band model, the mixing coefficients can be obtained by solving the triple-$\mathbf{q}$ version of the Hamiltonian Eq. (31a):

$$\begin{pmatrix}\xi^{c}_{\mathbf{k}+{\bf M}_{1}}&0&0&\Delta^{cv}_{\mathbf{k},1}\\ 0&\xi^{c}_{\mathbf{k}+{\bf M}_{2}}&0&\Delta^{cv}_{\mathbf{k},2}\\ 0&0&\xi^{c}_{\mathbf{k}+{\bf M}_{3}}&\Delta^{cv}_{\mathbf{k},3}\\ \Delta^{vc}_{\mathbf{k},1}&\Delta^{vc}_{\mathbf{k},2}&\Delta^{vc}_{\mathbf{k},3}&\xi^{v}_{\mathbf{k}}\end{pmatrix}\begin{pmatrix}t^{n}_{c\mathbf{k},1}\\ t^{n}_{c\mathbf{k},2}\\ t^{n}_{c\mathbf{k},3}\\ t^{n}_{v\mathbf{k}}\end{pmatrix}=E_{n(\mathbf{k})}\begin{pmatrix}t^{n}_{c\mathbf{k},1}\\ t^{n}_{c\mathbf{k},2}\\ t^{n}_{c\mathbf{k},3}\\ t^{n}_{v\mathbf{k}}\end{pmatrix},$$

(49)

where we neglect the pairing between states of $\mathbf{k}+{\bf M}_{I}$ and $\mathbf{k}+{\bf M}_{J}$ for which the momentum transfer exceeds the first BZ. It can be shown that the effect of triple-$\mathbf{q}$ will decouple from each other in the limit of $\Delta_{\mathbf{k}}\rightarrow 0$ and reduce to the single-$\mathbf{q}$ scenario in the linear order as Eq. (III.3). This property reflects that the inclusion of the triple-$\mathbf{q}$ effect only changes the finite $\Delta_{\mathbf{k}}$ behavior but maintains the critical temperature, $T_{c}$ where $\Delta_{\mathbf{k}}=0$. To correctly describe the behavior below $T_{c}$, in the implementation, we will construct and solve the gap equation with the triple-$\mathbf{q}$ effect included.
To comprehensively understand the role of each component in the gap equation, we first note that the temperature dependence in the EI/CDW gap equations lies in: electronic occupation, phonon frequency, electron-phonon coupling, and mass renormalization coefficients, and each of the factors is crucial to compute the critical temperature. To emphasize the importance, we first apply the MF gap equations, Eq. (37a) and Eq. (37b), as a baseline and use the DFT band structure with the phonon frequency near the experimental critical temperature, $T^{\rm exp.}_{c}=220~{}$K. However, the result gives an extreme high $T_{c}$, since the small frequency, $\Omega_{\bf M}(T_{c}^{\rm exp.})\approx 30~{}{\rm cm^{-1}}$, can boost the phonon thermal average with small potential $\Delta_{b}$ and enhance the coupling constant by Eq. (35). These two effects make the critical temperature higher than $10^{6}$ K.
To analyze how the phonon frequency affects the critical temperature, we restrict discussions in the MF gap equation and study three different scenarios. First, we restore the temperature dependence in the anharmonic phonon frequency without correction on Coulomb screened self-energy. A great improvement in the critical temperature is obtained (see blue line in Fig. 4), which gives $T_{c}\sim 3300$ K with the critical phonon frequency $\Omega^{c}_{\bf M}=308~{}{\rm cm^{-1}}$. Compared to the fixed-frequency case, a ten-times enhancement on $\Omega_{\bf M}$ can suppress the critical temperature by three orders of magnitude. This modification reflects that in a BCS-like mean-field theory, $T_{c}\sim\exp(1/U)$, where the pairing potential $U$ is characterized by $|g_{\bf M}^{2}|/\Omega_{\bf M}$ here [92, 93].
Next, we add the GW correction to increase the relative energy between occupied and unoccupied states. The result with a modified band structure differs from the calculation with DFT bands only at temperatures higher than 1200 K (see the orange line in Fig. 4); on the other hand, at temperatures lower than 1200 K, the coupling $\Delta_{\mathbf{k}}^{cv}$ is much larger than 1 eV, which makes the GW correction irrelevant. Last, fixed on the GW band, we use the phonon frequency dressed by the temperature-dependent Coulomb screening, which can suppress the lattice distortion potential in the whole temperature ranges and lower the $T_{c}$ to 1950 K. We found that although the $T_{c}$ is significantly reduced, the critical frequency $\Omega^{c}_{\bf M}$ is merely changed (see the table in Fig. 4). Thus, we conclude that in the mean-field approach, the major temperature-dependent factor in EI/CDW phase transition is dominated by the phonon frequency; disregarding the temperature, as long as the phonon has a frequency lower than $\sim 275~{}{\rm cm^{-1}}$, the lattice structure becomes unstable and transit into the distorted $2\times 2$ superlattice.