Pump driven normal-to-excitonic insulator transition: Josephson oscillations and signatures of BEC-BCS crossover in time-resolved ARPES

TL;DR

This paper models a pump-induced transition from a band insulator to a nonequilibrium excitonic insulator, revealing oscillations, spectral shape changes, and BEC-BCS crossover signatures observable via time-resolved ARPES and photocurrent oscillations.

Contribution

It provides the first ab initio band structure of LiF in NEQ-EI states and links spectral and photocurrent oscillations to Josephson-like phenomena in a pump-driven system.

Findings

Spectral convex-to-concave transition during pump driving.

Observation of BEC-BCS crossover signatures in ARPES spectra.

Photocurrent oscillations with pump-tunable frequency post-pump.

Abstract

We consider a ground-state wide-gap band insulator turning into a nonequilibrium excitonic insulator (NEQ-EI) upon visiting properly selected and physically relevant highly excited states. The NEQ-EI phase, characterized by self-sustained oscillations of the complex order parameter, neatly follows from a Nonequilibrium Green's Function treatment on the Konstantinov-Perel' contour. We present the first {\em ab initio} band structure of LiF, a ground-state bulk insulator, in different NEQ-EI states and show that these states can be generated by currently available pump pulses. We highlight two general features of time-resolved ARPES spectra: (1) during the pump-driving the excitonic spectral structure undergoes a convex-to-concave shape transition and {\em concomitantly} the state of the system goes through a BEC-BCS crossover; (2) attosecond pulses shone after the pump-driving at different times $t_{\rm delay}$ generate a photocurrent which {\em oscillates} in $t_{\rm delay}$ with a pump-tunable frequency -- we show that this phenomenon is similar to the AC response of an exotic Josephson junction.

Figures

Equations

$$\hat{H}$$

$$G^{\alpha\beta}_{k}(z,z^{\prime})=\frac{1}{i}\frac{{\rm Tr}\left[\mbox{$\mathcal{T}$}\left\{e^{-i\int d\bar{z}\hat{H}(\bar{z})}\hat{\psi}_{\alpha k}(z)\hat{\psi}_{\beta k}^{{\dagger}}(z^{\prime})\right\}\right]}{{\rm Tr}\left[\mbox{$\mathcal{T}$}\left\{e^{-i\int d\bar{z}\hat{H}(\bar{z})}\right\}\right]}$$

$${\hat{H}}^{\mathrm{M}}=\hat{H}-{\mu }_v{\hat{N}}_v-{\mu }_c{\hat{N}}_c,$$

$$\left[i\frac{d}{dz}-h_k(z)\right]G_k(z,z^{\prime })=\delta (z,z^{\prime })+\int d\bar{z}{\Sigma }_k(z,\bar{z})G_k(\bar{z},z^{\prime })$$

$$h_k(z)$$

$$\equiv$$

$$V_k^{\alpha \alpha }(z)$$

$$V_k^{\alpha \bar{\alpha }}(z)$$

$$h_k-m\equiv h_k^{\mathrm{M}},$$

$$G_k^{\mathrm{M}}(\tau ,{\tau }^{\prime })$$

$$G_k^{\alpha \beta ,\mathrm{M}}(\tau ,{\tau }^{+})=i\sum _{\lambda }f(e_k^{\lambda }){\phi }_{\alpha k}^{\lambda }{\phi }_{\beta k}^{\lambda },$$

$$V_k^{\alpha \alpha }$$

$$V_k^{\alpha \bar{\alpha }}$$

$${\widetilde{ϵ}}_{vk}$$

$${\widetilde{ϵ}}_{ck}$$

$${\Delta }_k$$

$$e_k^{\lambda }=\frac{{\widetilde{ϵ}}_{vk}+{\widetilde{ϵ}}_{ck}+\lambda R}{2},$$

$$\varphi_{k}^{\lambda}=\left(\begin{array}[]{c}\lambda\sqrt{\frac{1}{2}(1+\lambda\frac{\tilde{\epsilon}_{vk}-\tilde{\epsilon}_{ck}}{R})}\\ {\rm sign}(\Delta_{k})\sqrt{\frac{1}{2}(1-\lambda\frac{\tilde{\epsilon}_{vk}-\tilde{\epsilon}_{ck}}{R})}\end{array}\right),$$

$$G_k^{\alpha \beta ,<}(t,t^{\prime })=i\sum _{\lambda }f(e_k^{\lambda }){\phi }_{\alpha k}^{\lambda }(t){\phi }_{\beta k}^{\lambda \ast }(t^{\prime })$$

$$i\frac{d}{dt}\varphi_{k}^{\lambda}(t)=\big{[}h_{k}+V_{k}(t)\big{]}\;\varphi_{k}^{\lambda}(t),$$

$${\phi }_k^{\lambda }(t)=e^{-i(e_k^{\lambda }+\mu -{\sigma }_z\frac{\delta \mu }{2})t}{\phi }_k^{\lambda },$$

$$\mu \equiv \frac{{\mu }_c+{\mu }_v}{2}$$

$$\delta \mu \equiv {\mu }_c-{\mu }_v$$

$$\big{[}h_{k}+V_{k}(t)\big{]}\;\varphi_{k}^{\lambda}(t)=e^{-i(e^{\lambda}_{k}+\mu-\sigma_{z}\frac{\delta\mu}{2})t}\big{[}h_{k}+V_{k}\big{]}\varphi_{k}^{\lambda}.$$

$$\displaystyle\big{[}h_{k}+V_{k}\big{]}\varphi_{k}^{\lambda}$$

$$\big{[}h_{k}+V_{k}(t)\big{]}\;\varphi_{k}^{\lambda}(t)=\big{[}e_{k}^{\lambda}+\mu\mathbbm{1}-\frac{\delta\mu}{2}\,\sigma_{z}\big{]}\varphi_{k}^{\lambda}(t)$$

$${\Delta }_k(t)={\Delta }_k{e}^{-i\delta \mu t},$$

$${ϵ}_{\alpha k}=(-{)}^{\alpha }[2T(1-\cos (k))+{ϵ}_g/2]$$

$$E=-\frac{i}{2}\sum _k\int \frac{d\omega }{2\pi }\mathrm{Tr}\left[(\omega +h_k)G_k^{<}(\omega )\right]$$

$$A_k^{<}(\omega )=-i\mathrm{Tr}[G_k^{<}(\omega )]$$

$$A_k^{>}(\omega )=A_k^{<}(-\omega ).$$

$$|{\Psi }_{\mathrm{x}}⟩=\prod _p({\phi }_{vp}^{-}{\hat{v}}_p^{†}+{\phi }_{cp}^{-}{\hat{c}}_p^{†})|0⟩,$$

$${\hat{b}}_q^{†}=\sum _k{Y}_k^{(q)}{\hat{c}}_{k+q}^{†}{\hat{v}}_k$$

Full text

I Introduction

At low enough temperature small gap insulators, semimetals as well as metals with overlapping bands may turn into Excitonic Insulators (EI) due to the Coulomb electron-hole attraction Blatt et al. (1962); Keldysh and Kopaev (1965); Kozlov and Maksimov (1965); Jérome et al. (1967); Halperin and Rice (1968); Keldysh and Kozlov (1968); Comte and Nozières (1982). Theoretical and experimental works on different materials, e.g., quantum Hall bilayers Spielman et al. (2001); Eisenstein and MacDonald (2004), graphene bilayers Zhang and Joglekar (2008); Min et al. (2008); Lozovik and Sokolik (2008); Kharitonov and Efetov (2008), carbon nanotubes Varsano et al. (2017); Hellgren et al. (2018) and chalcogenide-based structures Cercellier et al. (2007); Monney et al. (2012); Zenker et al. (2013, 2014); Kogar et al. (2017), have addressed signatures and properties of this equilibrium EI phase, including dynamical responses to external laser pulses Golež et al. (2016); Murakami et al. (2017); Mor et al. (2017); Werdehausen et al. (2018a, b); Tanabe et al. (2018) and adiabatic switching of electronic correlations Tuovinen et al. (2019). One of the most interesting signatures of the EI phase is the flattening of the valence band, recently revealed in Ta2NiSe5 using Angle-Resolved Photoemission Spectroscopy (ARPES) Wakisaka et al. (2009); Seki et al. (2014).

The last decade has seen a renewed interest in systems which are semimetals Triola et al. (2017); Pertsova and Balatsky (2018) or insulators Szymańska et al. (2006); Hanai et al. (2016, 2017, 2018); Becker et al. (2019); Yamaguchi et al. (2012, 2013) in the ground state but exhibit an EI phase in some (possibly pump-induced) nonequilibrium (NEQ) excited state. The optical properties of a NEQ-EI have been calculated by several authors in the past Schmitt-Rink et al. (1988); Haug (1985); Glutsch and Zimmermann (1992); Chu and Chang (1996). However, it was not until 1993 that Östreichand and Schönhammer pointed out a crucial difference between a ground state EI and a NEQ-EI Östreich and Schönhammer (1993): in a NEQ-EI the macroscopic polarization after the pump has been switched off exhibits persistent, self-sustained oscillations with a finite, time-independent amplitude and a frequency which depends on the absorbed energy. A similar dynamical excited state of matter has been independently shown to emerge in mean-field approximations using an ansatz for the time-dependent order parameter Szymańska et al. (2006). Here we provide an alternative derivation based on nonequilibrium Green’s functions (NEGF) and show the equivalence between the pump-induced state Östreich and Schönhammer (1993) and the dynamical excited state Szymańska et al. (2006).

The stability of the NEQ-EI phase plays of course a crucial role in applications to, e.g., optoelectronics Jun et al. (2017). Hanai et al. have pointed out the existence of regimes for which small fluctuations of the order parameter destroy the phase (dynamical instability) Hanai et al. (2017) and also studied the effects of an exciton-boson coupling Hanai et al. (2016). Hannewald et al. studied the NEQ-EI life-time due to dephasing Hannewald et al. (2000). A beyond mean-field analysis including light-matter interaction and damping has confirmed the existence of the NEQ-EI phase Becker et al. (2019). Real-time first-principles studies are foreseeable in the near future for quantitative, material-dependent predictions.

In this work we apply the NEGF formalism to systems which exhibit a NEQ-EI phase but are normal Band Insulators (BI) in the ground state, and highlight the NEQ-EI fingerprints visible in ARPES spectra. We first show that pump-pulses with a properly chosen subgap frequency generate the same dynamical excited state Östreich and Schönhammer (1993) as an excited self-consistent calculation. The proposed self-consistent NEGF theory is then implemented in the ab-initio Yambo code Marini et al. (2009) and the spectral function of a LiF bulk insulator in the NEQ-EI phase is presented and discussed. Our main findings are: (1) for sufficiently long and non-resonant pump-pulses the system undergoes a BEC-BCS crossover Eagles (1969); Strinati et al. (2018) which causes a convex-to-concave shape transition in the excitonic structure of time-resolved ARPES spectra; (2) in the NEQ-EI phase the photocurrent generated by attosecond pulses shone at different times $t_{\rm delay}$ oscillates in $t_{\rm delay}$ with a pump-tunable frequency: we relate this phenomenon to the AC response of an exotic Josephson junction.

The paper is organized as follows. In Section II we illustrate the self-consistent NEGF approach to excited-state systems and derive the equations to calculate the Green’s function. In Section III we discuss the NEQ-EI phase diagram of a two-band model Hamiltonian. We determine the critical conduction density at the BEC-BCS crossover and the boundary between the BI and NEQ-EI phases. We further illustrate how the NEQ spectral function changes in these different regimes. The simplifications of the two-band model are relaxed in Section IV where we present the spectral function of a LiF bulk insulator in the NEQ-EI phase. In Section V we perform real-time simulations and obtain the optimal pump pulse which brings the system closest to the self-consistent NEGF results. The time-resolved ARPES spectrum is calculated in Section VI for different probe durations, highlighting the signatures of the BEC-BCS crossover and the Josephson-like oscillations. A summary of the main results and concluding remarks are contained in Section VII.

II Model Hamiltonian and excited-state Green’s function

To describe the main qualitative aspects of the NEQ-EI phase we initially discard several complications of realistic materials. We consider a one-dimensional spinless insulator with one valence band and one conduction band separated by a direct gap of magnitude $\epsilon_{g}$, and assume that the intra-band repulsion between the electrons is negligible. A similar model has been proposed in Ref. Yang et al., 2012 to discuss excitons within TDDFT and more recently in Refs. Golež et al. (2016); Mor et al. (2017); Murakami et al. (2017) to drive out of equilibrium systems which are EI in the ground state. In Section IV we apply the NEGF theory to LiF and relax all the simplications of the model; multiple bands and valleys, intra-band and inter-band repulsion, band anistropies and degeneracies as well as spin-exchange effects will be all taken into account.

The explicit form of the model Hamiltonian reads

[TABLE]

where $\hat{v}_{k}$ ($\hat{c}_{k}$) annihilates an electron of momentum $k$ in the valence (conduction) band, $U_{q}$ is the Coulomb interaction and $\mathcal{N}$ is the number of discretized $k$-points. The last term in the first row describes the interaction of (conduction) electrons with the positive background. The Hamiltonian in Eq. (1) commutes with the number operator of valence and conduction electrons; hence each many-body eigenstate is characterized by a well defined number of electrons $N_{v}$ and $N_{c}$ in these bands. A nonvanishing value of the average $\langle\hat{c}^{{\dagger}}_{k}\hat{v}_{k}\rangle$ is possible only provided that several many-body eigenstates with different $N_{c}$ and $N_{v}$, constrained by (N_{c}+N_{v})/\mbox{\mathcal{N}}=1, become degenerate in the thermodynamic limit. This can happen for the ground-state multiplet (in this case the exciton energy $\epsilon_{\rm x}\lesssim 0$) and/or for excited-state multiplets. In this work we are especially interested in the latter scenario.

We consider $\epsilon_{g}/U_{0}$ large enough for the exciton energy to satisfy $0\ll\epsilon_{\rm x}\lesssim\epsilon_{g}$. Then, the ground state is nondegenerate and it consists of a fully occupied valence band and a completely empty conduction band (consequently the ground-state average of the interaction with the positive background vanishes). To investigate the existence of a NEQ-EI phase we use the NEGF formalism. In NEGF the fundamental quantity is the Keldysh-Green’s function which for our model Hamiltonian reads Stefanucci and van Leeuwen (2013)

[TABLE]

where $\hat{\psi}_{\alpha k}=\hat{v}_{k},\hat{c}_{k}$ for $\alpha=v,c$. In Eq. (2) the arguments $z,z^{\prime}$ as well as the integral over $\bar{z}$ run on the Konstantinov-Perel’ contour Konstantinov and Perel’ (1961) consisting of a forward and backward branch $(0,\infty)\cup(\infty,0)$ joined to a vertical imaginary track $(0,-i\mbox{\textbeta})$ with \textbeta the inverse temperature; $\mathcal{T}$ is the contour ordering operator. In the absence of external fields $\hat{H}(z)=\hat{H}$ for $z$ on the forward or backward branches. If the system is initially in thermal equilibrium at chemical potential $\mu$ then the Hamiltonian on the imaginary track is $\hat{H}(z)\equiv\hat{H}^{\rm M}=\hat{H}-\mu\hat{N}$ with $\hat{N}=\hat{N}_{v}+\hat{N}_{c}$ and $\hat{N}_{\alpha}\equiv\sum_{k}\hat{\psi}^{{\dagger}}_{\alpha k}\hat{\psi}_{\alpha k}$ the number operator for electrons in band $\alpha$. This choice of $\hat{H}^{\rm M}$ does indeed correspond to averaging with the thermal equilibrium density matrix \hat{\rho}=e^{-\scriptsize{\mbox{\textbeta}}(\hat{H}-\mu\hat{N})}/\mbox{\mathcal{Z}}, $\mathcal{Z}$ being the partition function. The resulting $G^{\alpha\beta}_{k}(z,z^{\prime})$ in Eq. (2) is the Matsubara Green’s function for both $z$ and $z^{\prime}$ on the vertical track and the equilibrium real-time Green’s functions for $z$ and $z^{\prime}$ on the forward or backward branches Stefanucci and van Leeuwen (2013).

To calculate the Green’s function in some excited state we must average with an excited density matrix \hat{\rho}=e^{-\scriptsize{\mbox{\textbeta}}\hat{H}^{\rm M}}/\mbox{\mathcal{Z}}, where $\hat{H}^{\rm M}$ is a self-adjoint operator with the property $[\hat{H}^{\rm M},\hat{H}]=0$. Here we consider

[TABLE]

with $\mu_{v}\neq\mu_{c}$. This $\hat{H}^{\rm M}$ commutes with $\hat{H}$ since $[\hat{N}_{\alpha},\hat{H}]=0$. For $U_{q}=0$ and zero temperature averaging with $\hat{\rho}$ is equivalent to averaging over a state with the valence (conduction) band populated up to the chemical potential $\mu_{v}$ ($\mu_{c}$). Hereafter we use $\hat{H}^{\rm M}$ in Eq. (3) for the Hamiltonian on the imaginary track.

The exact equation of motion for the Green’s function in Eq. (2) reads (in a $2\times 2$ matrix form)

[TABLE]

where $\Sigma_{k}$ is the sum of the Hartree-Fock (HF) and correlation self-energies and $h_{k}(z)$ is given by

[TABLE]

with $\theta_{|}(z)=0$ if $z$ is on the forward or backward branches and $\theta_{|}(z)=1$ if $z$ is on the imaginary track. In the Hartree-Fock (HF) approximation we have $\Sigma_{k}^{\alpha\beta}(z,z^{\prime})=\delta(z,z^{\prime})V^{\alpha\beta}_{k}(z)$, where the HF potential reads

[TABLE]

with $\bar{\alpha}=c,v$ for $\alpha=v,c$, and $z^{+}$ the contour-time infinitesimally later than $z$. As we have discarded the intraband interaction the diagonal elements of $V$ are only due to the Hartree diagram whereas the off-diagonal ones are only due to the exchange (Fock) diagram. This simplification is relaxed in Section IV.

II.1 Matsubara Green’s function

Choosing $z=-i\tau$ and $z^{\prime}=-i\tau^{\prime}$ on the imaginary track we get the Matsubara Green’s function $G^{\rm M}_{k}(\tau,\tau^{\prime})\equiv G_{k}(-i\tau,-i\tau^{\prime})$. Since the HF self-energy is local in time the structure of $G^{\rm M}$ is the same as in the noninteracting case, the only difference being that the one-body noninteracting matrix

[TABLE]

see Eq. (10), is replaced by the one-body HF matrix $h_{k}^{\rm M}+V_{k}$, with $V_{k}\equiv V_{k}(-i\tau)$ the HF potential on the imaginary track Stefanucci and van Leeuwen (2013). Therefore

[TABLE]

where $f(\omega)=1/(e^{\scriptsize{\mbox{\textbeta}}\omega}+1)$ is the Fermi function and $\bar{f}(\omega)=1-f(\omega)$. As $V_{k}$ depends on $G_{k}^{\rm M}$, see Eqs. (11,12), the problem has to be solved self-consistently. We point out that for the Matsubara Green’s function to be correctly antiperiodic in $\tau$ and $\tau^{\prime}$ with period \textbeta, the population $n_{\alpha k}=-iG^{\alpha\alpha,\rm M}_{k}(\tau,\tau^{+})$ of band $\alpha=v,c$ is not an input of the self-consistency cycle.

Let $e_{k}^{\lambda}$ and \varphi_{k}^{\lambda}=\left(\begin{array}[]{c}\varphi_{vk}^{\lambda}\\ \varphi_{ck}^{\lambda}\end{array}\right), $\lambda=\pm$, be the two eigenvalues and eigenvectors of $h_{k}^{\rm M}+V_{k}$. Without any loss of generality we choose $\varphi_{\alpha k}^{\lambda}$ real. Then, from Eq. (14),

[TABLE]

and therefore the matrix elements of the HF potential in Eqs. (11-12) can be rewritten as

[TABLE]

In terms of the quantities

[TABLE]

it is straightforward to find the eigenvalues

[TABLE]

with $R=\sqrt{(\tilde{\epsilon}_{vk}-\tilde{\epsilon}_{ck})^{2}+4\Delta_{k}^{2}}$, and eigenvectors

[TABLE]

which inserted in Eqs. (16,17) provide a close system of equations for the HF potential. Henceforth we refer to $\Delta_{k}$ as the order parameter since a nonvanishing $\Delta_{k}$ implies that $\langle\hat{c}_{q}^{{\dagger}}\hat{v}_{q}\rangle\neq 0$ at least for some $q$, and hence a spontaneous symmetry breaking.

II.2 Real-time Green’s function

Once the Matsubara Green’s function is known we can calculate the Green’function with both arguments on the real axis. In the HF approximation the lesser Green’s function which solves Eq. (4) reads Stefanucci and van Leeuwen (2013)

[TABLE]

where the time-dependent vectors satisfy

[TABLE]

with boundary conditions $\varphi_{k}^{\lambda}(0)=\varphi_{k}^{\lambda}$ [we recall that $h_{k}(z)=h_{k}$ for $z$ on the forward/backward branch, see again Eq. (10)]. The time-dependent HF potential appearing in Eq. (24) can be calculated from Eqs. (11-12) by taking into account that for $z=t$ on the forward/backward branch we have $G(z,z^{+})=G^{<}(t,t)$.

Let us show that the solution of Eq. (24) is

[TABLE]

where $\sigma_{z}$ is the $2\times 2$ Pauli matrix,

[TABLE]

is the center-of-mass chemical potential and

[TABLE]

is the relative one. Inserting Eq. (25) in $G^{<}_{k}$ to evaluate the HF potential we find $V_{k}(t)=e^{i\sigma_{z}\frac{\delta\mu}{2}t}\,V_{k}\,e^{-i\sigma_{z}\frac{\delta\mu}{2}t}$. Taking into account that $h_{k}$ is diagonal we get

[TABLE]

Next we observe that $h_{k}=h_{k}^{\rm M}+m$, see Eq. (13), and that the $2\times 2$ matrix $m$ can be written as $m=\mu\mathbbm{1}-\frac{\delta\mu}{2}\,\sigma_{z}$. Hence

[TABLE]

where in the last row we used that $\varphi_{k}^{\lambda}$ are eigenvectors of $h_{k}^{\rm M}+V_{k}$. Thus, we can rewrite Eq. (28) as

[TABLE]

which coincides with the time-derivative $i\frac{d}{dt}\varphi_{k}^{\lambda}(t)$. We emphasize that all quantities needed to evaluate the lesser Green’s function have been previously calculated to determine the Matsubara Green’s function. The greater Green’s function can be derived similarly; the final result looks like Eq. (23) after replacing $if(e_{k}^{\lambda})\to-i\bar{f}(e_{k}^{\lambda})$.

Due to the presence of $\sigma_{z}$ in Eq. (25) only the diagonal elements of $G^{<}_{k}$ depend on the time-difference. The off-diagonal ones oscillate monochromatically at the frequency $\delta\mu/2$. In particular, the order parameter $\Delta_{k}(t)=V^{cv}_{k}(t)=[V^{cv}_{k}(t)]^{\ast}$ oscillates monochromatically at the frequency $\delta\mu$, i.e.,

[TABLE]

in agreement with Ref. Szymańska et al. (2006). The time-depencence of the order parameter in the NEQ-EI phase resembles the behavior of the superconducting order parameter in a Josephson junction. Josephson oscillations between two superconducting electrodes separated by a thin insulator arise by applying a DC voltage to the equilibrium junction, thus effectively introducing a difference in the electrochemical potentials. This mechanism has recently been investigated by replacing the superconductors with EI electrodes Rontani and Sham (2009); Joglekar et al. (2005); Hsu and Su (2015); Apinyan and Kopeć (2019). The NEQ-EI oscillations in Eq. (31) are of a slightly different nature. They can be understood in terms of an exotic Josephson junction where Cooper pairs are formed by electrons in different electrodes. In our system the electrodes are the valence and conduction band and the Cooper pairs are the bound excitons. As we shall see in Section VI.2, the oscillating behavior of $\Delta_{k}$ has interesting implications in time-resolved photocurrent spectra.

III Phase diagram

To illustrate the possible solutions of the self-consistent calculation we consider the band structure

[TABLE]

where $(-)^{v}\equiv-1$ and $(-)^{c}\equiv 1$, $T>0$, and a short-range Coulomb interaction $U_{q}=U$ independent of $q$. The wavevectors $k,q,\ldots$ vary in the first Brillouin zone $(-\pi,\pi)$. Then the HF potential $V^{\alpha\beta}_{k}$ is independent of $k$, see Eqs. (16,17), and hence $\Delta_{k}=\Delta$ is independent of $k$ too. We consider the system at zero temperature and express all energies in units of the bare gap $\epsilon_{g}$. We fix the bandwidths $W=4T=2$, the center-of-mass chemical potential $\mu=0$ and vary $U\geq 0$ and $\delta\mu\geq 0$. For $\delta\mu=0$ we recover the ground-state Green’s function. Figure 1 shows the color-plot of the self-consistent $\Delta$. In the white regions we have the trivial (and unique) solution $\Delta=0$. These are the BI regions characterized by $N_{c}=0$ or $N_{v}=0$, and hence by a gap $\epsilon_{g}$ between occupied an unoccupied single-particle states. In the colored region we have two solutions, the trivial one ($\Delta=0$) and the nontrivial one ($\Delta\neq 0$). For both solutions we have estimated the total energy according to Stefanucci and van Leeuwen (2013)

[TABLE]

and found that the energy of the nontrivial solution is always the lowest. In Section III.1 we calculate the spectral function and show that the unoccupied and occupied single-particle states are separated by an energy gap $\Delta$. Hence the system is still an insulator. However, since $N_{c}$ and $N_{v}$ are both nonvanishing this insulating phase can only be due to the Coulomb attraction between a conduction electron and a valence hole. We then say that the system is in the EI phase. The EI phase exists in the ground state ($\delta\mu=0$) for large enough $U$ as well as in excited states ($\delta\mu\neq 0$) even for rather small $U$’s. The region below the blue diamonds is characterized by a nontrivial solution with small $\Delta$ (in this region $\Delta\leq 10^{-2}$); here the system behaves essentially like a normal metal (NM).

III.1 Spectral function

From the NEGF soluton of Section II.2 we can extract the occupied part of the spectral function according to

[TABLE]

where $G^{\alpha\alpha,<}_{k}(\omega)$ is the Fourier transform of $G^{\alpha\alpha,<}_{k}(t,t^{\prime})$ with respect to the time-difference. In Fig. 2 we show results for $U=1$ and three different $\delta\mu=0.82,\,0.93,\,1.05$ corresponding to a density of electrons in the conduction band n_{c}=-i\sum_{k}G^{cc,<}_{k}(t,t)/\mbox{\mathcal{N}}=0.07,0.15,\,020. Notice that for $U=1$ the nontrivial solution $\Delta\neq 0$ exists only for values of $\delta\mu\gtrsim 0.76$, see again Fig. 1. Interestingly, for our parameters the energy of the zero-momentum exciton is $\epsilon_{\rm x}\simeq 0.76$. The fact that $\epsilon_{\rm x}$ coincides with the value of $\delta\mu$ at the BI-EI boundary is not a coincidence Yamaguchi et al. (2012). For completeness we provide a proof in Appendix A. For $\delta\mu=0.82$ (low densities $n_{c}$) an excitonic structure appears at roughly $\epsilon_{g}-\delta\mu$ below the conduction band minimum (CBM). This spectral function agrees with the results of Ref. Rustagi and Kemper (2018) where a system with one single exciton, i.e., n_{c}=1/\mbox{\mathcal{N}}, was considered.

By increasing further $\delta\mu$ (and hence the density in the conduction band) the excitonic structure changes its convexity, see middle and right panel in Fig. 2. This phenomenon is distinct from the one reported in Ref. Rustagi and Kemper (2018) where the change of convexity is obtained by averaging spectral functions each with one single exciton of different center-of-mass momentum. In our case the convex-to-concave shape transition develops with increasing the density of excitons, see also Section III.2.

The unoccupied part of the spectral function can be calculated similarly: $A^{>}_{k}(\omega)=-i{\rm Tr}[G^{>}_{k}(\omega)]$. Since $\mu=0$ the system is particle-hole symmetric and therefore

[TABLE]

As anticipated there is a gap of order $\Delta$ between the unoccupied and occupied bands.

III.2 BEC-BCS crossover

The EI phase is similar to the BCS superconducting phase when the width $\xi$ of the excitons is larger than the average electron distance (or Wigner-Seitz radius) $r_{s}$. In the opposite regime, i.e., $\xi\ll r_{s}$, the excitons behave like point-like bosons in a Bose-Einstein condensate (BEC). To determine the BEC-BCS crossover Comte and Nozières (1982); Nozières and Schmitt-Rink (1985); Yamaguchi et al. (2012); Strinati et al. (2018) we have calculated the width $\xi$ of the zero-momentum excitonic wavefunction and compare it with the average distance $r_{s}$ of electrons in the conduction band. The excitonic wavefunction in real space is given by (modulo a normalization constant) $Y(x)=\sum_{k}e^{ikx}Y_{k}$, with $Y_{k}$ the excitonic solution of Eq. (56), whereas (for one-dimensional systems) r_{s}\simeq\mbox{\mathcal{N}}/N_{c}. The green triangles in the phase diagram correspond to the values of $U$ and $\delta\mu$ for which $\xi=r_{s}$, where $\xi$ has been estimated from $Y(\xi)/Y(0)=0.1$.

An alternative way to find the values of $\delta\mu$ at the BEC-BCS crossover consists in determining the onset of the broadening of the excitonic population in momentum space – in the BEC regime the excitonic population is peaked around zero momentum. Let us write the many-body state of the system in the NEQ-EI phase. Since the minus branch $e_{k}^{-}$ is entirely below the plus branch $e_{k}^{+}$ we have

[TABLE]

where $|0\rangle$ is the state with no electrons. We define the creation operator for an exciton of momentum $q$ according to

[TABLE]

where the excitonic amplitude $Y_{k}^{(q)}$ can be calculated by solving the analogous of Eq. (56) for finite-momentum excitons. It is a matter of simple algebra to show that

[TABLE]

Therefore, the excitonic density in the NEQ-EI phase is given by

[TABLE]

In Fig. 3 we show $n_{q}^{\rm ex}$ for the interaction strength $U=1$ and different densities in the conduction band. The broadening starts for n_{c}=N_{c}/\mbox{\mathcal{N}}=1/r_{s}\simeq 0.15. This density is obtained when $\delta\mu\simeq 0.93$, in agreement with the abscissa of the green triangle in Fig. 1 for the same value of $U$.

Remarkably, the convex-to-concave shape transition in the excitonic structure of the spectral function roughly occurs at the BEC-BCS crossover. In Fig. 2 the Coulomb strength is $U=1$ and the transition occurs for $\delta\mu\simeq 0.93$, in agreement with the previously estimated value. We have verified (not shown) that this property remains true for all values of $U<U_{c}\simeq 2.3$ (for larger $U$’s the ground state is an EI). We infer that time-resolved ARPES spectra may provide a tool to detect the BEC-BCS crossover in NEQ-EI’s, see also Section VI.1.

IV LiF bulk insulator

We now apply the NEGF theory for excited states, see Section II, to a LiF bulk insulator. In the ground state this material has an experimental gap of about 14.5 eV and the optical spectrum exhibits a $1s$ exciton with a binding energy of $1.9$ eV Roessler and Walker (1967).

To describe the system we have chosen the Kohn-Sham (KS) basis of a self-consistent calculation performed with the QuantumEspresso package Giannozzi et al. (2009) at the level of the local density approximation (LDA). The ground-state has been converged using a $6\times 6\times 6$, $\Gamma$ centered grid with a cutoff of $80$ Ry. Then the wave-functions have been computed with a non self-consistent calculation on a $8\times 8\times 8$ grid. Our equilibrium bands well reproduce those of Ref. Wang et al. (2003) with a KS gap of $8.45$ eV, which is quite far from the experimental value. To improve the description we have constructed the Hartree plus Screened EXchange (HSEX) potential using the Yambo code Marini et al. (2009). Multiple bands and valleys, intra-band and inter-band repulsion, band anistropies and degeneracies as well as spin-exchange effects are all built in. A first approximated screened interaction has been evaluated in RPA with KS energies and wave-functions for the polarization. The HSEX spectrum is found to be well approximated by the original KS spectrum if the KS valence (conduction) energies are stretched by a factor 1.65 (1.4). The stretched KS energies with a scissor of 6 eV (necessary to reach the experimental gap) have been used to improve the polarization and hence the screened interaction. We verified that a further iteration of the screened interaction did not change substantially the stretching factors. The HSEX with improved interaction widens the gap from 8.5 eV to 12.8 eV; we have therefore applied a scissor of 1.7 eV to match the experimental value of the gap. The improved screened interaction has also been implemented in the Bethe–Salpeter equation with kernel given by the functional derivative of the HSEX potential. The optical spectrum well reproduces the experimental results. In particular the $1s$ excitonic peak is located at about $\epsilon_{\rm x}=12.48$ eV with respect to the VBM, hence its binding energy is $2$ eV.

Having a satisfactory description of LiF in equilibrium we have carried out the self-consistent Matsubara procedure of Section II using the HSEX potential with ground-state RPA screened interaction. For both equilibrium and excited-state calculations we have explicitely considered 50 bands. The screened interaction and the screened contribution to the exchange term has been calculated with $8$ Ry ($=59$ RL vectors) whereas the Hartree term and the unscreened contribution to the exchange term has been calculated with $64$ Ry ($965$ RL vectors). In the self-consistency cycles we have fixed the total charge per unit cell (charge neutrality) instead of the center-of-mass chemical potential $\mu$. Convergence paramenters for the HSEX calculation are chosen to have an error $<0.1~{}$eV.

In Fig. 4 we show the order parameter $\Delta_{0}$ and the macroscopic polarization $P$ versus $\delta\mu$. In LiF we have defined the order parameter as the following zero-momentum average

[TABLE]

where $h^{\rm HSEX}$ is the HSEX Hamiltonian, and the sum extends over all valence ($v$) and conduction ($c$) bands. Notice that Eq. (40) reduces to Eq. (20) in the two-band model if the HSEX hamiltonian is replaced by the HF one. By construction the order parameter $\Delta_{k}$ is basis-dependent. However, the transition to the NEQ-EI phase can also be detected by looking at the macroscopic polarization $P$, which is a measurable physical quantity and, as shown in Fig. 4, a good order parameter ($\Delta$ and $P$ are proportional in the two-band model system). For $\delta\mu<\epsilon_{\rm x}$ both $\Delta$ and $P$ vanish, in agreement with the findings of the previous Section.

LiF shares other important qualitative features with the model Hamiltonian. In Fig. 5 we display the occupied part of the spectral function for three different values of $n_{c}$. Both the development of the excitonic structure and the convex-to-concave shape transition with increasing $n_{c}$ are clearly visible. We also observe that excitons distribute among different bands giving rise to multiple excitonic peaks at fixed parallel momentum $k$.

The results presented in this Section support the use of the two-band model for understanding other general features of the NEQ-EI phase.

V Pump driven BI-EI transition

In this Section we demonstrate that the self-consistent NEQ-EI phase is accessible by shining suitable laser pulses on the BI ground state.

We consider again the Hamiltonian of Eq. (1) at zero temperature and set $\delta\mu=0$ (ground-state) and $U<U_{c}$ (BI phase with $\Delta=0$). The system is driven out of equilibrium by a time-dependent electric field $E(t)$ coupled to the valence-conduction dipole moments $D_{k}$. The driving Hamiltonian reads

[TABLE]

This light-matter coupling has been used to study changes in the optical properties due to a strong monochromatic electric field (dynamical Stark effects) Comte and Mahler (1986, 1988) and more recently weak monochromatic fields as those generated in quantum optical cavities Latini et al. . Here we are interested in the electronic properties of the system when the electric field is a pump-pulse centered around frequency $\omega_{P}$ with finite duration $T_{P}$:

[TABLE]

All results in this Section refer to a system with interaction strength $U=1$, dispersion of Eq. (32) with bandwidth $W=4T=2$, center-of-mass chemical potential $\mu=0$ and, for simplicity, dipole moments $D_{k}=D$ independent of $k$ (energies are in units of the bare gap $\epsilon_{g}$). The real-time simulations have been performed with the CHEERS code Perfetto and Stefanucci (2018). Since the driving Hamiltonian depends on $E_{P}$ and $D$ only through their product we define the Rabi frequency $\Omega_{P}\equiv E_{P}D$.

We evolve the system in the time-dependent HF approximation and express times in units of $1/\epsilon_{g}$. The main result of the simulations is that for $\omega_{P}$ larger than the exciton energy $\epsilon_{\rm x}$, see Eq. (56), and smaller than $\epsilon_{g}$ the modulus of the order parameter attains a constant and finite value after the pump, in agreement with Ref. Östreich and Schönhammer (1993). In Fig. 6(a) we show the steady values of the density in the conduction band $n_{c}$ and of the modulus of the order parameter $|\Delta|$ versus the pump frequency $\omega_{P}$ (we recall that for the chosen parameters $\epsilon_{\rm x}=0.76$) for a pump-pulse of duration $T_{P}=100$ and Rabi frequency $\Omega_{P}=0.05$, see Fig. 6(d). Independently of the values of $T_{P}$ and $\Omega_{P}$ the trend of Fig. 6(a) is general: there exists an optimal frequency $\omega_{\rm opt}\in(\epsilon_{\rm x},\epsilon_{g})$ (of course depending on $T_{P}$ and $\Omega_{P}$) at which the steady values of $n_{c}$ and $|\Delta|$ are simultaneously maximized. We point out that in the noninteracting case, $U=0$, the steady value of $|\Delta|$ is zero for all $\omega_{P}$, $T_{P}$ and $\Omega_{P}$ due to perfect dephasing.

In Fig. 6(b-c) we show $n_{c}(t)$ and ${\rm Re}[\Delta(t)]$ during the action of the pump-pulse with $\omega_{P}=\omega_{\rm opt}\simeq 0.87$. We observe the occurrence of persistent monochromatic oscillations in ${\rm Re}[\Delta(t)]$ after $T_{P}=100$. We therefore conclude that

[TABLE]

where $\omega_{E}$ is a pump-dependent frequency. The time-dependent behavior of the order parameter is the same as that of the self-consistent solution, see Eq. (31). We have extracted the frequency $\omega_{E}$ from the real-time solution ${\rm Re}[\Delta(t)]$, set $\delta\mu=\omega_{E}$ in the self-consistent calculation of Section II.1 and found that the self-consistent values of $n_{c}$ and $\Delta$ are identical to the steady-state values of $n_{c}$ and $|\Delta|$. This remains true for several different pump intensities (not shown). We therefore conclude that the pump-driven state is precisely the NEQ-EI state analyzed in detail in Section II.

A second evidence in favor of our conclusion is presented in Fig. 7 where we show the values of $\omega_{E}$ (empty circle) resulting from pumping at the optimal frequency (black triangle) for $T_{P}=100$ and different Rabi frequencies $\Omega_{P}$. In the limit of vanishing $\Omega_{P}$ (hence vanishing intensities of the electric field) the density pumped in the conduction band approaches zero and both $\omega_{E}$ and $\omega_{\rm opt}$ approach the energy $\epsilon_{\rm x}=0.76$ of the zero-momentum exciton (dashed green line). This is the same behavior of $\delta\mu$ versus $n_{c}$. Indeed, a vanishingly small value of $n_{c}$ in the self-consistent solution implies a vanishingly small value of $\Delta$ which we know to occur for $\delta\mu=\epsilon_{\rm x}$, see Appendix A.

VI Time-resolved ARPES

The time-resolved ARPES signal is proportional to the number of electrons $N_{k}(\omega)$ with energy $\omega$ and parallel momentum $k$ ejected by a probing pulse ${\mathbf{e}}(t)$. For arbitrary probe pulses we have Perfetto et al. (2016); Freericks et al. (2009)

[TABLE]

where the ionization self-energy reads

[TABLE]

and ${\mathbf{D}}_{k,\omega}^{\alpha}$ is the dipole matrix element between a band state $\alpha k$ and a continuum LEED state of energy $\omega$ and parallel momentum $k$.

VI.1 BEC-BCS crossover with femtosecond probe-pulses

From Eq. (44) we see that $N_{k}(\omega)$ is a complicated two-times convolution of $G^{<}(t,t^{\prime})$. Let us introduce the relative time $\tau\equiv t-t^{\prime}$ and the center-of-mass time $t_{\rm CM}\equiv(t+t^{\prime})/2$. Since the lesser Green’s function yields the probability amplitude for a hole to propagate freely from $t^{\prime}$ to $t$ the $G^{<}(t,t^{\prime})$ varies in $\tau$ no slower than the inverse of the energy of the highest occupied state. For a system in equilibrium this energy is the ionization potential; in our nonequilibrium system this energy can be estimated as the CBM, i.e., $\epsilon_{c0}$. The time-scale of the variation of $G^{<}$ in the center-of-mass time $t_{\rm CM}$ can be inferred from the rate of change of the occupations, see Fig. 6(b), and can be estimated as a fraction of the pump duration $T_{P}$. For band gaps in the eV range our pump pulse has a duration $T_{P}\simeq 10^{2}$ fs. Therefore, a probe pulse of duration $T_{p}\simeq T_{P}/10\simeq 10^{1}$ fs and central frequency $\omega_{p}\gg 2\pi/T_{p}$ is enough to resolve all removal energies provided that $T_{p}\gg 2\pi/\epsilon_{c0}$. If such a probe impinges on the system at time $t_{\rm delay}$ then $N_{k}(\omega)\propto A^{<}_{k}(t_{\rm delay},\omega-\omega_{p})$ where the transient spectral function is defined according to

[TABLE]

and $t_{\rm delay}=t_{\rm CM}$. Of course, for $t_{\rm delay}\gg T_{P}$ the transient spectral function becomes independent of $t_{\rm delay}$.

To obtain $A^{<}(t_{\rm delay},\omega)$ we observe that the time-dependent HF equations give access to the “time diagonal” $G_{k}^{<}(t,t)$ from which we can calculate the retarded HF Green’s function

[TABLE]

where $\sigma_{x}$ is the Pauli matrix, see Eq. (41). From $G^{\rm R}_{k}(t,t^{\prime})=[G^{\rm A}_{k}(t^{\prime},t)]^{\ast}$ we can extract the “time off-diagonal” $G_{k}^{<}(t,t^{\prime})$ as Stefanucci and van Leeuwen (2013)

[TABLE]

and calculate the transient spectral function from Eq. (46).

We consider again the pump-pulse in Fig. 6(d). The green dots are drawn at times $t=50,\,70,\,100$ for which $n_{c}=0.07\,,0.15,\,0.20$ respectively; these values of the conduction density are the same as in the three panels of Fig. 2. In Fig. 8 we show the color plot of $A^{<}(t_{\rm delay},\omega)$ at $t_{\rm delay}=50,70,100$. For $t_{\rm delay}\lesssim T_{P}$ the pump is still active and the transient spectral function is affected by nonadiabatic effects Myöhänen et al. (2010). However, these are very small and the agreement between Fig. 2 and Fig. 8 is unexpectedly good. Incidentally, we observe that this agreement supports the adiabatic approximation employed in Ref. Perfetto et al. (2016) to calculate time-resolved ARPES spectra. For $t_{\rm delay}\gg T_{P}$ the transient spectral function becomes independent of $t_{\rm delay}$ (as it should be) and almost indistinguishable from the right panel of Fig. 2 (not shown). This corroborates further the equivalence between the pump-driven state and the NEQ-EI state of Section II.

From Fig. 8 we also see that the convex-to-concave shape transition of the excitonic structure (signaling the BEC-BCS crossover) can be revealed by time-resolved ARPES spectra. It is worth emphasizing that linear response theory, i.e., $n_{c}\propto\Omega_{P}^{2}$, is totally inadequate to describe the pump-driven evolution of the system toward the BCS domain. In Fig. 9 we show the time-dependent evolution of $n_{c}(t)$ for $\omega_{P}=0.87$, $T_{P}=100$ and different values of the Rabi frequency $\Omega_{P}$. For $\Omega_{P}=0.02,\,0.04,\,0.05$ we display the corresponding steady-state ($t_{\rm delay}\gg T_{P}$) spectral functions. The concavity changes for $\Omega_{P}=0.05$ and for this Rabi frequency the steady density of conduction electrons is clearly not proportional to $\Omega_{P}^{2}$.

VI.2 Josephson oscillations with attosecond probe-pulses

The ARPES signal changes dramatically for ultrashort probe pulses since the probe has no time to wash out the oscillatory contribution of the off-diagonal $G^{<}$. To highlight the main qualitative difference we consider a $\delta$-like pulse of the form

[TABLE]

Then, the number of electrons in Eq. (44) becomes

[TABLE]

where we have defined the Rabi frequencies

[TABLE]

At the steady state $G^{\alpha\alpha,<}_{k}(t,t)$ is independent of $t$ but $G^{cv,<}_{k}(t,t)=-[G^{vc,<}_{k}(t,t)]^{\ast}\sim e^{-i\delta\mu t}$, see Eq. (23). Consequently $N_{k}(\omega)$ consists of a DC part $N^{\rm(DC)}_{k}(\omega)$ and an AC part of amplitude $N^{\rm(AC)}_{k}(\omega)$ and frequency $\delta\mu$:

[TABLE]

We therefore predict oscillations in $N_{k}(\omega)$ versus $t_{\rm delay}$ even for $t_{\rm delay}\gg T_{P}$. These oscillations should be observed in time-resolved ARPES provided that the duration of the probe pulse is much smaller than the period $2\pi/\delta\mu$. For values of $\delta\mu$ in the eV range (like in LiF) probe durations of the order of a hundred of attoseconds are sufficient.

The AC response in Eq. (52) can again be explained using the analogy with the exotic Josephson junction, see Section II.2. The application of a pump pulse coupling electrons in different electrodes (bands in our case) is equivalent to the application of a DC bias between the electrodes. If the DC bias is kept on for a time $T_{P}$ then a macroscopic number of electrons is transferred from one electrode (valence band) to the other (conduction band) thus generating a difference in the electrochemical potentials of the electrodes. Once the DC bias (pump in our case) is switched off this difference does not damp to zero since electrons cannot hop back to the original electrode (band). By all means this difference in electrochemical potentials is equivalent to the difference $\delta\mu$ discussed in Section II.

VII Summary and outlook

We presented a derivation of the NEQ-EI phase using the NEGF formalism on the Konstantinov-Perel’ contour. The NEQ-EI phase can be physically induced by pump-pulses with properly chosen subgap frequency, and the nature (BEC or BCS) of the excitonic condensate can be tuned by the pump duration. The most remarkable difference between a ground-state EI and a NEQ-EI is the time-dependent behavior of the order parameter and the total polarization. In fact, they both exhibit persistent, self-sustained monochromatic oscillations even at vanishing pump. Interestingly, these oscillations have the same nature of the AC oscillations in an exotic Josephson junction, where Cooper pairs are formed by electrons at the opposite side of the junction.

The NEQ-EI phase leaves clear fingerprints in time-resolved ARPES spectra. For probe durations long-enough to resolve the band structure the oscillatory contribution of the order parameter is washed out and an excitonic “band” appears inside the gap. Depending on the BEC or BCS nature of the excitonic condensate this band can be either convex or concave. Time-resolved ARPES experiments could in principle monitor the convex-to-concave shape transition using, e.g., pump-pulses of hundreds of femtoseconds and overlapping probe-pulses of few femtoseconds shone at different delays. Ultrafast probe pulses of, e.g., a few hundreds of attoseconds, do instead broaden the band structure and unveil the time-dependent contribution. The resulting photocurrent oscillates in time even after the end of the pump. The period of the oscillations is the same as the period of the order parameter, which is no smaller than the exciton energy $\epsilon_{\rm x}$.

The equivalence between the time-dependent approach and the self-consistent NEGF scheme for excited states allows for studying the NEQ-EI phase in realistic materials using different strategies. The inclusion of electronic correlations at, e.g., the GW level is more feasible in the self-consistent scheme. On the other hand, the effects of phonon-induced coherence-losses and exciton recombination on the self-sustained oscillations of the total polarization could be investigated using the real-time Kadanoff-Baym equations Stefanucci and van Leeuwen (2013) in the GKBA framework Lipavský et al. (1986). Of course accurate calculations are needed for quantitative predictions on specific materials. However, the highlighted qualitative features of the NEQ-EI phase are general and provide a useful guide for the interpretion of time-resolved ARPES spectra.

Appendix A Bethe-Salpeter equation and the BI-EI phase boundary

The equation for the BI-EI boundary marked with red crosses in Fig. 1 can be obtained analytically. Infinitesimally close to the boundary $\Delta_{k}$ is small. Then the minus band $e^{-}_{k}$ is entirely below the plus band $e^{+}_{k}$ and only $\lambda=-$ contributes in the sum of Eq. (17). Expanding Eq. (22) to lowest order in $\Delta_{k}$ and taking into account that $\tilde{\epsilon}_{vk}-\tilde{\epsilon}_{ck}<0$ we find

[TABLE]

Using this expansion the diagonal elements of the HF potential can be approximated to first order in $\Delta_{k}$ as $V^{vv}\simeq 0$ and $V^{cc}_{k}\simeq U_{0}$, see Eqs. (16). Consequently, from Eqs. (18,19) we get $\tilde{\epsilon}_{vk}\simeq\epsilon_{vk}-\mu_{v}$ and $\tilde{\epsilon}_{ck}\simeq\epsilon_{ck}-\mu_{c}$, which implies that $\tilde{\epsilon}_{vk}-\tilde{\epsilon}_{ck}=\epsilon_{vk}-\epsilon_{ck}+\delta\mu$. Using this result in Eq. (53), inserting the resulting expressions into Eq. (17) and recalling that, by the definition in Eq. (20), $V^{vc}_{k}=\Delta_{k}$ we get

[TABLE]

The lowest value of $\delta\mu$ for which this equation admits nontrivial solutions defines the BI-EI boundary.

Interestingly, Eq. (54) coincides with the eigenvalue equation of a many-body eigenstate with a single exciton of vanishing momentum. Let $|\Psi_{0}\rangle=\prod_{k}\hat{v}^{{\dagger}}_{k}|0\rangle$ be the ground state of energy $E_{0}$ in the BI phase and let us write the one-exciton state as

[TABLE]

where we introduced the eh states $|\Phi_{k}\rangle\equiv\hat{c}^{{\dagger}}_{k}\hat{v}_{k}|\Psi_{0}\rangle$. Since $\hat{H}$ preserves the number of electrons in each band $\hat{H}|\Psi\rangle$ is again a linear combination of the $|\Phi_{k}\rangle$’s. The possible excited state energies $E=E_{0}+\epsilon_{\rm x}$ are found by solving the eigenvalue problem $\hat{H}|\Psi\rangle=(E_{0}+\epsilon_{\rm x})|\Psi\rangle$. Inserting the expansion in Eq. (55) and taking the sandwich with $\langle\Phi_{k}|$ we get the Bethe-Salpeter equation for the coefficients $Y_{k}$

[TABLE]

Renaming $\epsilon_{\rm x}=\delta\mu$ and setting $Y_{q}=\Delta_{q}/(\epsilon_{vq}-\epsilon_{cq}+\delta\mu)$ we see by inspection that Eq. (56) is identical to Eq. (54). Thus, for a given $U$ the BI-EI boundary is found at $\delta\mu=\epsilon_{\rm x}$, in agreement with Ref. Yamaguchi et al. (2012). For $U=0$ the BI-EI boundary is at $\delta\mu=\epsilon_{g}$, as it should, whereas for $U\geq U_{c}\simeq 2.3$ (in this case the ground state is an EI) the solutions of Eq. (54) occur for $\delta\mu\leq 0$. The $\delta\mu\leq 0$ criterion is indeed used to establish the occurrence of an equilibrium EI phase Halperin and Rice (1968).