Second order Josephson effect in excitonic insulators
SSecond order Josephson effect in excitonic insulators
Zhiyuan Sun, Tatsuya Kaneko, Denis Goleˇz,
2, 3 and Andrew J. Millis
1, 2 Department of Physics, Columbia University, 538 West 120th Street, New York, New York 10027 Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY, 10010 Joˇzef Stefan Institute, Jamova 39, SI-1000, Ljubljana, Slovenia
We show that in electron-hole bilayers with excitonic order arising from conduction and valencebands formed by atomic orbitals that transform differently under inversion, nonzero interlayertunneling leads to a second order Josephson effect. This means the interlayer electrical current isrelated to the phase of the excitonic order parameter as J = J c sin 2 θ instead of J = J c sin θ , andthat the system has two degenerate ground states that can be switched by an interlayer voltage. Ina three dimensional stack of alternating electron-hole planes or a two dimensional stack of chains,the second order Josephson coupling can lead to a Weyl semimetal or a quantum anomalous hallinsulator, respectively. A generic order parameter steering effect is demonstrated, whereby electricfield pulses perpendicular to the layers and chains can steer the order parameter phase between thetwo degenerate ground states. The steering is also applicable to the excitonic insulator candidateTa NiSe . Excitonic condensation [1–4] has been experimentallyrealized in electron-hole bilayers (EHB) [5–11] where elec-trons in one layer pair with holes in the other layer toform excitons that condense into a single macroscopicstate. The ideal excitionic insulator (EI) transition in-volves an order parameter that triggers hybridization be-tween the two electronic bands and spontaneously breaksthe U (1) symmetry due to the conservation of their rela-tive charge. In real situations, intrinsic interband (layer)hybridization generically occurs and explicitly breaks the U (1) symmetry. In 1976, Kulik and Shevchenko [12]noted that nonzero interlayer tunneling endowes the EHBwith a Josephson effect similar to that found in super-conductors. The Josephson effect was extended by Refs.[13, 14], was observed in 2000 by Spielman et al. inquantum hall bilayers [15, 16] and explained in detailin Refs. [17–19].If the electron and hole bands are formed by atomicorbitals that transform differently under crystal symme-tries, the intrinsic hybridization vanishes at high symme-try points of the Brillouin zone and is very small nearby.In this case, the excitonic transition typically breaks adiscrete symmetry [4, 20–23]. In this paper, we show thatthis situation can have remarkable consequences if the or-bitals lie at different spatial locations as shown in Fig. 1.Specifically, if one of the relevant orbitals is odd under in-version while the other is even (e.g. p and d orbitals), theordered state sustains a second order Josephson effect. Inan isolated EHB the ground state breaks parity. In threedimensional (3D) stacks of coupled planes or two dimen-sional (2D) stacks of coupled chains (Fig. 2), the EI statebreaks time reversal symmetry to form a Weyl semimetal[24] or a quantum anomalous Hall insulator [25], respec-tively. In all cases the excitonic order parameter maybe ‘steered’ by applied interlayer/chain electric fields viathe AC Josephson effect, enabling controlled switching ofdegenerate ground states. This order parameter steeringapplies as well to the EI candidate Ta NiSe [22, 26–30]. The electron-hole bilayer shown in Fig. 1(a) consists of 𝑥𝑦𝑧𝑗 ! 𝜃𝜋Re[Δ]Im[Δ]𝜃 (b) (c)(a) 𝜙 " 𝑡 Δe $% d 𝜃 𝜙 " 𝑡(ps) FIG. 1. (a) Schematics of the electron-hole bilayer showingelectrons ( − ) and holes (+) and the interlayer current-phaserelation in the excitonic insulating phase. (b) False color rep-resentation of the free energy on the plane of complex orderparameter where lower energy appears bluer. (c) Solid curve:time dependence of order parameter phase after a voltagepulse φ a = − φ e − ( t − t ) /T (dashed curve) computed fromEq. (4) with ∆ p = 14 meV, T = 0 . γ = 0 . p , D = 2and C = 1. two planes labelled 1 and 2 with two-component electroncreation operator ψ † = ( ψ † , ψ † ) from the two bands. TheHamiltonian is H EHB = (cid:88) k ψ † k (cid:18) ξ ( k + A ) + φ e idA z t k + A e − idA z t ∗ k + A ξ ( k + A ) + φ (cid:19) ψ k + (cid:90) drdr (cid:48) V ( r − r (cid:48) ) ρ ( r ) ρ ( r (cid:48) ) (1)where ψ k = (cid:82) dre ikr ψ ( r ), ρ ( r ) = ψ † ( r ) ψ ( r ) is the den-sity, ξ , ( k ) is the kinetic energy describing in-plane mo-tion with ξ dispersing upwards from a minimum value − G/ ξ dispersing downwards from a maximumvalue G/ k = 0, both isotropic. a r X i v : . [ c ond - m a t . s t r- e l ] F e b We distinguish the ‘BCS’ case (
G >
0) where the twobands cross at a Fermi momentum k F with Fermi veloc-ity v F , and the ‘BEC’ case ( G <
0) where they don’toverlap. ( φ i , A i ) is the electromagnetic (EM) potentialat layer i , A = ( A + A ) / A z is the average out of planecomponent and we have set e = c = (cid:126) = 1. We assumethe Hamiltonian is invariant under time reversal ˆ T andin-plane inversion defined as ˆ P : r → − r, ( ψ , ψ ) r → ( ψ , − ψ ) − r where r = ( x, y ). These two symmetriesimply ξ , ( k ) = ξ , ( − k ) and that the intrinsic interlayertunneling satisfies t − k = t ∗ k = − t k . Thus one can write t k = i ∆ p f k where f k is odd under k → − k , ∆ p > k -oddnature. Without loss of generality, we set f k = c f sin k x where c f is a normalization constant chosen such that | f k F | = 1 in the BCS case.To study the excitonic order we write the modelas a path integral and decompose the interac-tion in the electron-hole pairing channel: Z = (cid:82) D [ ψ, ∆ k , A ] e (cid:82) dτdrL ( ψ, ∆ k ,A ) where ∆ k is the Hubbard-Stratonovich field. The excitonic state appears asa saddle point with the order parameter ∆ k = (cid:80) k (cid:48) V k − k (cid:48) (cid:104) ψ † k (cid:48) ψ k (cid:48) (cid:105) where V q is the Fourier transform of V ( r ). For physically reasonable interactions, the energet-ically favored order parameter ∆ e iθ has s -wave symmetry[31] so the k dependence may be neglected. The quasi-particle properties are described by replacing the term e idA z t k + A in Eq. (1) by ∆ e iθ + e idA z t k + A [32]. There isalways an odd parity phonon [20, 26, 33–35] (e.g., shearmotion between the two layers) that couples linearly to∆ and accompanies the order, but may be integrated out.Integrating out the fermions, phonons and the orderparameter amplitude fluctuations one obtains a low en-ergy effective Lagrangian for the order parameter phase: L = 12 ν (cid:34) − ( ∂ t θ + φ a ) + v g ( ∇ θ − A a ) − D ∆ cos(2( θ − A z d )) (cid:35) (2)where ( φ a , A a ) = ( φ − φ , A − A ) / L to second order in t k (assumed small relative to ∆ or temperature), observ-ing that only cross terms involving t and even pow-ers of sin θ survive (see Ref. [38] Sec. I). An inversioneven t k would change this term to ∝ t cos θ , giving riseto the usual Josephson effect [12, 13, 17–19]. The z-dipole density is ρ a = δL/δ ( ∂ t θ ) = − ν ( ∂ t θ + φ a ) andEq. (2) should be supplemented by the electric field en-ergy (cid:80) q φ a ( q ) / (2 V eff ( q )) representing the dipole-dipoleinteractions V eff ( q ) = (1 − e − dq ) V q [36, 39]. In the BCSweak coupling limit ∆ (cid:28) G at zero temperature, thecoefficients of Eq. (2) have simple ∆-independent forms: D = 2 is the space dimension, ν is the density of statesin the normal state at the band crossing energy and the bare phase mode velocity is v g = v F / √ t k is zero, Eq. (1) conserves the charge in eachplane and there is a continuous family of excitonic phasesparametrized by θ , as manifested by the U (1) symmetryunder transformation θ → θ + θ of the first two termsof Eq. (2). A non-zero t k gives rise to the third termwhich reduces the U (1) invariance to ˆ P , a Z symme-try and implies that there are two degenerate excitonicphases characterized by θ = 0 , π (Fig. 1(b)). The exci-tonic order breaks ˆ P , implying a non-vanishing in-planeelectrical polarization whose sign is opposite for θ = 0 , π [21]. Calculations similar to Ref. [31] show that in theBCS limit, the polarization density relative to the nor-mal state is P = P (cid:104) − tan (cid:16) ArcTan | ∆ p ∆ | (cid:17)(cid:105) Sgn[∆] / P = 2 k F /π . Since P reverses direction as ∆ re-verses sign, a measurement of P by an electrical circuitcan distinguish the two ground states.In spinful systems both singlet and triplet excitoniccondensates may be defined. The triplet case exhibitsspin instead of charge polarization. In the pure electronicsystem the two phase are degenerate at the Hartree-Focklevel, but electron-lattice coupling favors the singlet state[4, 20, 33] (see Ref. [38] Sec. V). We focus on the morecommonly studied singlet phase here. Second order Josephson effect and order parametersteering—
The interplane current − J z = δL/δ ( dA z ) = νD ∆ sin(2 θ ) ≡ J c sin(2 θ ) (3)is periodic under θ → θ + π in contrast to the usualJosephson effect where the current is periodic only under θ → θ + 2 π ; the former is thus referred to as a secondorder Josephson effect. Assuming a quadratic band witheffective mass 0 . m e and ∆ p = 10 meV, the critical cur-rent is estimated as J c ≈ /µ m in the BCS limit.To observe the DC Josephson effect, one can source acurrent at one layer and drain it on the other layer, bothon the left side of the device where the in-plane counterflow current J a = νv g ∂ x θ is fixed as the boundary con-dition [15]. From the static limit of the Euler-Lagrangeequation (charge continuity equation) implied by Eq. (2), νv g ∂ x θ = J c sin(2 θ ), the phase decays to the right witha decay length l d = (cid:113) νv g /J c ∼ √ Dv g / ∆ p [40]. Thus ina long junction, only the region within a distance l d tothe contact contributes to the Josephson current [13].To treat the dynamical regime, we focus on spatiallyuniform dynamics which applies to a short EHB withside contacts or to a device with gates covering the wholesample. Eq. (2) in the gauge A = 0 implies1 C ∂ t ( ∂ t θ + φ a ) + γ∂ t θ + 1 D ∆ sin 2 θ = 0 (4)where a C (cid:54) = 1 expresses the effect of dipole-dipole in-teractions (charging energy) and we have added a phe-nomenological damping γ . Thus the time derivative ofan interlayer voltage φ a provides a force that pushesthe phase to increase, meaning that a suitable voltagepulse can switch the system between ground states asin Fig. 1(b)(c). If φ a is applied by side contacts or bygates in immediate adjacent to the bilayer, the externalelectrical circuit controls φ a which is already the totalvoltage across the layers, and one has C = 1 in Eq. (4).If the voltage is applied by gates which do not completelyscreen out the dipolar interactions, C will depend on thedevice geometry and is also affected by the nonlinear dy-namics. To climb the potential hill at θ = π/ ν ∆ /
4, the threshold voltage required for a typi-cal pulse φ a = φ e − ( t − t ) /T is φ c ∼ T ∆ C/D , giving φ c ∼
25 mV in the BCS limit for T = 1 ps, ∆ p = 10 meVand C = 1. In the limit of strong drive ( φ a (cid:29) φ c ), theequation of motion becomes ∂ t θ = − φ a , recovering thefamiliar AC Josephson effect. Note that the switchingfrequency scale 1 /T is upper bounded by the gap ∆.We have assumed that lattice distortions, if present,can dynamically follow the order parameter. In the op-posite limit of slow lattice dynamics (very low phononfrequency), one should fix the lattice distortion. For weakelectron lattice coupling (EPC), the only change is thatthe Z symmetry remains broken and the second min-imum is at higher energy. For larger EPC the secondminimum no longer exists. Thus fast phase steering canreveal the strength of EPC. Beyond bilayers—
The combination of excitonic orderand odd parity tunneling leads to rich physics also inthe 3D/2D systems resulting from stacking the electron-hole bilayers/chains as shown in Figs. 2(a),(b). Thestacking is along z and the conjugate wavevector is k z ∈ ( − π, π ] / (2 d ). The model is invariant under trans-lations by the z-direction lattice constant 2 d and reflec-tion z ↔ − z with respect to a plane containing eitherthe electron or holes. We specialize to short rangeddensity-density interaction g such that excitonic order∆ i / only links adjacent layers as shown in Fig. 2, andconsider mean field solutions where the amplitude ∆ isspatially uniform but allow for the phases θ , on thetwo bonds to be different. We define the symmetric andantisymmetric phase combinations θ s,a = ( θ ± θ ) / θ s ∈ ( − π, π ] , θ a ∈ [0 , π ). In the mo-mentum basis of field operators ψ † k = (cid:16) ψ † k , ψ † k (cid:17) = (cid:82) dr (cid:80) j e i ( k ⊥ r + k z j d ) (cid:0) ψ j ( r ) , e ik z d ψ j ( r ) (cid:1) where k ⊥ isthe momentum along the planes/chains, the Lagrangianreads L = (cid:80) k ψ † k ( ∂ τ + H k ) ψ k + g | ∆ | with the mean fieldHamiltonian H k = (cid:18) ξ ( k ⊥ ) ∆( k ) − i ∆ p f k cos dk z ∆( k ) ∗ + i ∆ p f k cos dk z ξ ( k ⊥ ) (cid:19) (5)where the ∆ p term is the intrinsic interlayer tunneling t k and the order parameter is∆( k ) = e iθ a ∆ cos( dk z + θ s ) . (6)The spatially uniform electric field enters through k → k + A , including the ∆( k ) term, meaning that we are Δ !" 𝑒 ! !" Δ !$ 𝑒 ! ! 𝑡 % 𝑡 % 𝑖1𝑖2𝑖 − 1, 2𝑡−𝑡 (a)(b) 𝑘 & 𝑘 ’ 𝐸 % (c)(d) Δ ! 𝑡 𝜃 ( 𝐸 ’ 𝑗/𝑗 ) 𝐺d𝐸 ’ Δ * Δ FIG. 2. (a) Schematic of the 3D stack of alternating elec-tron (blue) and hole (blank) planes with pairing order pa-rameters labeled. (b) An example of the 2D stack of alter-nating electron and hole chains. The green and blue dotsrepresent atomic orbitals forming the conduction and valencebands. Their different parities lead to asymmetric inter chainhoping t/ − t . Arrows represent the spontaneous circulat-ing currents. (c) The ground state band dispersion of the2D stack. (d) The order parameter phase dynamics (blackcurve) and the Josephson current (blue curve) induced by anelectric field pulse E z ( t ) = E max e − ( t − t ) /T (red curve) im-plied by Eq. (10), with ∆ = 10∆ p , E max = 3 . p /d and T = 0 . / ∆ p . using a different gauge than in our analysis of the singlebilayer.At ∆ p = 0, the energy is independent of θ and θ andthe symmetry broken excitonic phase has nodal points(2D) or line (3D) at k z = ( π/ − θ s ) /d in its quasi par-ticle spectra. A nonzero intrinsic tunneling ∆ p reducesthe symmetry to ˆ T and ˆ P , and the excitonic ground stateturns out to spontaneously break ˆ T instead of ˆ P , corre-sponding to ( θ a , θ s ) = (0 , ∓ π/ θ i = θ i = ± π/ p (see Ref. [38] Sec. II). Fixing θ a = 0 and inthe gauge φ = 0, one finds: L = K [ ˙ θ s + d ˙ A z , A x ] + c ν ∆ cos 2 θ s + F (7)where K is the kinetic term that vanishes in the staticlimit, F ( | ∆ | ) is the ground state free energy withoutinterlayer tunneling, and we have neglected constant O (∆ ) terms. The cos 2 θ s term means a ‘Josephson’ cur-rent j z = j c sin 2 θ s where j c = 2 dc ν ∆ . In the BCSlimit, c ν ∼ ν where ν is the density of states in the nor-mal states at the band crossing energy. In the deep BECregime ( G (cid:28) − ∆), c ν ∼ ν ∆ /G where ν is a char-acteristic density of states of the bands in the normalphase. In the ground state, the total electrical polar-ization is zero. But since ˆ T is broken, there are alter-nating in-plane currents j inter ,a = (cid:104) (cid:80) k ( ∂ k t k ) sin( dk z ) σ (cid:105) ( ∼ nv F ∆∆ p ε F ln Λ∆ in the BCS limit where Λ is an energycutoff) as shown in Fig. 2(b). The quasi particle disper-sion is E k = ± (cid:113) ξ k + | ∆( k ) | + ∆ f k cos ( dk z ) where wehave set ξ ( k ) = − ξ ( k ) = ξ k for simplicity. Note thatthis ground state is linearly stable to lattice distortion.We next study the Gaussian fluctuations of the phase θ s around the ground state: θ ≡ θ s − θ where θ = ± π/ θ s and the EM fields A x/z . In the low energyregime ω (cid:28) ∆ p and long wavelength limit q = 0, it reads S s = − (cid:88) ω c ( ω ) (cid:0) θ + dA z (cid:1) − ω (cid:0) θ + dA z (cid:1) ω + (cid:90) dtdr (cid:16) c θ + σ h ( θ + dA z ) ∂ t A x /d (cid:17) + S A x (8)neglecting terms sub leading in ∆ p . The first two termsare the kinetic and potential energy costs of phase fluctu-ations where c ( ω ) is the kinetic kernel that vanishes inthe static limit and c = 2 c ν ∆ for ∆ p (cid:28) ∆. The thirdterm means the anomalous hall conductivity σ h whichcan also be written into an ‘Axion’ form [24]. The lastterm is the bare optical response in x direction. Inte-grating out the phase in Eq. (8), one obtains the EMaction S EM = − i (cid:80) ω ωσ ij A i ( − ω ) A j ( ω ) with the opticalconductivity σ = 1 c − c (cid:18) i σ d ω + ( c − c ) σ x − c σ h c σ h ic c d /ω (cid:19) + σ h (cid:15) xz (9)where i/j takes the value of x, z . In the DC limit, thestacks do not exhibit superconductivity along z since c ( ω ) vanishes faster than ω , but have a hall response σ h due to broken ˆ T . As frequency increases, the firstterm takes effect which contains a resonance at the barephase mode frequency set by c ( ω ) = c .Deep in the BEC regime ( G (cid:28) − ∆), the sys-tem is topologically trivial with the minimal gap being √ G + 4∆ and σ h vanishes, although there is a nonzeroAC hall response ∼ ∆ p which can be measured by Kerrrotation (neglected in Eq. (8)). The kinetic energy kernel c ∼ ν ω ∆ /G gives the phase mode gap ω ∼ ∆ p . TheBCS regime ( G >
0) is topological nontrivial which wediscuss separately for the stacks of chains and layers.In the 2D stack of electron-hole chains, the quasiparti-cle dispersion is fully gapped with massive Dirac points at( k x , k z ) = ( ± k F ,
0) with mass ∆ p , as shown in Fig. 2(c).The Chern number of the valence band is C = Sgn[ θ s ]with each massive Dirac point contributing C/
2, so thatthe system is a quantum anomalous Hall insulator [25]with quantized Hall conductivity σ h = Ce /h and chiraledge states (see Ref. [38] Sec. II). If one views k z as anadiabatic variable in time, this system can be viewed asa topological charge pumping process [31]. The kinetickernel c = ν ω ∆ / ∆ p implies that the bare phase modefrequency ω ∼ ∆ p (cid:112) ∆ p / ∆ is below the minimum gap.In a 3D stack of electron-hole layers, the quasiparti-cle dispersion has Weyl nodes at k = (0 , ± k F , c ∼ ν ∆∆ p ω (cid:16) i Sign( ω ) + π ln ∆ p | ω | (cid:17) . The anomalous hallconductivity is σ h = Sgn[ θ s ] k F π e /h whose value canbe understood by noting that there are k F /π quantumanomalous hall k x − k y planes between the two Weylpoints (see Ref. [38] Sec. V for the effect of spin degreesof freedom). Order parameter steering by light—
In all the systemsconsidered above, the order parameter can steered be-tween the two degenerate ground states by electric fieldsperpendicular to the layers/chains: from ground state | g (cid:105) to ˆ P | g (cid:105) for the bilayer and to ˆ T | g (cid:105) for the 3D/2D stacks.The order parameter steering follows the spirit of ACJosephson effect: the phase θ s enters the kinetic term inEq. (7) together with the vector potential as θ s + dA z .This term has different forms in different regimes. Forexample in the 2D stacks, it behaves as K ∼ ν | ∆∆ p | ˙ θ s inthe slow limit of ˙ θ s (cid:28) ∆ p and as K ∼ ν | ∆ | θ s ˙ θ s in themoderately fast case of ∆ p (cid:28) ˙ θ s (cid:28) ∆ where we havesuppressed A z for notational simplicity. Nevertheless,upon strong electric field E z such that the free energypotential cos 2 θ s can be neglected, the equation of mo-tion all reduces to ˙ θ s = d ˙ A z = − dE z , i.e., it provides aforce to rotate the phase θ s so as to switch the systembetween the two ground states θ s = ± π/
2, as shown byFig. 2(d). The pulse needed to exactly deliver such aswitch is thus d (cid:82) E z ( t ) dt = π . For a pulse duration of1 ps and d = 1 nm, the field needed is E z ∼ × V / cm.For weaker fields such that the free energy potentialmatters, the dynamics depends on the time scale. In thecase of ∆ p (cid:28) ˙ θ (cid:28) ∆, the equation of motion for θ s implied by Eq. (7) is simply˙ θ s = ∆ | ∆ | sin(2 θ s ) − dE z . (10)The threshold field to climb over the potential barrierand switch the ground states is about E c ∼ ∆ | ∆ | d whichreads E c ∼ V / cm for ∆ p = 10 meV , | ∆ | = 100 meVand d = 1 nm. Interestingly, a static electric field candrag a continuous phase winding if the self consistentcondition ˙ θ s (cid:29) ∆ p is satisfied such that Eq. (10) holds,requiring a threshold field E dc ∼ ∆ p /d = 10 V / cm, un-like the case of superconductor stacks where no thresh-old field is needed. The resulting AC Josephson current j = j c sin 2 θ s thus oscillates with frequency 2 dE z [41].The combination of switchability and readout makesthese systems potential memory devices that can be writ-ten by electric field pulses either in the THz regime fromphotons or in the GHz regime from electrical circuits. Material realizations—
The bilayer could be homo-bilayers made of, e.g., phosphorene [42, 43] which hasa direct bandgap at Γ point whose conduction and va-lence bands transform differently under certain reflection.Hetero-bilayers may also be considered. Gating bilayerphosphorene [44] could bring the conduction band of onelayer and valence band of the other layer closer in energyand even lead to their crossing, possibly achieving theEI phase. The 3D and 2D stacks may be either naturalcrystals such as WTe (in monolayer form, a 2D stack ofchains) [45–47] or artificial structures. Identifying mate-rial realizations of these topological excitonic insulators[47–52] is an important future research direction.The light induced order parameter steering also ap-plies to the EI candidate Ta NiSe . Its excitonic physicsis captured by a Ta-Ni-Ta chain [22, 26] along x wherethe conduction bands from Ta d -orbitals are even un-der reflection σ ⊥ : x → − x while the valence band fromNi d -orbitals is odd. The EI state breaks σ ⊥ sponta-neously. Under applied electric field vertical to the chain,the induced phase dynamics is still described by Eqs. (4)and (7) which predict that a strong photon pulse canswitch the system between the two degenerate groundstates (See Ref. [38] Sec. IV). This may have alreadybeen observed in a recent experiment [30]. From theinterband hybridization ∆ p ≈
36 meV and the chain width d ≈ E c ∼ T ∆ /d ∼ V / cm for 0 . ACKNOWLEDGMENTS
Z. S. and A. J. M. acknowledge support from the En-ergy Frontier Research Center on Programmable Quan-tum Materials funded by the US Department of Energy(DOE), Office of Science, Basic Energy Sciences (BES),under award No. DE-SC0019443. T. K. is supported byGrants-in-Aid for Scientific Research from JSPS (GrantNos JP18K13509) and by the JSPS Overseas ResearchFellowship. D. G. is supported by Slovenian ResearchAgency (ARRS) under Program J1-2455. The FlatironInstitute is a division of the Simons Foundation. Wethank S. Zhang, Y. Murakami and M. M. Fogler for help-ful discussions. [1] N. F. Mott, Philos. Mag. , 287 (1961).[2] L. V. Keldysh and Y. V. Kopaev, Soviet Phys. Solid State , 2219 (1965).[3] D. J´erome, T. M. Rice, and W. Kohn, Phys. Rev. ,462 (1967).[4] B. I. Halperin and T. M. Rice, Rev. Mod. Phys. , 755(1968).[5] L. V. Butov, A. Zrenner, G. Abstreiter, G. B¨ohm, andG. Weimann, Phys. Rev. Lett. , 304 (1994).[6] L. V. Butov, A. C. Gossard, and D. S. Chemla, Nature , 751 (2002).[7] L. Du, X. Li, W. Lou, G. Sullivan, K. Chang, J. Kono,and R. R. Du, Nat. Commun. , 1 (2017).[8] J. I. Li, T. Taniguchi, K. Watanabe, J. Hone, and C. R.Dean, Nat. Phys. , 751 (2017).[9] Z. Wang, D. A. Rhodes, K. Watanabe, T. Taniguchi,J. C. Hone, J. Shan, and K. F. Mak, Nature , 76(2019).[10] J. Eisenstein, Annu. Rev. Condens. Matter Phys. , 159(2014).[11] M. M. Fogler, L. V. Butov, and K. S. Novoselov, Nat.Commun. , 4555 (2014).[12] I. O. Kulik and S. I. Shevchenko, Sov. J. Low Temp.Phys. , 1405 (1976).[13] Y. E. Lozovik and A. V. Poushnov, Physics Letters A , 399 (1997).[14] X.-G. Wen and A. Zee, Phys. Rev. Lett. , 1811 (1992).[15] I. B. Spielman, J. P. Eisenstein, L. N. Pfeiffer, and K. W.West, Phys. Rev. Lett. , 5808 (2000).[16] I. B. Spielman, J. P. Eisenstein, L. N. Pfeiffer, and K. W.West, Phys. Rev. Lett. , 036803 (2001).[17] M. M. Fogler and F. Wilczek, Phys. Rev. Lett. , 1833(2001).[18] A. Stern, S. M. Girvin, A. H. MacDonald, and N. Ma,Phys. Rev. Lett. , 1829 (2001).[19] Y. N. Joglekar and A. H. MacDonald, Phys. Rev. Lett. , 196802 (2001).[20] B. Halperin and T. Rice, Solid State Physics , edited byF. Seitz, D. Turnbull, and H. Ehrenreich, Vol. 21 (Aca- demic Press, 1968) pp. 115 – 192.[21] T. Portengen, T. ¨Ostreich, and L. J. Sham, Phys. Rev.B , 17452 (1996).[22] G. Mazza, M. R¨osner, L. Windg¨atter, S. Latini,H. H¨ubener, A. J. Millis, A. Rubio, and A. Georges,Phys. Rev. Lett. , 197601 (2020).[23] T. Kaneko, Z. Sun, Y. Murakami, D. Golez, andA. J. Millis, “Bulk photovoltaic effect driven by col-lective excitations in a correlated insulator,” (2020),arXiv:2012.09786 [cond-mat.str-el].[24] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev.Mod. Phys. , 015001 (2018).[25] C.-X. Liu, S.-C. Zhang, and X.-L. Qi, Annual Review ofCondensed Matter Physics , 301 (2016).[26] T. Kaneko, T. Toriyama, T. Konishi, and Y. Ohta, Phys.Rev. B , 035121 (2013).[27] Y. F. Lu, H. Kono, T. I. Larkin, A. W. Rost,T. Takayama, A. V. Boris, B. Keimer, and H. Takagi,Nat. Commun. , 1 (2017).[28] D. Werdehausen, T. Takayama, M. H¨oppner, G. Al-brecht, A. W. Rost, Y. Lu, D. Manske, H. Takagi, andS. Kaiser, Sci. Adv. , 1 (2018).[29] K. Sugimoto, S. Nishimoto, T. Kaneko, and Y. Ohta,Phys. Rev. Lett. , 247602 (2018).[30] H. Ning, O. Mehio, M. Buchhold, T. Kurumaji, G. Re-fael, J. G. Checkelsky, and D. Hsieh, Phys. Rev. Lett. , 267602 (2020).[31] Z. Sun and A. J. Millis, Phys. Rev. Lett. , 027601(2021).[32] This t k also couples linearly to a p -wave order parameterand would induce the latter given arbitrarily weak inter-action in the p -wave channel [31]. Therefore, t k should beunderstood as the intrinsic hybridization together withany induced p -wave mean field.[33] T. Kaneko, B. Zenker, H. Fehske, and Y. Ohta, Phys.Rev. B , 115106 (2015).[34] D. Goleˇz, Z. Sun, Y. Murakami, A. Georges, and A. J.Millis, Phys. Rev. Lett. , 257601 (2020). [35] Y. Murakami, D. Goleˇz, T. Kaneko, A. Koga, A. J. Millis,and P. Werner, Phys. Rev. B , 195118 (2020).[36] Z. Sun and A. J. Millis, Phys. Rev. B , 041110(R)(2020).[37] Due to the U (1) breaking intrinsic tunneling, the exci-tonic insulating state is not a perfect superfluid, meaningthat system won’t display in-plane counterflow supercon-ductivity [14]. However, there is a counter flow current j a = n s /mA a in response to an anti symmetric vectorpotential which corresponds to a nonzero in plane mag-netic field [13], leading to diamagnetism, given that themagnetic field is smaller than a threshold value H c de-termined by the strength of interlayer tunneling.[38] See Supplemental Material at [URL will be inserted bypublisher] for details.[39] Taking into account the dipole-dipole interaction of theexcitons, neglecting screening from the gates, the phasemode (density fluctuation of excitons) has the dispersion ω q = (cid:113) (1 + 2 πνd ) (cid:0) ∆ p /D + v g q (cid:1) with ω q =0 being the‘Josephson plasma frequency’.[40] For ∆ p = 10 meV and v g = 10 m / s, the length scale is l d ∼ . µ m in the BCS limit.[41] Note that as θ s crosses the regimes of 0 , π , the systemhas electrical polarizations in x direction similar to theground states of the bilayer. The resulting energy costdue to the surface charge is a marginal effect in 2D butmay be an obstacle in 3D. However, one can use metallicgates to screen out this effect.[42] L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng,X. H. Chen, and Y. Zhang, Nature Nanotechnology ,372 (2014).[43] A. Carvalho, M. Wang, X. Zhu, A. S. Rodin, H. Su, andA. H. Castro Neto, Nature Reviews Materials , 16061(2016).[44] J. Kim, S. S. Baik, S. H. Ryu, Y. Sohn, S. Park, B.-G.Park, J. Denlinger, Y. Yi, H. J. Choi, and K. S. Kim,Science , 723 (2015).[45] Y. Jia, P. Wang, C.-L. Chiu, Z. Song, G. Yu, B. J.,S. Lei, S. Klemenz, F. A. Cevallos, M. Onyszczak,N. Fishchenko, X. Liu, G. Farahi, F. Xie, Y. Xu,K. Watanabe, T. Taniguchi, B. A. Bernevig, R. J.Cava, L. M. Schoop, A. Yazdani, and S. Wu, “Ev-idence for a monolayer excitonic insulator,” (2020),arXiv:2010.05390 [cond-mat.mes-hall].[46] Y. H. Kwan, T. Devakul, S. L. Sondhi, and S. A.Parameswaran, “Theory of competing excitonic orders ininsulating wte monolayers,” (2020), arXiv:2012.05255[cond-mat.str-el].[47] D. Varsano, M. Palummo, E. Molinari, and M. Rontani,Nature Nanotechnology , 367 (2020).[48] Q. Zhu, M. W.-Y. Tu, Q. Tong, and W. Yao, Sci. Adv. , 1 (2019).[49] R. Wang, O. Erten, B. Wang, and D. Y. Xing, Nat.Commun. , 1 (2019).[50] L.-H. Hu, R.-X. Zhang, F.-C. Zhang, and C. Wu, Phys.Rev. B , 235115 (2020).[51] E. Perfetto and G. Stefanucci, Phys. Rev. Lett. ,106401 (2020).[52] Z.-R. Liu, L.-H. Hu, C.-Z. Chen, B. Zhou, and D.-H. Xu, “Topological excitonic corner states and nodalphase in bilayer quantum spin hall insulators,” (2021),arXiv:2101.00923 [cond-mat.mes-hall].[53] A. Altland and B. D. Simons, Condensed Matter Field Theory (Cambridge University Press, Cambridge, 2010).[54] Z. Sun, M. M. Fogler, D. N. Basov, and A. J. Millis,Phys. Rev. Research , 023413 (2020).[55] D. N. Basov, M. M. Fogler, and F. J. Garc´ıa de Abajo,Science , 195 (2016).[56] Z. Sun, ´A. Guti´errez-Rubio, D. N. Basov, and M. M.Fogler, Nano Lett. , 4455 (2015). Supplemental Material for ‘Second order Josephson effect in excitonic insulators’
CONTENTS
Acknowledgments 5References 5I. The electron-hole bilayer 1A. Circulating currents in the state θ = ± π/ U (1) invariant case 42. With P -type intrinsic tunneling 5C. Optical conductivity and hyperbolic phase polaritons 7III. Stacks with s -type tunneling 7IV. Three chain model for Ta NiSe I. THE ELECTRON-HOLE BILAYER
In this section, we show the detailed derivation of the effective low energy Lagrangian (Eq. (2) of the main text)for the order parameter phase of the electron-hole bilayer. We reproduce the Hamiltonian in real space H EHB = (cid:90) drψ † (cid:18) ξ ( p + A ) + φ e idA z t p + A e − idA z t ∗ p + A ξ ( p + A ) + φ (cid:19) ψ + (cid:90) drdr (cid:48) V ( r − r (cid:48) ) ψ † ( r ) ψ ( r ) ψ † ( r (cid:48) ) ψ ( r (cid:48) ) (S1)here for convenience. The Ginzburg-Landau action for order parameter fields and the EM field is obtained byintegrating out the fermions in the Hubbard-Stratonovich decoupled action ( e − S [∆ ,A ] ≡ (cid:82) D [ ¯ ψ, ψ ] e − S [ ψ, ∆ ,A ] ), resultingin S [∆ , A ] = Tr ln [ ∂ τ + H m ] + (cid:90) drdτ | ∆ | g ≡ (cid:90) drdτ L (∆ , A ) . (S2)where H m = (cid:18) ξ ( p + A ) + φ e idA z t p + A + ∆ e − idA z t ∗ p + A + ∆ ∗ ξ ( p + A ) + φ (cid:19) . (S3)The Tr ln means trace of logarithm of the infinite dimensional matrix where k should be interpreted as the spatialderivative − i ∇ acting on the fermion fields, i.e., the matrix is just the kernel ∂ τ + H m for all Fermion fields at all ( r, t )[31, 53]. Performed in Fourier basis, it involves a summation over momenta k , the fermion Matsubara frequencies ω n = (2 n + 1) πT ( n ∈ Z , T is the temperature and we have set the Boltzmann constant to be 1) and a trace oflogarithm of the 2 × L is just the static free energy function.Using a local gauge transformation ( ψ ( r ) , ψ ( r )) → ( e iθ ( r ) ψ ( r ) , e − iθ ( r ) ψ ( r )) that shifts the phase in Eq. (2) ofthe main text to the diagonal terms, it is straightforward to see that the gradients of phase always appear togetherwith the asymmetric EM fields as ∂ t θ + φ a and ∇ θ − A a , formally analogous to superconductors [36, 53, 54]. Sameas superconductors, the coefficients of the leading quadratic terms L = K µν ( ∂ µ θ + A aµ )( ∂ ν θ + A aν ) are just K µν =diag( − ν, n/m ) where n/m = νv g and v g = v F / √ θ in Eq. (2) of the main text is beyond O ( θ ) but comes from expanding the static freeenergy F to second order in ∆ p . In the BCS weak coupling case, F = 1 g | ∆ | − ν π (cid:90) dφ | ∆ + i ∆ p f ( φ ) | ln 2Λ | ∆ + i ∆ p f ( φ ) | = 1 g | ∆ | − ν π (cid:90) dφ (cid:0) | ∆ | + ∆ p f ( φ ) + 2 i ∆ ∆ p f ( φ ) (cid:1) (cid:18) ln 2Λ −
12 ln (cid:0) | ∆ | + ∆ p f ( φ ) + 2 i ∆ ∆ p f ( φ ) (cid:1)(cid:19) ≈ g | ∆ | − ν | ∆ | ln 2Λ | ∆ | − ∆ p ν π (cid:90) dφf ( φ ) (cid:18) ln 2Λ | ∆ | − − sin θ (cid:19) ≈ F ( | ∆ | ) −
14 ∆ p ν (cid:18) | ∆ | − θ (cid:19) (S4)where ∆ is the imaginary part of ∆ and we have made use of f ( φ ) = cos φ and (cid:82) dφf ( φ ) = π .At temperatures close to T c , the free energy reads [54] F = − ν π (cid:82) dφ (cid:0) − (cid:0) ln Λ T (cid:1) | ∆ + i ∆ p f ( φ ) | + c T | ∆ + i ∆ p f ( φ ) | (cid:1) + g | ∆ | where c is an O (1) constant, and the third term in Eq. (3) of the main text becomes − ν ∆ T ∆ p cos(2( θ + A z d )).In the deep BEC regime ( G (cid:28) −| ∆ | ), this term is replaced by − ν G ∆ p cos(2( θ + A z d )) where ν is a characteristicdensity of state in the normal state. A. Circulating currents in the state θ = ± π/ Note that during the dynamics, the system would inevitably pass θ = π/
2, the imaginary order state. In this state,the polarization is zero but there are microscopic circulating currents with zero uniform component. Assuming ∆ ismomentum independent for simplicity, the mean field Hamiltonian is H k = ξ k σ − (∆ p f k + ∆) σ (S5)and the in plane current operator is j x = (cid:88) k ψ † k ( ∂ k H k ) ψ k = (cid:88) k ψ † k ( v p σ − ∆ p ∂ k f k σ ) ψ k = j intra + j inter (S6)where j intra is the current within each layer and j inter is the current between the layers. Integrated over space, itsexpectation value is j x ∝ (cid:88) k (cid:104) ψ k | ∂ k H k | ψ k (cid:105) = (cid:88) k ∂ k (cid:104) ψ k | H k | ψ k (cid:105) = 0 . (S7)However, the j inter component is nonzero since ∂ k f k is even in k and (cid:104) σ (cid:105) k has the same sign at both directions ofmomentum. In the BCS limit and taking f ( k ) = k/k F , it can be directly verified that the currents are j inter = Tr[j inter G] = (cid:88) k,iω n Tr[j inter , k G k , i ω n ] = − ∆ p k F (cid:88) k,iω n Tr[ (i ω n + H k ) σ (i ω n ) − E ]= 2 π ∆ p ∆ k F (cid:88) k E k ∼ nv F ∆∆ p ε F ln Λ∆ (S8)where n is the carrier density in the normal state, ε F = G/ ∼ G is an UV cutoff and we haveassumed ∆ (cid:29) ∆ p . Therefore, there is a nonzero current j inter between the orbitals lying on adjacent layers, whose in 𝑖Δ𝑖Δ−𝑡 𝑡 FIG. S1. Circulating current pattern in the θ = π/ i ∆) drawn on the x − z cross section. Shown isa tight binding example with blue (black) arrows indicating the directions of interlayer (intralayer) currents. plane component must be compensated by j intra = − j inter within the layers to satisfy Eq. (S7). In other words, thisstate has microscopic circulating currents between orbitals while there is no uniform current, as shown in Fig. S1.If the order parameter depends on momentum, the mean field state might appear to have a nonzero uniform currentsince the mean field Hamiltonian is no longer consistent with the current operator ( j (cid:54) = (cid:80) k ψ † k ( ∂ k H k ) ψ k ), i.e., localgauge invariance under ψ ( r ) → ψ ( r ) e iφ , A µ → A µ + ∂ µ φ appears broken. Eq. (S5) will be such an example if the ∆ p term comes from a p -wave excitonic order instead of from intrinsic hybridization: the only contribution to the currentshould now be j intra (cid:54) = 0. To fix this problem, one should note that there is another phase degree of freedom θ p in the p -wave order parameter ( f k → f k − θ p ) [31] that couples linearly with the vector potential, and the true ground statehas nonzero θ p which finally renders a zero net current. II. STACK OF ELECTRON-HOLE LAYERS/CHAINSA. Static free energy and the ground state
We represent dk z by k z wherever possible for notational simplicity. Eq. (7) of the main text is obtained by expandingthe zero temperature static free energy to second order in ∆ p : F = F ( | ∆ | ) −
14 ∆ p (cid:88) k f k cos ( k z ) E k Tr (cid:20) σ − H k σ H k σ E k (cid:21) = F ( | ∆ | ) −
12 ∆ p (cid:88) k f k cos ( k z ) E k (cid:20) E k (cid:0) ξ k + | ∆ | cos 2 θ a cos ( k z + θ s ) (cid:1)(cid:21) = F ( | ∆ | ) − ∆ p (cid:88) k f k cos ( k z ) E k (cid:20) ξ k + 12 | ∆( k ) | (cos 2 θ a + 1) (cid:21) . (S9)where ∆( k ) = e iθ a ∆ cos( k z + θ s ), E k = (cid:112) ξ k + | ∆( k ) | and H k = ξ k σ + ∆ cos θ a cos( k z + θ s ) σ − ∆ sin θ a cos( k z + θ s ) σ . In the BCS weak coupling case | ∆ | , | ∆ p | (cid:28) G , the momentum summation leads to F BCS −−−→ F ( | ∆ | ) − ∆ p ν D (cid:90) dk z π cos ( k z ) (cid:90) dφ k π f φ k (cid:18) (cos 2 θ a + 1) + 2 ln 2Λ | ∆( k ) | (cid:19) = F ( | ∆ | ) − ∆ p ν (cid:18) (cos 2 θ a + 1) + 4 π (cid:90) dk z cos ( k z ) ln 2Λ | ∆( k ) | (cid:19) = F ( | ∆ | ) − ∆ p ν θ a + 1) + h ( θ s )) . (S10)Note that the coefficient ν holds for 3D and should be replaced by ν in 2D. Here we have defined the function h ( θ s ) which is obviously odd in θ s and periodic in θ s with period π . It can thus be Fourier expanded as h ( θ s ) = (cid:80) n ∈ Z a n cos(2 nθ s ), where a = 4 π (cid:90) π dθdk z cos ( k z ) ln 2Λ | ∆ cos( k z + θ ) | = 2 ln 4Λ | ∆ | ,a = 8 π (cid:90) π dθdk z cos(2 θ ) cos ( k z ) ln 2Λ | ∆ cos( k z + θ ) | = − . (S11)Therefore, to leading Fourier expansion the free energy reads F = F ( | ∆ | ) + ν p (cid:18) − cos(2 θ a ) + cos(2 θ s ) − − | ∆ | (cid:19) (S12)in the BCS weak coupling case.For systems close to the transition temperature T c such that T c (cid:29) | ∆ , ∆ p | , or in the BEC case G (cid:28) −| ∆ , ∆ p | , thefree energy can be expanded in powers of ∆ k = e iθ a ∆ cos( k z + θ s ) − i ∆ p f k cos k z as F = (cid:88) k (cid:2) c | ∆ k | + c | ∆ k | + O ( | ∆ k | ) (cid:3) . (S13) (a) (b) (c) FIG. S2. The free energies (a) f ( θ a , θ s = π/
2) (b) f ( θ a = 0 , θ s ) in arbitrary units for several values of (cid:15) F for the 2D stack.(c) Blue curve is the optimized order parameter ∆( k F , k z ) as a function of k z . Red curve is the imaginary intrinsic tunneling i ∆ p f k F cos k z . The band dispersion is ξ = − ξ = − cos k x +1 − (cid:15) F . The s -wave gap magnitude is fixed at ∆ = 0 .
1, the intrinsichybridization is ∆ p = 0 .
01. The dependence on the phases are similar in both BCS ( (cid:15) F >
0) and BEC ( (cid:15) F <
0) regimes.
One has c ∼ − /T, c ∼ /T for ξ k (cid:28) T and c ∼ /ξ k , c ∼ /ξ k for ξ k (cid:29) T . The first term of Eq. (S13) reads | ∆ k | = ∆ cos ( k z + θ s ) − p f k sin θ a cos( k z + θ s ) cos k z + ∆ p f k cos k z . After being summed over k , it has no θ a or θ s dependence since the cross product term ∼ ∆∆ p sums to zero due to f k being odd. The second term ofEq. (S13) gives phase dependence: F = F + 2∆ ∆ p (cid:0) θ a + 1 (cid:1) π ) D (cid:90) dk z dk D − ⊥ c ( k ⊥ ) f k ⊥ cos ( k z + θ s ) cos k z = F + c θ ∆ ∆ p (cid:0) θ a + 1 (cid:1) (2 cos θ s + 1) (S14)where c θ = π ) D − (cid:82) dk D − ⊥ c ( k ⊥ ) f k ⊥ . In the BCS limit and close to T c , one has c θ ∼ ν/T . In the BEC regime atzero temperature, one has c θ ∼ ν /G where ν is the characteristic density of states of the band ξ ( k ) in the normalphase.In any case, the intrinsic tunneling has reduced the symmetry group from U (1) × U (1) to Z × Z , i.e., inversion ˆ P that maps ( θ a , θ s ) to ( θ a , − θ s + π ) and time reversal ˆ T that maps ( θ a , θ s ) to ( − θ a , − θ s ). Since the free energy minimalie at ( θ a , θ s ) = (0 , ∓ π/ T . Numerically exact results (non perturbative in∆ p ) for the 2D stack are shown in Fig. S2.We take the BCS weak coupling limit to describe the ground state and edge states. In the 2D stack, the groundstate is a quantum anamolous hall insulator with chiral edge states. For example, along the edges parallel to z , theedge states have dispersion E ( k z ) = ± ∆( k z ) and wave function ψ ± ( k z ) = (1 , ±
1) sin( k F x ) e ∓ ∆ p cos( k z ) x e ik z z with the ± sign denoting the left/right edge [31]. In the 3D stack of layers, the ground state is a Weyl semi metal with fermiarc surface states. On the y − z surface, the surface states have dispersion E ( k y , k z ) = ± ∆( k z ) and wave function ψ ± ( k y , k z ) = (1 , ±
1) sin( k F x ) e ∓ ∆ p f k cos( k z ) x e i ( k y y + k z z ) . Here f k means the function f evaluated on the fermi surfaceat ( k x , k y ) = ( (cid:113) k F − k y , k y ).Note that in the limit of only one bilayer with periodic boundary condition, k z can only take the value 0, and theground state favors the order parameter ∆( k ) = ± ∆ cos( k z ) which does not break ˆ T . As number of layers is increasedbeyond two bilayers, the ˆ T breaking order parameter ∆( k ) = ± ∆ sin( k z ) starts to have lower energy. B. Kinetic terms U (1) invariant case In this subsection, we discuss the ∆ p = 0 limit where the system has U (1) × U (1) invariance corresponding tovarying θ and θ , and the ground state has ∆ i = ∆ i = ∆. Without the EM field and assuming ξ = − ξ = ξ ,the quasiparticle dispersion is E = ± (cid:112) ξ k + ∆ cos ( k z ) with Dirac nodal lines at ( k ⊥ , k z ) = ( k F , ± π/ L a = 12 ν (cid:2) − ( ∂ t θ a + φ a ) + v g ( ∇ θ a − A a ) (cid:3) (S15)where ∇ means the in plane gradient and ( φ a , A a ) = ( φ − φ , A − A ) / A a = 0, including the electrical field energy L EM = πd φ a (charging energy of the bilayers viewed as capacitors) and integrating it out, one obtains the Lagrangianof Coulomb renormalized phase modes: L a = −
12 11 /ν + 2 πd ( ∂ t ϕ a ) + 12 νv g ( ∇ ϕ a ) . (S16)The low energy Lagrangian for symmetric fields is L s ∼ νv F E s + ν | ∆ | ( θ s − dA z ) ∂ t ( θ s − dA z ) (S17)where E s is the in-plane electric field. The coefficient of the first term is determined by in-plane optical conductivityneglecting dissipative terms. As expected, there is no z-direction current response since A z drops out after integratingout θ s .
2. With P -type intrinsic tunneling To second order in ∆ p , we write the action for quadratic fluctuations of θ ≡ θ s − θ around the ground state as S s = − (cid:88) ω c ( ω ) (cid:0) θ + dA z (cid:1) − ω (cid:0) θ + dA z (cid:1) ω + (cid:90) dtdr (cid:16) c θ + 2 χθ∂ t A x /d + σ h A x ∂ t A z (cid:17) + S A x . (S18)The kinetic kernel c ( ω ) in the first term is the ω dependent part of the correlation function c ( ω ) ≡ χ θ,θ = ∆ χ cos( k z ) σ , cos( k z ) σ ( ω ) = − (cid:88) k cos k z E − ∆ sin k z ( ω − E ) E . (S19)For the 3D stack c D ( ω ) − c D (0) ≈ πi ∆ (cid:88) k E − ∆ sin ( k z ) E ( δ ( ω − E ) − δ ( ω + 2 E )) + Real Part ≈ ν ∆∆ p ω (cid:16) i Sign( ω ) + π ln ∆ p | ω | (cid:17) ω (cid:28) ∆ pπν ( i ∆ ω + π ω ) ∆ p (cid:28) ω (cid:28) ∆ ν ∆ ( i Sign( ω ) + π ln ∆ | ω | ) ω (cid:29) ∆ (S20)where the real part is obtained from the imaginary part through Kramers-Kronig relation and we have neglectedsubleading terms. For the 2D stack, there is no dissipation when ω < ∆ p and one can expand to O ( ω ): c D ( ω ) − c D (0) O ( ω ) −−−−→ ∆ (cid:88) k cos ( k z ) ω E − ∆ sin ( k z )4 E (cid:28) ε F −−−−→ ∆ ω ν dπ (cid:90) dk z cos ( k z ) (cid:18) / k ) − ∆ sin ( k z ) 1 / k ) (cid:19) = ∆ ω ν π (cid:90) − dt (cid:112) − t (cid:32) / t + ∆ p (1 − t ) − ∆ t / (cid:0) ∆ t + ∆ p (1 − t ) (cid:1) (cid:33) = ∆ ω ν π (cid:32) π − ∆ p + | ∆∆ p | ∆ p (∆ − ∆ p ) − π − ∆ p + 2 | ∆∆ p | + | ∆∆ p | (∆ − ∆ p )∆ ∆ p (∆ − ∆ p ) (cid:33) ∆ p (cid:28) ∆ −−−−→ ν | ∆ || ∆ p | ω (S21)In the regime 2∆ p (cid:28) ω (cid:28) ∆, the minimal gap ∆ p can be neglected such that there are Dirac nodes at ( k x , k z ) =( ± k F , − θ s + π/
2) which renders the dynamics of θ s dissipative. The kinetic kernel for θ s becomes c D ( ω ) − c D (0) = ∆ πi (cid:88) k ξ E ( δ ( ω − E ) − δ ( ω + 2 E )) + Real Part= 18 v F d ( i ∆ ω + ω ) = πν i ∆ ω + 1 π ω ) (S22)where the real part is obtained through Kramer-Kronig relation by noting that the iω behavior of the imaginary parthas an UV cutoff ∼ ∆ beyond which it approaches a constant.The second term in Eq. (S18) comes from expanding Eq. (S12). For example for the 2D stack, the static free energyas a function of θ s is F (∆ , ∆ p , θ s ) = − νπ (cid:90) π/ − π/ dk z | ∆ sin( k z + θ s ) + i ∆ p cos( k z ) | ln 2Λ | ∆ sin( k z + θ s ) + i ∆ p cos( k z ) | (S23)where k z means the z direction momentum times the interlayer thickness. It can be expanded to O ( θ s ) as c θ s = ν π (cid:90) π/ − π/ dk z (cid:20)(cid:18) − ln 4Λ | ∆( k ) | + 1 (cid:19) (cid:0) − ∆ p cos(2 k z ) θ s (cid:1) + 12 1 | ∆( k ) | p sin ( k z ) cos ( k z ) θ s (cid:21) = ν π θ s (cid:90) π/ − π/ dk z (cid:34) ∆ p ln 4Λ | ∆( k ) | cos(2 k z ) + ∆ p ∆ sin ( k z ) + ∆ p cos ( k z ) sin (2 k z ) (cid:35) = ν (cid:32) − ∆ p ∆ ∆ p ∆ ∆ p + ∆ p ∆(∆ + ∆ p ) (cid:33) θ s O (∆ p ) −−−−→ ν p θ s (S24)The third term is due to the fact that θ s fluctuations are accompanied by x direction polarization. In the BCS limit,each k x chain at k z in momentum space contributes a polarization density of P k z = P D π ( π/ − ϕ k z ) in the 2D stack,and each k x − k y surface at k z contributes P k z = P D [1 − tan( ϕ k z / ϕ k z = ArcTan ∆ p cos k z ∆( k z ) .Noting that ∆( k z ) = ∆ sin( k z + θ ), the change of polarization is related to θ as χ = ∂ θ P x = P D π (cid:90) π/ dk z π ∆ cos( k z )∆ p cos( k z ) sin ϕ = P D π ∆∆ p (cid:90) π/ dk z
11 + ∆ ∆ p tan ( k z ) = P D π ∆ / ∆ p / ∆ p (S25)for the 2D stack and as χ = ∂ θ P x = 14 P D (cid:90) π/ dk z π ∆ cos( k z )∆ p cos( k z ) sin ϕ cos ( ϕ/
2) = P D π ∆∆ p (cid:90) π/ dk z sin ϕ ϕ = P D π ∆∆ p (cid:90) π/ dk z cos k z / ( ∆ ∆ p sin k z + cos k z )1 + ∆∆ p sin k z / (cid:113) ∆ ∆ p sin k z + cos k z = P D π f (cid:18) ∆ p ∆ (cid:19) (S26)for the 3D stack where f ( x ) → x → f ( x ) → π/ (2 x ) as x → ∞ . Note that χ reduces to σ h in the limit of∆ p (cid:28) ∆ which we assume in Eq. (10) of the main text for notational simplicity. In general, χ (cid:54) = σ h because θ s doesnot act completely in the same way as A z since it does not enter the ∆ p term in the Hamiltonian.Finaly, the S A x term is determined by the bare optical conductivity along x . In the 3D stack, the bare opticalconductivity along x direction in the low frequency regime is controlled by the Weyl nodes with the Hamiltonian forone of them being H W = v F k y σ + v x k x σ + v z k z σ where v x = ∆ p /k F and v z = ∆ d . The optical conductivity inthe regime ω (cid:28) ∆ p reads σ x = v x v z d iω c D ( ω ) = 13 νk F ∆ p ∆ (cid:18) | ω | − iπ ω ln ∆ p | ω | (cid:19) . (S27)In the 2D stack of chains, the optical conductivity is controlled by the massive Dirac nodes with Hamiltonian H D = v F k x σ + v z k z σ + ∆ p σ . The optical conductivity reads σ x = v F v z d iω c D ( ω ) = v F ∆ ν (cid:40) − | ∆ || ∆ p | iω ω (cid:28) ∆ pπ (∆ − iπ ω ) ∆ p (cid:28) ω (cid:28) ∆ . (S28) C. Optical conductivity and hyperbolic phase polaritons
Integrating out the phase in Eq. (S18), one obtains the EM Lagrangian: L EM ∼ − c − c (cid:0) c c d A z + χ ( ∂ t A x ) /d + 2 c χA z ∂ t A x (cid:1) + σ h A x ∂ t A z + S A x (S29)in the gauge φ s = 0, rendering the optical conductivity Eq. (9) of the main text. Note that for the 2D stack, theconductivity simplifies to σ = 1 ω − ω (cid:18) χ iω/c + ( ω − ω ) σ x − ω χω χ iωc (cid:19) + σ h (cid:15) xz (S30)To compare, a layered superconductor with interlayer Josephson tunneling has the Lagrangian L ∼ − ν ( ∂ t ϕ s + φ s ) + v z ( ∂ z ϕ s − A z ) along z direction which corresponds to superconductivity: j z ∼ v z A z .Note that in the BEC case and in the limit of ∆ p →
0, one has χ → σ h = 0. Thus there is no Hall response even atnonzero frequency, consistent with the fact that time reversal symmetry is effectively unbroken if ∆ p = 0 (tunnelingbetween the layers is not possible). In general, ∆ p (cid:54) = 0 and there is nonzero AC Hall response.Being optically active, the bulk phase mode hybridizes with photons to form hyperbolic phase polaritons whosedispersion is determined by 1 + πiω σ ij q i q j /q = 0 for the 3D stack [55, 56] and by 1 + πiω σ ij q i q j / | q | = 0 for the 2Dstack. Due to the hall response, there are also chiral surface phase polaritons on the side surfaces/edges circulatingthe x − z plane which we leave for future study. III. STACKS WITH s -TYPE TUNNELING In this section we assume s -type intrinsic hybridization t p = t > t . The ground state order parameter is a real ∆ with the same sign as t which reads ∆( k ) = ∆ cos( k z ) inmomentum space, meaning ( θ a , θ s ) = (0 , E = ± (cid:112) ξ k + 4(∆ + t ) cos ( k z ) withDirac nodes/lines at ( k ⊥ , k z ) = ( k F , ± π/ (2 d )) in the BCS case. We focus on the physics related to Josephson effect,meaning we neglect spatial fluctuations. In the continuous limit, the relevant degrees of freedom are the θ s and theuniform A z . The low energy Lagrangian valid for ω (cid:28) ∆ is L s ∼ ν ∆ (cid:2) ( θ s + dA z ) ∂ t ( θ s + dA z ) + ω θ s (cid:3) . (S31)where ω = t ln Λ∆ . Thus the phase mode is overdamped due to the gapless quasi particle dispersion. Integrating outthe phase θ s , one obtains the optical conductivity along z σ z = ν ∆ ω d − iω + ω (S32)which has a Drude form with width ω . This is a surprising result since simple Dirac nodes/lines do not give such aDrude conductivity at zero temperature since there are no quasi particles. The Drude behavior is the result of couplingto the collective phase mode θ s . Note that the DC conductivity σ DC = ν ∆ d does not depend on t although therecannot be interlayer tunneling without t . The reason is that as electric field increases beyond E c ∼ ω /d ∼ t/d ln Λ∆ ,the system won’t have a fixed θ s but enters the AC Josephson effect regime ˙ θ s = dE z and the current oscillates withfrequency dE z . Therefore, as t increases from zero, it expands the field regime for linear DC response while the linearresponse DC conductivity is a constant value. IV. THREE CHAIN MODEL FOR TA NISE Ta NiSe can be viewed as a 3D stack of its basic element, a composite Ta-Ni-Ta chain [22, 26] shown in Fig. S3.Its excitonic physics can be reasonably captured by this 1D chain, with the excitonic order as bond variables ∆ and∆ . The Lagrangian within the three band model is L = (cid:80) k ψ † k ( ∂ τ + H k ) ψ k + g ( | ∆ | + | ∆ | ) with the mean fieldHamiltonian H k = ξ c ( k ) φ ( k ) 0 φ ∗ ( k ) ξ v ( k ) φ ∗ ( k )0 φ ( k ) ξ c ( k ) , φ ( k ) = e iA ( it sin k + ∆ cos k ) φ ( k ) = e − iA ( it sin k + ∆ cos k ) . (S33) te !" + Δ e !(% ! &") te (!" + Δ ) e !(% " (") 𝐸 NiTaTa t−t t−t
FIG. S3. Schematic of the three chain model for Ta NiSe . There is a mirror symmetry σ ⊥ of reversing the direction of thechain. The Ni chain with active d orbitals (odd under σ ⊥ , giving one valence band) is sandwiched by two Ta chains with actived orbitals (even under σ ⊥ , giving two degenerate conduction bands). Here and following we suppress the half lattice constants a for notational simplicity. We also set ξ c = − ξ v = ξ withoutchanging the qualitative physics. The resulting dispersions for three quasi particle bands are E k = (cid:18) ξ k , ± (cid:113) ξ k + | φ ( k ) | + | φ ( k ) | (cid:19) . (S34)At t = 0, as a result of the mirror symmetry of interchanging the two Ta chains, Eq. (S33) has an U (2) symmetrycorresponding to (∆ , ∆ ) → (∆ , ∆ ) ˆ U where ˆ U is an arbitrary unitary matrix. The intrinsic tunneling t reducesthis symmetry, rendering the low energy manifold to satisfy | ∆ | = | ∆ | . Fixing ∆ = ∆ e iθ , ∆ = ∆ e iθ , the lowenergy Lagrangian in terms of the two phase angles θ and θ is L = K [ θ + A, θ − A ] + c tns (sin θ + sin θ ) + F (S35)where K is the kinetic energy term that vanishes in the static limit, F is the ground state free energy at t = 0 and wehave expanded to O ( t ) terms. The coefficient c tns is c tns = c tns1 t ∆ in the BEC regime. This term can be computedby noting that the condensation energy is δF = (cid:88) k ( | ξ k | − | E k | ) = (cid:88) k (cid:20) − | φ ( k ) | + | φ ( k ) | | ξ k | + ( | φ ( k ) | + | φ ( k ) | ) | ξ k | + O ( | φ | ) (cid:21) , | φ i ( k ) | = t sin k + ∆ cos k + t ∆ sin θ i sin 2 k . (S36)Considering the summation over k , the leading θ i dependence comes from the O ( φ ) terms and reads F = t ∆ (sin θ + sin θ ) (cid:88) k sin k | ξ k | ≡ c tns1 t ∆ (sin θ + sin θ ) . (S37)Therefore, the ground states lie at the lines satisfying θ = − θ or θ = θ + π in the ( θ , θ ) space. From the equationof motion implied by Eq. (S35), the z direction current is j z = − ∂ A L = ( ∂ θ − ∂ θ ) F = 4 c tns sin (cid:0) θ + θ (cid:1) sin ( θ − θ ).Mean field analysis considering effects beyond the three chain model has found that the actual symmetry is Z , andthere are only two degenerate ground states lying at ( θ , θ ) = (0 , π ) and ( θ , θ ) = ( π,
0) [22].With a fast electric field E = − ∂ t A along z (perpendicular to the chain), as shown in Fig. S3, only the kinetic termin Eq. (S35) matters in the equation of motion, giving ∂ t θ = E = − ∂ t θ . (S38)In other words, the phase θ i has to adjust to cancel the Peierls phase A , and the order parameter won’t catch upwith the dynamics in our choice of gauge. Therefore, an electric field pulse can rotate θ counterclockwise from 0to π while rotating θ clockwise from π to 0, switching the two ground states. In the mean time, there is nonzero z direction current oscillating between the three chains. The bulk inversion ˆ P is unbroken so there is no spontaneouspolarization to identify the change of state, but the z direction Josephson current could lead to measurable radiations. V. EFFECT OF SPIN
In this section, we show that the main conclusion of this paper is not affected by spin degrees of free-dom of the electron and hole bands which can be represented by the fermion creation operators ψ † =( ψ † c ↑ , ψ † c ↓ , ψ † v ↑ , ψ † v ↓ ). After decomposing the interband part of the density-density interaction V = (cid:82) drdr (cid:48) V ( r − r (cid:48) ) (cid:80) s ,s ψ † c,s ( r ) ψ c,s ( r ) ψ † v,s ( r (cid:48) ) ψ v,s ( r (cid:48) ) in the s -wave electron-hole pairing channel (the dominant channel), oneobtains the spinful Lagrangian L = ψ † ( ∂ τ + H k ) ψ + 1 g Tr[ ˆ∆ † ˆ∆] + 1 g X (cid:18) X + 1 ω ˙ X (cid:19) (S39)where the mean field Hamiltonian is H k = (cid:18) ξ c ˆ I ˆ∆( k ) + i ∆ p f k ˆ I + X ˆ I ˆ∆ † ( k ) − i ∆ p f k ˆ I + X ˆ I ξ v ˆ I (cid:19) , ˆ∆ = ∆ ˆ I + −→ ∆ · −→ σ . (S40)Here we have added the shear phonon X which does not flip the spins. The Lagrangian Eq. (S39) has an SU(2)symmetry of spin rotations, whose element is represented as ˆ U g : ( ψ †↑ , ψ †↓ ) → ( ψ †↑ , ψ †↓ ) ˆ U , ˆ∆ → ˆ U † ˆ∆ ˆ U = ∆ ˆ I +( ˆ R −→ ∆) ·−→ σ where ˆ U is a 2 × R is a 3 × componentof the order parameter is invariant under spin rotation, and thus corresponds to a spin singlet condensate. The −→ ∆ transforms as a 3-vector and corresponds to a triplet condensate. The singlet and triplet form two irreduciblerepresentations of the symmetry. The phonon appears in the ˆ I channel and couples linearly only to the singletcondensate [4, 20].Without the phonons, i.e., at g X →
0, the singlet and triplet states are degenerate in energy. The pure singlet stateis simply two copies of the ferroelectric states (or time reversal broken states in the stacks) discussed in the main text,both having identical orbital properties such as the electrical polarization. Therefore, all conclusions are the same asspinless models discussed in the main text. The pure triplet state is formed by two copies of the ferroelectric stateswith opposite electrical polarization, corresponding to a spin polarized state. For example, the −→ ∆ = (0 , , ∆) statehas nonzero (cid:104) s z (cid:105) on one side of the sample and opposite (cid:104) s z (cid:105) on the other side.In the bilayer system in Fig. 1 of the main text, the phonon X couples linearly and cooperates withthe singlet condensate, lowering its energy than the triplet state by some energy δF which is about δF ∼ Λ (cid:16) g + g X e − ν ( g + gX ) − g e − νg (cid:17) in the BCS weak coupling case. Therefore, as long as the shear phonon mode ofthe bilayer is not infinitely rigid, the ground state is the singlet state which is also ferroelectric.In the 3D stacks of alternating electron-hole layers and 2D stacks of electron-hole chains, the shear phonon does notcouple linearly with the grounds state order parameters, and thus does not lift the degeneracy of singlet and tripletstates. Among the effects not considered in this paper, direct exchange interaction favors the triplet state while kineticexchange favors the singlet [20]. The singlet state is made of two identical copies of time reversal breaking ones withspin up and down, respectively. The triplet state has opposite time reversal breaking for spin up and down, exhibitingmicroscopically circulating spin currents instead of charge currents [4]. As a result, the 2D stack is in a quantum spinhall state, the 3D stack is a dirac semimetal to which spin orbit coupling might open a gap, and the anomalous hallresponse vanishes due to time reversal symmetry. However, the order parameter steering and Josephson current isorbital physics which works universally for both singlet and triplet states. VI. CORRELATION FUNCTIONS
This section defines the two point correlation functions that are coefficients of the quadratic terms in the effectiveaction after integrating out the fermions. The correlation function χ σ i σ j is defined as χ σ i σ j ( q ) = (cid:68) ˆ T (cid:0) ψ † σ i ψ (cid:1) ( r,t ) (cid:0) ψ † σ j ψ (cid:1) (cid:69) (cid:12)(cid:12)(cid:12)(cid:12) q = (cid:88) ω n ,k Tr [ G ( k, iω n ) σ i G ( k + q, i ( ω n + Ω)) σ j ] . (S41)where ˆ T is the time order symbol, Tr is over the spinor basis, q = ( q , i Ω) and G ( k, iω n ) = (cid:68) ˆ T ψ ( x ) ψ † (0) (cid:69) (cid:12)(cid:12)(cid:12)(cid:12) k,iω n = 1 iω n − H k (S42)is the electron Green’s function. At zero temperature, rotating i Ω to ω , Eq. (S41) reads χ σ i σ j ( ω, q ) = 12 (cid:88) k ω − ( E + E (cid:48) ) (cid:40) ( E + E (cid:48) )Tr (cid:20) σ i σ j − H k σ i H k (cid:48) σ j EE (cid:48) (cid:21) + ω Tr (cid:20) σ i H k (cid:48) σ j E (cid:48) − H k σ i σ j E (cid:21) (cid:41) (S43)where H kk