Electrical Breakdown of Excitonic Insulator

In this paper, we propose a new electrical breakdown mechanism for exciton insulators in the BCS limit, which differs fundamentally from the Zener breakdown mechanism observed in traditional band insulators. Our new mechanism results from the instability of the many-body ground state for exciton condensation, caused by the strong competition between the polarization and condensation energies in the presence of an electric field. We refer to this mechanism as ``many-body breakdown''. To investigate this new mechanism, we propose a BCS-type trial wave function under finite electric fields and use it to study the many-body breakdown numerically. Our results reveal two different types of electric breakdown behavior. If the system size is larger than a critical value, the Zener tunneling process is first turned on when an electrical field is applied, but the excitonic gap remains until the field strength reaches the critical value of the many-body breakdown, after which the excitonic gap disappears and the system becomes a highly conductive metallic state. However, if the system size is much smaller than the critical value, the intermediate tunneling phase disappears since the many-body breakdown happens before the onset of Zener tunneling. The sudden disappearance of the local gap leads to an ``off-on'' feature in the current-voltage ($I-V$) curve, providing a straightforward way to distinguish excitonic insulators from normal insulators.


I. INTRODUCTION
The excitonic insulator is an insulating phase where electron-hole pairs condensate [1][2][3][4].For semiconductors with a direct gap, this phase appears when the binding energy of excitons exceeds the band gap, as the spontaneous generation and condensation of excitons leads to a lowering of the system's energy.On the other hand, in semimetals, where the Fermi surfaces formed by electrons and holes match each other, the attractive Coulomb interaction binds the free electrons and holes at the Fermi level leading to BCS-type paring instability of the Fermi surfaces (FS).Such a BCS-like condensation in momentum space opens an energy gap at the FS, resulting in an excitonic insulating state.
The short lifetime of excitons in semiconductors often impedes the realization of exciton condensation.To overcome this challenge, quantum wells or van der Waals heterostructures are often used [5][6][7] where electrons and holes are separated in different layers.The real space separation of electrons and holes in these systems results in a longer lifetime for excitons, and the weak screening in 2D makes it easier for the binding to occur.So far, a large body of experimental [8][9][10][11][12][13] and theoretical [14][15][16][17][18] studies have been conducted to investigate the properties of excitonic insulators in such 2D bilayer systems.
Although excitonic insulators have been discussed in the literature for over half a century, very few material systems have been confirmed experimentally to exhibit such exotic states.This is because the exciton condensation only breaks the particle-hole U (1) symmetry, resulting in charge-neutral superfluidity, which is very hard to detect directly.In this study, we propose that the excitonic insulator in BCS limit may possess a unique breakdown mechanism, which can serve as a critical "smoking gun" type of experimental evidence, helping to distinguish an excitonic insulator from ordinary narrow-gap semiconductors.
We treat the breakdown problem of excitonic insulators by a simplified theoretical model, which contains a 2D bilayer system with a non-zero inter-layer distance, as shown in Fig. 1(a).The application of a vertical displacement field E ⊥ will result in the charging of the two layers by electrons and holes.If the interaction is absent, the charged bilayer would be expected to exhibit metallic behavior.However, the presence of an attractive interaction U (r) between electrons and holes will drive the system into an excitonic insulator state at the charge neutral point (CNP).When an in-plane electric field E is gradually added, the excitonic insulator is expected to be polarized and eventually broken down.
The most well-known intrinsic breaking down mechanism for band insulators is attributed to inter-band Zener tunneling [19][20][21][22][23].In an infinite system, the total energy becomes unbounded below when a uniform electric field is applied, resulting in the absence of a ground state.However, a finite system can still maintain an insulating stationary state at low electric fields [24,25].If we take the rigid band assumption and only include the electric field by a positional dependent chemical potential, the single particle Zener tunneling process can occur when the in-plane bias voltage eEL becomes comparable to the band gap ∆ as shown in Fig. 1(b).This means the Zener critical field is inversely proportional to the system size L. To go beyond the rigid band picture, Souza et al. [25] consider the polarization of the occupied bands and they find the 1/L behavior of the Zener field still stands.
We would emphasize that this field denotes the onset of Zener tunneling when a current proportional to the tunneling probability starts to flow.Under WKB approximation, the tunneling probability could be expressed as e − /ξ [26], where ξ is the correlation length determined by the gap ∆ and the tunneling length = ∆/eE is the width of the classically forbidden region for the Zener tunneling process.For an excitonic insulator in the BCS limit, ξ ≈ 4v F /π∆ is just the coherence length of the exciton condensate (details could be found in Appendix F).At fixed voltage, the current is exponentially small as the system size increases [27,28] due to its tunneling nature, and this current-carrying state is indeed a longlived resonance state as discussed by Souza et al. [29].
Setup of the bilayer system.A vertical displacement field E ⊥ is applied to equally charge the two layers with electrons and holes.At the CNP, attractive interaction U (r) between electrons and holes will drive the system into an excitonic insulator phase.When an in-plane electric field E is applied, the insulating system will be polarized and even broken down.(b) Under rigid band assumption, inter-band Zener tunneling happens only when the in-plane bias voltage exceeds the band gap, i.e. eEL > ∆.For any energy-allowed tunneling process, there exists a classically forbidden region (from B to A) with width = ∆/eE where the wavefunction decays.The correlation length of the gap ξ characterizes the penetration depth of the wavefunction into the classically forbidden region.(c) The rigid band assumption is valid only when the polarization effect is negligible.This is true for the BEC limit where the coherence length ξ is much smaller than the distance between exciton (roughly 1/ √ nex in 2D) and the inter-band Zener tunneling could also be viewed as single particle tunneling of each electron-hole pair.However, in the BCS limit, a macroscopic polarization P will be induced by an external electrical field and the competition between the polarization and the condensation energies will lead to the instability of the many-body ground state, resulting in the many-body breakdown.
In the BEC limit, the Zener tunneling current experiences a sharp increase when the electric field reaches ∆/eξ.This phenomenon can be explained as the single particle tunneling of each electron-hole pair in the center of mass system from a bound state to an extended state, similar to what happens in a normal insulator.However, in the BCS limit, the insulating behavior in an excitonic insulator is caused by the BCS type pairing wave function, which spontaneously breaks the particlehole U (1) symmetry and gives rise to a breakdown mechanism that is specific to BCS type excitonic insulators.This mechanism arises from the instability of the manybody ground state caused by the competition between polarization and condensation energies.
In the present paper, we will show that the critical field strength for this many-body breakdown is much smaller than that of the Zener breakdown (estimated roughly by ∆/eξ) in the BCS limit.Therefore, this unique electric breakdown feature can be considered as an important experimental signal for excitonic insulators, serving as a "smoking gun" to identify their presence.

II. POLARIZED MEAN FIELD THEORY
The actual breakdown scenario in excitonic insulators is complex since these two mechanisms could take effect at the same time.To better understand the breakdown of excitonic insulators, we will utilize a self-consistent mean field theory to analyze the interplay between Zener tunneling and the many-body breakdown.
For simplicity, we will limit our analysis to the lowest conduction band in the electron layer and the highest valence band in the hole layer which are described by the electron creation operators c † ek and c † hk , respectively.Within the k • p approximation, the many-body Hamiltonian can be written as where V is the area of the 2D system and Ψ(r) is field operator defined as Ψ(r) = V −1/2 sk e ik•r u s (r)c sk .The single particle Hamiltonian h 0 k is taken as where m e,h are effective masses of these two kind of electrons and E g > 0 is the original band gap.The exciton chemical potential µ ex is tuned by the vertical displacement field E ⊥ .Vanishing of the off-diagonal term in Eq.
(2) means direct inter-layer hopping is forbidden.This assumption is made because we are concerned with the breakdown of an excitonic gap.In real materials, this could be realized by symmetry constraints or just a large separation between layers.The inter-and intra-layer interaction are taken as the Coulomb ones: V (r) ≡ V s=s = e 2 / r and U (r) ≡ V s =s = e 2 / √ r 2 + d 2 whose Fourier transformations are V (q) = 2πe 2 / q, U (q) = V (q)e −qd .
Although an in-plane field breaks translation symmetry, to describe an insulating ground state, we can always take a trial state that keeps translation symmetry as long as the field is adiabatically added (the proof is in Appendix A).A trial HF state with translation symmetry at the CNP is |G = k c † vk |vac., where the valence band c † vk = α k c † ek + β k c † hk is a linear combination of the electron and hole band with constraints |α| 2 + |β| 2 = 1.By using Dirac notation |vk = [α k , β k ] T , energy per area becomes a functional of |vk , i.e.
where ρ ≡ ρ − ρ 0 is the density matrix relative to the initial uncharged state ρ 0 ss = δ ss δ sh and ρ is calculated as ρ ss k ≡ G|c † s k c sk |G = (|vk vk|) ss .A general form of this functional could be found in Appendix B.
The four terms in Eq. ( 3) could be viewed as kinetic, polarization, Hartree, and Fock energies separately.The Hartree energy is just the charging energy of the twolayer capacitor with the electron/hole density (exciton density) n ex = 1/V k ρ eek .The relative density matrix ρ is used in the Fock energy expression to avoid the double counting problem [14].For numerical convenience, a periodic boundary condition is assumed, and the polarization energy is calculated with the help of the expectation value of many-body position operators defined on a ring geometry [30], which is just a discrete form of Berry phase [31][32][33].This form of polarization energy functional has already been used to calculate the electrical properties of insulators in the literature [25,34,35].On the other hand, for the open boundary problem, the polarization energy functional should be written in real space by Wannier functions [36,37].
The local minimum is found by requiring the first order derivative of ε tot to be zero, i.e. δε tot /δ vk| = 0 (details are presented in Appendix C.).This gives the mean-field as well as the self-consistent equation

III. RESULTS
In the phase diagram depicted in Fig. 2(a)(b), the abscissas represent the system size 1/L x and exciton density n ex separately, and the vertical axis is the in-plane electric field strength E. The zero-field band gap ∆ 0 (black line, left axis) and the correlation length ξ (purple line, right axis) estimated by Eq. (F21) are also plotted as functions of system size and exciton density separately in Fig. 2(c)(d).The critical field Ec (solid blue lines) firstly divides the entire region into a locally gapped phase and a metallic phase.The Zener field Ez (dashed orange lines) solved by eEzL = ∆(Ez) marks the onset of Zener tunneling and further separates the locally gapped phase into an excitonic insulating phase and tunneling phase.(c)(d) Zero field band gap ∆ 0 (black lines, left axis) and the correlation length ξ (purple lines, right axis) as functions of system size and exciton density.The red labels above the bottom axis of (c) mark the number of k points used for the corresponding system size.
The parameters in our model are set as m e,h = 0.3m 0 (m 0 is the electron bare mass), = 8 and d = 100a.u.≈ 5.3nm.The momentum space summation in Eq. 3 is restricted in the region |k x,y | < k c = 0.05a.u.≈ 0.94nm −1 .The numerical results are nearly independent of the cutoff k c when k c k F since the BCS-type condensation only occurs in a small range around k F .The size of the system is defined by the spacing of k-mesh as L = 2π/∆k, so the varying of system size is realized by using different sizes of k-mesh.The electrical field is applied in the x direction, and the length of the system perpendicular to it is fixed at L y = 2πN ky /2k c ≈ 266nm (N ky = 80) for numerical convenience.
In Fig. 2(a)(c), the exciton density is fixed at about n ex ≈ 0.84 × 10 4 µm −2 and the number of k points in the x direction is taken as N kx = 40M (M is an integer and some used N kx are marked by red texts above the bottom axis in Fig. 2(c)).On the contrary, in Fig. 2(b)(d), the system size is fixed (k-mesh is fixed at 120 × 80) and the exciton density varies.
As is shown in Fig. 2(c)(d), the correlation length ξ is about 10nm within the range of the parameters we consider.The correlation length ξ is much smaller than the system size L x along the direction of the electrical field, which means tunneling current at the onset of Zener tunneling I ∝ e −Lx/ξ is negligible.
To overcome the Zener instability of the energy functional for the electrical field in the range ∆/eL x ∼ ∆/eξ, the polarization Hamiltonian h P k and the polarization energy are always evaluated on the coarse 40×80 mesh.For an original 40M ×80 k-mesh, this is equivalent to dividing the system into M copies with size L x = L 0 ≈ 133nm.Thus the Zener tunneling process whose tunneling length satisfies M L 0 > > L 0 > ξ is ignored.This approximation is reasonable since the tunneling probability e − /ξ for such process is smaller than e −L0/ξ ≈ 10 −3 .
The blue lines in Fig. 2(a)(b) represent the critical field E c accounting for the many-body breakdown of the excitonic gap, which divides the entire region into a metallic phase and a locally gapped phase.By solving eE z L x = ∆(E z ), the minimum field E z required for Zener tunneling is obtained and plotted by the orange lines and further separates the locally gapped phase into an excitonic insulating phase and a tunneling phase.In the excitonic insulating phase, the system is fully gapped, and no current flows.In the tunneling phase, an exponentially small Zener tunneling current appears while the system is still locally gapped.In the metallic phase, the excitonic gap is destroyed, the system becomes highly conductive and the resistivity-temperature (R − T ) curve becomes typical metallic.
To understand the breakdown phase transition, let's examine the stability of the local minimum, which is realized by calculating the second-order derivatives (Hessian matrix) of the energy functional.
Assume we are in the region of insulating state, so the local minimum |vk; E could be found by our selfconsistent procedure.The self consistent equation at the mean field solution reads h M F k [|vk; E ; E]|ik; E = ξ ik,E |ik; E , where |ck; E , |vk; E are conduction and valance bands with band energies ξ ck,E > ξ vk,E .At the local minimum, the trial HF state could be reparameterized as This parametrization is complete and unconstrained (f k is an arbitrary complex-valued function).Then the total energy becomes a functional of f k , i.e.
, and the Hessian matrix at this point is The details and specific expression of H kk could be found in Appendix D.
If f k is small, |v k; E is approximated by |v k; E ∼ |vk; E + f k |ck; E .Such a form can be viewed as the low-energy excitation modes in the variational parameter space illustrated above.By diagonalizing the Hessian matrix, the eigenmodes for the low-energy excitations k H kk f λ k = λf λ k can be obtained.For convenience, the eigenmodes in the following text are normalized by On a 120 × 80 k-mesh with exciton density n ex ≈ 0.84 × 10 4 µm −2 (dashed gray line in Fig. 2(a)(b)), the total energy functional is analyzed, and the results are shown in Fig. 3.In Fig. 3(a), we plot the smallest few eigenvalues λ 0−4 of the Hessian matrix Eq. ( 7) as functions of field strength.By taking trial HF state as |v k; E; θ λi ∝ |vk; E + θ λi f λi k |ck; E , the total energy difference between the trial state and the HF ground state along the directions f λi k in the variational parameter space is evaluated as Using the lowest three eigenmodes f λ0,1,2 k for example, the total energy difference as a function of the electric field E and excitation amplitudes θ λi is plotted in Fig. 3(b)(c)(d).In these plots, the horizontal axes are the amplitudes of those eigenmodes, while different electric field strengths are represented by different color lines.
There is a consistent zero mode λ 0 for any electric field strength, as shown in Fig. 3(a).However, the behaviors of the total energy functional along the direction f λ0 k in Fig. 3(b) indicates that it's not a "breaking down mode" because the high-order derivatives of the total energy functional along this direction are always positive.Such a zero mode is exactly the Goldstone mode related to phase fluctuation of the exciton condensate and accounts for the exciton superfluidity (see details in Appendix E).
The real breaking down direction in parameter space is f λ1 k as shown in Fig. 3(c).When the electric field is small, all eigenvalues of the Hessian matrix(except the Goldstone mode λ 0 ) satisfy λ > λ 1 > 0, which means the solution is indeed a local minimum.As the electric field approaches the critical field strength E c , the eigenvalue of the breakdown mode λ 1 approaches 0 and the excitonic insulator ground state becomes unstable as the local minimum turns into a saddle point.Such a manybody breakdown mechanism is completely different from traditional Zener tunneling and the corresponding critical field strength can be much weaker than the one for Zener tunneling as discussed in the following section.

IV. DISCUSSION
The results in Fig. 2(a)(b) indicate that the critical field E c for the many-body breakdown is nearly independent of the system size and decreases dramatically with the increase in exciton density.This is reasonable since with the increase in exciton density, the binding between electron and hole becomes weaker and the excitonic insulator will eventually turn into a quantum electron-hole plasma state [1,38,39].Then the intersection point of the two critical fields E c and E z (intersection points of the blue and orange lines in Fig. 2(a)(b)) gives a critical length roughly estimated by which separates the n ex − L plane into two regions as illustrated in Fig. 4(a).And the I − V characteristic may behave differently in the two regions.
In the many-body breakdown region below the line of critical length L c (the green color region in Fig. 4(a)), the excitonic gap is disrupted before the onset of interband Zener tunneling.As the electrical field increases, the BCS-type exciton condensation wave function will lose stability and exhibit a typical first-order transition feature.After this transition, the system becomes gapless and highly conductive, and a metallic current I m ∝ eV will flow in the system.
In the Zener region where L L c and Ẽz Ẽex (the orange color region in Fig. 4(a)), a tunneling current will first appear when the in-plane bias voltage exceeds the band gap.Fot gate voltage in the range ∆ ∼ eE z L eV < eE c L, this current is in the form of The exponential factor e −∆L/eV ξ is the WKB tunneling probability and power term (eV − ∆) 3/2 arises from the density of states of the tunneling channels in 2D systems (details could be found in Appendix F).Different from the metallic current, the tunneling current exhibits a different R − T characteristic, i.e. the current increases when the temperature rises.
The tunneling current persists until the field strength reaches the critical field of many-body breakdown, after which the excitonic gap disappears and a metallic current I m ∝ eV appears replacing the Zener tunneling current I z .However, even at the critical field E c , the tunneling current I z (E = E c ) ∝ e −Lc/ξ in the BCS limit is still exponentially small as the critical length L c is nearly two orders larger than correlation length ξ as is shown in Fig. 4(a).This means a switching phenomenon of the I − V curve is still observable even in the Zener region.
From the discussion above, the typical I − V characteristics in the two regions are schematically illustrated in Fig. 4(b)(c).
We notice that Sugimoto et al. [26] also proposed a breaking down mechanism in correlated insulators which has a threshold field much smaller than that for Zener breakdown.However, the mechanism in their work is distinct from the many-body breakdown mechanism proposed in our work.The many-body breakdown is intrinsic for an excitonic insulator while the critical field in their work is related to the extrinsic relaxation time.Besides, the typical I −V curve for an excitonic insulator as illustrated in Fig. 4(b)(c) has size dependence which is already observed by the experiments of Yang et al. [40].
At last, the many-body breakdown mechanism is a breakdown of the electronic band structure and has nearly no influence on the lattice, which means the breaking-down process is reversible and the switching phenomenon of the I − V characteristic is promising for practical usage.We first prove that under a uniform electric field, the many-body state will keep its lattice translation symmetry at all times.
A many-body state |Ψ ; t is said to have lattice translation symmetry if and only if the wavefunction satisfies where N e is the total number of electrons and R 0 is arbitrary lattice vector.
The many-body Schrödinger equation in length gauge (using a scalar field ϕ = eE • r to include electric field) is written as which seems to break lattice translation symmetry.However, by taking gauge transformation of the electric field ∂ t A(t) = −E and defining

A2) we find that the Schrödinger equation for |Ψ
which keeps the lattice translation symmetry.So starting from a many-body state |Ψ 0 with lattice translation symmetry, the many-body state |Ψ A ; t as well as |Ψ E ; t will have lattice translation symmetry at any time: When treating a static uniform electric field, as long as the field is adiabatically turned on, a trial HF state with lattice translation symmetry could be safely assumed.For insulators, this state is written as where n e is electron per cell and |vac. is vacuum state.c † nk is creation operators of Bloch electron with wavefunction where V = N v c is the volume of the system, N is the number of unit cells and v c is cell volume.
As electron creation operators, c † nk should satisfy which forces the corresponding Bloch functions to be orthonormal, i.e.
ψ mk |ψ nk = dr ψ * mk (r)ψ nk (r) = δ mn δ kk , (A8a) Appendix B: Polarized HF Energy Functional In this part, a general form of the polarized HF energy as a functional of occupied Bloch states will be derived.
Using field operator Ψ(r), the second quantization form of the single particle (kinetic and potential energy), polarization, and interaction Hamiltonians are written as Matrix elements of the single-particle density operator ρ under position basis are calculated as Then its k-dependent counterpart is defined by Need to notice that the single particle Hilbert space H of ρ is all kinds of functions while the Hilbert space H k of ρk is only the cell-periodic functions.That's why the prefactor N , number of cells, appears in the definition of ρk in Eq. (B5).And we will see the single-particle and interaction energies could be expressed as functionals of ρk and therefore functionals of occupied states |u nk .The single-particle part is where ĥ0 k = e −ik•r ĥ0 ( p, r)e ik• r = ĥ0 ( p + k, r) is the k-dependent single particle Hamiltonian acting on cellperiodic functions with matrix elements Similarly, the interaction part is evaluated with the help of Wick's theorem q) drdr e iq•(r−r ) ρ(r, r)ρ(r , r ) where V (q) ≡ dr V (r)e −iq•r is the Fourier transformation of V (r).The first part in Eq. (B8) is the Hartree energy and is simplified as where G is reciprocal vector.The e iG• r term in Eq. (B9) should be understood as a single particle operator that acts on |u nk as (B10) The second part in Eq. (B8) is the Fock energy The polarization energy can't be expressed by density operator ρk but is still a functional of occupied states This result is consistent with the Berry phase definition of polarization.For a finite-size system with periodic boundary conditions, the polarization, and the polarization energy should be written with the discrete form of Berry phase as [30] where |∆k | = 2π/L and is along the direction of electric field.The overlap matrix S is defined as Appendix C: Mean Field Hamiltonian and Self-consistent Equation The total energy as a functional of occupied bands {|u nk } ne n=1 is written as and the stationary state is found by minimize E tot with constraints By introduce Lagrange multipliers ξ nk , the constrained minimization of E tot is transformed into an unconstrained minimization of (C3) Let's calculate the unconstrained derivatives of F with respect to u nk |.We first show that The single-particle, Hartree, and Fock energy functionals all take this form thus are easily evaluated From the expression above, we could define the Hartree and Fock Hamiltonian as As functionals of gauge invariant single-particle density operator ρk , the Hartree and Fock Hamiltonian defined in Eq. (C8)(C9) are also invariant under k-space gauge transform of the occupied bands.
As for the polarization term, we start from the discrete form Eq. (B13) and take the thermodynamic limit later.
The unconstrained derivatives of E P is where abbreviation k σ = k + σ∆k is used for simplicity.Denote |D nk = δE P /δ u nk |.Easy to see that So polarization Hamiltonian could be defined as and satisfies Before processing, one should verify that this definition of polarization Hamiltonian is a gauge invariant.By denoting , the polarization Hamiltonian is written in a more neat form A k-space gauge transformation (U k ) ne×ne on occupied bands will transform Φ k into Φ k U k and the polarization Hamiltonian becomes which is invariant.It's easier to see this gauge invariance in the thermodynamic limit L → ∞ and dk = ∆k → 0. In this limit So and the polarization Hamiltonian in the thermodynamic limit is written as The thermodynamic limit expression Eq. ( C18) is only a functional of the gauge invariant ρk and thus is also a gauge invariant.Finally, minimization of F [|u nk ; E] gives us the selfconsistent equation where the mean-field Hamiltonian is Appendix D: The Hessian Matrix Assume E < E c , and the self consistent equation has solutions The valence band |vk; E is chosen as the one with lower band energy, i.e. ξ vk;E < ξ ck;E , The E label in wavefunctions and band energies means they are converged solutions.
At the converged point (local minimum of the total energy functional), the trial HF state could be reparameterized as where f k is arbitrary complex-valued function defined on Brillouin zone.This parametrization is unconstrained and complete, and the total energy then becomes functional of f k as By writing , where f k,r and f k,i are real variables, the Hessian matrix is defined as For simplicity, the E label will be omitted in the following derivations.
We first calculate the derivatives of |v k with respect to f k,r/i for further usage.
At f k = 0, they are simplified as The second order derivatives of |v k at f k = 0 are The first order derivative of ε tot defined by Eq. ( 3) is We use the definition of mean-field Hamiltonian h M F k [|vk ]|vk ≡ Vδε tot /δ vk| for the last equality in Eq. (D11).At f k = 0, the first-order derivative is just which is consistent with the fact that |vk is a local minimum.Then let's evaluate second-order derivatives of ε tot So the final task is to evaluate the derivatives of h M F k with respect to f k,r/i .The Hartree one is The Fock one is As for the polarization term, use the fact that we have and Appendix E: The Goldstone Mode The many-body Hamiltonian Eq. ( 1) is invariant under gauge transformations of the electron creation operators: c † ek → e iφe c † ek , c † hk → e iφ h c † hk .This U (1) × U (1) symmetry corresponds to the charge conservation in each layer.
After this gauge transformation, the valance band electron creation operator becomes which gives a new trial wavefunction where φ = (φ e + φ h )/2, φ ex = (φ e − φ h )/2 are related to the conservation of total charge and exciton number respectively.The relative density matrix ρ = ρ−ρ 0 transforms into or equivalently, ρ ss k = e i(φs−φ s ) ρss k .Besides, the overlap matrix S(k, k) = vk|vk becomes the valance band fluctuation.To see this, let's rewrite Compare Eq. (D2)(E5) we find the corresponding Goldstone mode in parameter space is just expressed as The overlap between the Goldstone mode f Gs k and the zero mode f λ0 k of the Hessian matrix Eq. ( 7) is calculated as and plotted in Fig. 5.The results show that I is equal to 1 in numerical precision, which means the zero mode of the Hessian matrix is indeed the Goldstone mode f Gs k discussed in this section.
Appendix F: The Inter-band Zener Tunneling Consider the inter-band tunneling problem of the 2D continuous model ĥ where the barrier potential V (x) is defined as For a given tunneling energy E, the Schrödinger is ĥ|Ψ; E = E|Ψ; E . (F3) Since the electrical field is applied only along x-direction, translation symmetry in y direction still holds and k y is a good quantum number.Following Zener [19], we could write the approximated WKB wavefunction as If k(x) is slow varying so that ∂ x k(x) could be neglected, substitute Eq. (F4) into the Schrödinger equation we find that where and µ ex (k y ) = µ 0 ex − k 2 y /m.Solving the secular equation (F5) gives the relation between the complex wavevector k(x) and position x Things are different for µ ex > 0 and µ ex < 0 and should be discussed separately.The condition µ ex (k y ) = 0 gives a critical k y as The tunneling scenario for µ ex (k y ) > 0 (or equivalently k 2 y ≤ k 2 y,c ) is illustrated in Fig. 6(a).As a tunneling state propagating to the right, |Ψ ky ; E should behave like a valance band electron in the region x −L/2 (k ± L states in Fig. 6(a)) and like a conduction band electron in the region x L/2 (k ± R states in Fig. 6(a)).This places a restriction on the tunneling energy (∆ − eEL)/2 ≤ E ≤ −(∆ − eEL)/2, which further demands that eEL ≥ ∆.In other words, the inter-band Zener tunneling only occurs when the in-plane bias voltage exceeds the band gap.
As the electron propagates to the right in the region |x| ≤ L/2, the complex wavevector k(x) will travel from k σ L to k σ R in the complex plane along the line [23] Im[k 2 (x) − mµ ex ] 2 = 0. (F9) Eq. (F9) is just the imaginary part of Eq. (F7) and is solved as The solutions of Eq. (F10) in the complex plane are represented by dashed gray lines in Fig. 6(b).The paths of k σ (x) in the complex plane are also illustrated by solid black arrow lines in Fig. 6(b).This analysis means that k + (x) and k − (x) are two independent tunneling channels.
Need to notice that, the tunneling channel k + (x) only exists for tunneling energy E ≥ −(∆ − eEL)/2 where ∆ = µ 2 ex + ∆ 2 .This is because there is no k + L state in the region x −L/2 when E < −(∆ − eEL)/2 as is shown in Fig. 6 In this case, the classical turning points are x ± = (±∆ /2 − E)/eE.In addition to the region x − ≤ x ≤ x + , the complex wavevector k(x) also has an imaginary part in the region x − ≤ x ≤ x − and x + ≤ x ≤ x + which is

FIG. 2 .
FIG. 2. (a)(b) Phase diagram as a function of the in-plane electrical field E, system size 1/Lx and exciton density nex.The critical field Ec (solid blue lines) firstly divides the entire region into a locally gapped phase and a metallic phase.The Zener field Ez (dashed orange lines) solved by eEzL = ∆(Ez) marks the onset of Zener tunneling and further separates the locally gapped phase into an excitonic insulating phase and tunneling phase.(c)(d) Zero field band gap ∆ 0 (black lines, left axis) and the correlation length ξ (purple lines, right axis) as functions of system size and exciton density.The red labels above the bottom axis of (c) mark the number of k points used for the corresponding system size.

= 2 2 (direction of f 2
FIG. 3. (a) The smallest five eigenvalues of the Hessian matrix Eq. (7) as a function of the electric field.(b)(c)(d) Total energy difference Eq. (8) as a function of the electric field and excitation amplitudes along the directions f λ 0,1,2 k in the variational parameter space.The excitation amplitudes are used as the horizontal axes while different field strengths are represented by different color lines.These data are generated on a 120 × 80 k-mesh with exciton exciton density nex ≈ 0.84 × 10 4 µm −2 (along the dashed gray line in Fig. 2(a)(b)).

FIG. 4 .
FIG. 4. (a) The critical length Lc where the Zener field Ez equals the critical field Ec for the many-body breakdown in Fig. 2(a) is plotted as a function of exciton density by the dotted blue line, which separates the nex − L plane into two regions, i.e. the many-body breakdown region (green color) and the Zener region (orange color).The correlation length of the excitonic gap ξ given by Eq. (F21) is also plotted by the purple line for reference.The two dashed gray lines mark the paths along which Fig. 2 is generated.(b)(c) I − V characteristics for the excitonic insulator in the many-body breakdown region and Zener region separately.

FIG. 6 .
FIG. 6.(a) Tunneling scenario for µex(ky) > 0. The tunneling channels k ± L → k ± R only exist when in-plane bias voltage overcomes the band gap, i.e. eEL > ∆.Under WKB approximation, the valance band k σ L states in the region x ≤ −L/2 will continuously turn into the conduction band k σ R states in the region x ≥ L/2 as propagating to the right.x± = (±∆/2 − E)/eE marks the classical turning points.(b)The paths of the complex wavevectors k σ (x) in the complex plane are indicated by the black arrow lines.

µ 2 1 0FIG. 7 .
FIG. 7. (a)Tunneling scenario for µex(ky) < 0. In this case, there is no band inversion and the band gap becomes ∆ = ∆ 2 + µ 2 ex .There exists one and only one tunneling channel kL → kR when tunneling energy satisfies |E| ≤ (eEL − ∆ )/2.(b) The path of the complex wavevector k(x) in the complex plane is indicated by the black arrow line.