Exciton condensation in bilayer spin-orbit insulator

We investigate the nature of the magnetic excitations of a bilayer single-orbital Hubbard model in the intermediate-coupling regime. This model exhibits a quantum phase transition (QPT) between a paramagnetic (PM) and an insulating antiferromagnetic (AFM) phase at a critical value of the coupling strength. By using the random phase approximation, we show that the QPT is continuous when the PM state is a band insulator and that the corresponding quantum critical point (QCP) arises from the condensation of preformed excitons. These low-energy excitons reemerge on the other side of the QCP as the transverse and longitudinal modes of the AFM state. In particular, the longitudinal mode remains sharp for the model parameters relevant to Sr$_{3}$Ir$_{2}$O$_{7}$ because of the strong easy-axis anisotropy of this material.


I. INTRODUCTION
Low-dimensional Mott insulators play a central role in correlated electron physics because of the novel states of matter and phase transitions that they can host. In the strongcoupling limit, U/|t| 1, the on-site repulsion U "freezes out" the charge degrees of freedom, turning the Mott insulator into a quantum magnet described by an effective spin Hamiltonian. The lattice connectivity can then be exploited to generate competing exchange interactions that induce quantum phase transitions (QPTs) between different states of matter. One of the simplest examples is provided by bilayer materials [1][2][3][4][5][6], where the competition between intralayer and interlayer hopping amplitudes, t and t z , can induce a transition between an antiferromagnetic (AFM) state and a quantum paramagnet (QPM) comprising local singlet states on the interlayer dimers. This QPT, which is typically induced by applying pressure, has been extensively studied in multiple quantum magnets to understand different properties of the QCP, such as the emergence and decay of the longitudinal mode (LM) that is present in the AFM state [7][8][9][10][11][12][13][14][15][16].
The discovery of low-dimensional intermediate-coupling 4d and 5d-electron correlated insulators introduces a knob, U/t, that can be used to unleash the charge degrees of freedom via pressure or strain. For instance, the iridate materials [17][18][19][20][21][22] have a charge gap ∆ that is comparable to the magnon bandwidth W . The reduction of U/t leads to the suppression of the AFM ordering in favor of a paramagnetic state. In bilayer materials, such as Sr 3 Ir 2 O 7 , the paramagnetic state can either be metallic or a band insulator depending on the spin orbit coupling. Similarly to the large U/|t| limit, where the transition from the QPM to the AFM state can be described as a triplon condensation, the QPT from the band insulator to the AFM state corresponds to condensation of preformed magnetic excitons. These bound states reemerge on the other side of the QCP as longitudinal and transverse modes (magnons) of the AFM state. In contrast, the longitudinal mode (LM) is absent in the AFM phase induced by a metal-insulator transition with perfect Fermi surface (FS) nesting.
In this article, we study the possible QPTs induced by re-ducing the U/|t| ratio in bilayer materials. Although the transition into the band insulator has some similarities with the transition into the QPM of pure spin systems (U/|t| 1) [7][8][9][10][11][12][13][14][15][16], there are also some important differences associated with the enhanced charge fluctuations. By applying our analysis to the easy-axis bilayer antiferromagnet Sr 3 Ir 2 O 7 , with ∆ = 130meV [23] and W = 70meV [24], we reveal the existence of a LM in some regions of the Brillouin zone, which arises from the band-insulating character of the noninteracting limit of the model: The Néel phase is induced by condensation of preformed S z = 0 excitons at a critical coupling strength U = U c . This exciton reemerges in the Néel phase (U > U c ) as a LM whose energy scale and dispersion are consistent with previous resonant inelastic x-ray scattering (RIXS) measurements [24][25][26][27]. It is important to note that the LM that emerges near the same QCP in the strong coupling limit of the Hubbard model [25] exists over the whole Brillouin zone and its energy is much lower than the charge gap because it is obtained from a pure spin model. In contrast, the LM that we are proposing for Sr 3 Ir 2 O 7 only exists in finite regions of momentum space, around the wave vectors q = 0 and q = (π, π), because it is induced by strong charge fluctuations (U ∼ U c ). Indeed, the mode disappears inside the particle-hole continuum for wave vectors that are far enough from q = 0 and q = (π, π).
Our results then suggest that Sr 3 Ir 2 O 7 is a realization of the long sought excitonic insulator that was predicted almost 60 years ago [28][29][30]. However, higher resolution RIXS experiments are needed to confirm this prediction.

II. MODEL
We consider a bilayer single-orbital Hubbard model H = −H K + H U , with H U = U r n r↑ n r↓ and where t, t z , α, U ∈ R, c † r ≡ [c † ↑,r , c † ↓,r ], is the Nambu spinor of the electron field [r ≡ (r ⊥ , l), l = 1, 2 is the layer index and r ⊥ = r 1 a 1 + r 2 a 2 ], and a ν (ν = 1, 2) are the primitive vectors of the square lattice of each layer. The sign r takes the values 1 (−1) for r ∈ A (B) sublattice of the bipartite bilayer system. This Hamiltonian is an effective model for a bilayer system with finite SOC, which is realized in ruthenates and iridates. For example, the large SOC of the bilayer iridates splits the 5d t 2g orbitals of the Ir 4+ ion into J = 1/2 and 3/2 multiplets [31]. Consequently, H becomes a lowenergy model for the 5d hole (5d 5 electronic configuration) of the strontium iridates after projecting the relevant multi orbital Hubbard model onto the lowest energy J = 1/2 doublet. The phase α arises from hopping matrix elements between d xz and d yz orbitals allowed by staggered octahedral rotations. The phase of the intralayer hopping is gauged away by applying a sublattice-dependent gauge transformation.

III. MEAN FIELD APPROXIMATION
For α = 0, the model (1) has an easy z-axis spin anisotropy and the ground state of the AFM phase can have Néel ordering, S µ r = (−1) γr M δ µz , where γ r = (1 + r )/2, S µ r = 1/2c † r σ µ c r (µ = x, y, z), and M is the magnetization. A mean field decoupling of H U leads to r c r is the electron density and the second term can be absorbed into the chemical potential.
By Fourier transforming annihilation and creation operators, where N u = 1 2 N s is the number of unit cells, and k ≡ (k ⊥ , k z ) runs over the first Brillouin zone (BZ): k ⊥ = k 1 b 1 + k 2 b 2 , with b 1 = (1/2, −1/2), b 2 = (1/2, 1/2), k 1 , k 2 ∈ [0, 2π) and k z = 0, π, we obtain the momentum space representation of with H MF k is diagonalized by a unitary 4 × 4 matrix U (k), where n ≡ (s, σ) with s = ±, σ =↑, ↓ and each column of U (k) is an eigenvector of H MF k : with and δ = U M . The corresponding eigenenergies, ε sσ (k) = s δ 2 +|b σk | 2 , are independent of the spin flavor, |b σk | 2 = b 2 k and z σk ≡ z k , because of the U(1) invariance of H under global spin rotations about the z axis. The noninteracting system is then metallic when the band gap at each k point, Because we are considering a bilayer system, the half-filled condition corresponds to an integer number of electrons per unit cell, implying that the ground state of the noninteracting (U = 0) system can either be a metal or a band insulator, as illustrated in Fig. 1(a). For α = 0, the metal to insulator transition occurs via a semimetallic state with small and identical electron and hole pockets that shrink into a quadratic Fermi point. For instance, for t z = 0, the FS is the square defined by the equations k 1 + k 2 = ±π and k 1 − k 2 = ±π. Upon increasing |t z |, while keeping α = 0, the square shrinks into a circular pocket that finally disappears for |t z | > 4|t|. In contrast, the noninteracting system is always a band insulator for α = 0.
The on-site repulsion U induces AFM ordering via a QPT that depends on the nature of the noninteracting state. In both cases, noninteracting metal and insulator, the longitudinal AFM susceptibility (q = 0) at ω = 0 is given by where γ 1,2 = A, B are sublattice indices. Note that the antiferromagnetic susceptibility can diverge at q = 0 because the nonmagnetic unit cell contains one site of each sublattice. From the mean field calculation, the critical interaction strength is given by The resultant phase boundaries for several values of α are shown in Fig. 1(b). For α = 0 and U = 0, the system is metallic for |t z /t| < 4. In this case, the perfect nesting due to the coincidence of the particle and hole Fermi surfaces leads to a logarithmic divergence of χ (γ1,z);(γ2,z) 0 (q, ω = 0) at q = 0. Note that the divergence becomes ln 2 q for a square FS (t z = 0) because of (from top to bottom). For α = 0 and U = 0, the system is a metal for |tz/t| < 4 and a band insulator for |tz/t| > 4. The metal becomes an AFM insulator for an infinitesimally small U . For α = 0, the present model (1) exhibits a continuous phase transition from the band insulator to the AFM insulator at a finite U/t. The solid circle symbol indicates the estimated parameter set for the bilayer iridate Sr3Ir2O7 with α = 1.45. (c) Energy diagram as a function of U for the parameter set relevant to Sr3Ir2O7. The red dashed line indicates the energy of the S z = ±1 exciton that becomes the transverse mode for U > Uc. The blue dotted line indicates the energy of the S z = 0 exciton that becomes a longitudinal mode for U > Uc. The gray area indicates the particle-hole continuum. From a comparison with RIXS data [24][25][26][27], we estimate that U = 0.33 eV for Sr3Ir2O7.
the Van Hove singularity in the density of states at the Fermi level (the corners of the square correspond to saddle points of the dispersion relation). In other words, the metal becomes an AFM insulator for an infinitesimally small value of U , implying that the critical interaction strength U c is equal to zero. In contrast, the system becomes a band insulator (finite charge gap) for α = 0 because the sublattice symmetry is no longer present. As shown in Fig. 1(b), the metal-insulator transition at U = 0 is then replaced by a continuous QPT between the band and the AFM insulators at a finite U value, U = U c , because the integral (10) is now convergent. This phase diagram clearly shows the significance of the spin-orbit coupling: The proximity to the critical point is caused by a finite α even for small t z /t. The physical interpretation of the difference between the metallic and the band-insulating cases will become clearer upon analyzing the behavior of the magnetic excitations.
The order parameter is determined by solving the selfconsistent mean field equation Near the critical point U = U c , we obtain M ∝ (t/U )e −1/(ρ0U ) for the metal, where ρ 0 the density of states at the Fermi level, and

IV. EXCITON CONDENSATION
In view of the U (1) invariance of H and the AFM ground state, the transverse and longitudinal spin fluctuations are decoupled from each other. Within the random phase approximation (RPA), the magnetic susceptibilities of the transverse and the longitudinal modes are given by respectively, where τ 0 is the 2 × 2 identity matrix. Here χ +− and χ zz refer to 2 × 2 matrices in the sublattice space, while χ +− 0 and χ zz 0 are the bare magnetic susceptibilities. See Appendix A for details of the RPA calculation.
The eigenfrequencies ω q of the collective transverse modes (magnons) can be extracted from the poles of the transverse RPA susceptibility: det[τ 0 − U χ +− 0 (q, ω q )] = 0. The spectrum is fully gapped because the U(1) symmetry of H is not spontaneously broken by the AFM ordering along the z-axis. As shown in Fig. 1(c), the transverse modes remain gapped at U = U c and they become the S z = ±1 exciton modes of nonmagnetic band insulator for U < U c . Note that this gap closes in the absence of spin-orbit coupling (α = 0) because the Hamiltonian becomes SU(2) invariant.
The most interesting feature is the emergence of a LM below the particle-hole continuum around the Γ point of the second Brillouin zone [q = (π, π, π)]. The origin of this mode can be understood by analyzing the excitation spectrum of the band insulator for U < U c . The bare magnetic susceptibility at q = 0, which is equivalent to (π, π, π), is where τ x is the Pauli matrix and The function Φ(ω) is real, and it increases monotonically with ω up to the lower edge of the particle-hole continuum ω ph = 2min(b k ), where it diverges: Φ(ω = 0) = 1/U c and lim ω→ω ph Φ(ω) = ∞. The putative pole of the longitudinal susceptibility is determined by the condition which implies that a pole must exist in the energy window 0 ≤ ω ≤ ω ph for 0 < U < U c . As long as U c is finite, i.e., the noninteracting system is a band insulator, the pole appears just below the gap ∆ = ω ph that signals the onset of the particle-hole continuum. This pole corresponds to the formation of an S z = 0 exciton. Upon examining the transverse susceptibility, we also find a doubly degenerate pole associated with the formation of S z = ±1 excitons [see Fig. 1(c)], whose energy is higher than the energy of the S z = 0 exciton because of the easy-axis anisotropy. We note that the three excitonic states become degenerate in the absence of spin-orbit interaction because H is SU(2) invariant in that limit. The binding energy exhibits the singular behavior E b ∝ ω ph exp(−κω ph /U ), characteristic of 2D systems in the weak-coupling limit (κ is a nonuniversal number). As shown in Fig. 1(c), the S z = 0 exciton becomes soft at The condensation of this mode signals the onset of the Néel phase with magnetic moments pointing along the z axis due to the effective easy-axis anisotropy.
In the AFM phase (U > U c ), the bare magnetic susceptibility at q = 0 is given by If the noninteracting state is a band insulator (finite U c ), it holds that b k > 0 for any k, implying that a pole must then exist in the window 0 ≤ ω ≤ ω ph because Φ(ω = 0) < 1/U and Φ(ω ph ) > 1/U . For U U c , the exciton energy scales as ω L = 4 (U − U c )/(DU 2 c ). Upon further increasing U , ω L increases quickly and approaches the lower edge of the particle-hole continuum asymptotically. This behavior can be understood by considering the large-U limit. In this limit, the effective particle-hole interaction that provides the "glue" for the formation of the LM is the exchange interaction, on the order of t 2 z /U , along the vertical bonds, implying that the binding energy of the LM must vanish for U → ∞. The particle and the hole break only one AFM link when they occupy the two sites of the same vertical bond, while they break two vertical links when they occupy two different vertical bonds. In contrast, the binding energy of the low-energy transverse modes is on the order of U because the particle and the hole occupy the same site. See Appendix B for the analysis of the exciton wave function.
The situation is qualitatively different for a metallic noninteracting system with perfect FS nesting. If min(b k ) = 0, the condition (15) cannot be fulfilled because Φ(ω < ω ph ) ≤ Φ(ω ph ) = 1/U . Therefore, the LM of the AFM phase is absent in this case.

V. LONGITUDINAL MODE OF BILAYER IRIDATE
From the above analysis, we predict that the bilayer iridate Sr 3 Ir 2 O 7 should exhibit a LM in a relatively small region around the center of the Brillouin zone. This predic-tion is qualitatively different from a previous interpretation of the RIXS data [25], based on a pure spin model (i.e. strong-coupling limit), that reports the existence of a LM over the whole Brillouin zone. The effective Hamiltonian (1) for Sr 3 Ir 2 O 7 can be obtained by constructing a tight-binding model of the t 2g orbitals and projecting it onto J eff = 1/2 lowest energy doublet [32]. The resulting parameters of the effective single-orbital Hubbard model can be optimized to reproduce the experimentally observed magnon dispersion [24][25][26][27]: α = 1.45 and t = 0.11, t z = 0.09, U = 0.33 in units of eV. This value of α ≈ π/2 produces a strong easy-axis anisotropy, realized in a tetragonal elongation of octahedra consistent with |t| > |t z |, as estimated for Sr 3 Ir 2 O 7 [32].
The dynamical spin structure factor where S µ q = 1 √ N r S µ r e iq·r , is obtained from the dynamical spin susceptibility given in Eqs. (12) and (13). Figure 2 shows the out-of-phase (q z = π) transverse (OT) response, S xx (q, ω) = S yy (q, ω), and the out-of-phase longitudinal (OL) response, S zz (q, ω), for U = 0.33, U c ≈ 0.28, and 0.23 eV. The in-phase, or q z = 0, transverse (IT) response is not shown because it is practically identical to the OT response [24]. This is a direct consequence of the strong easyaxis effective interlayer exchange that suppresses the single magnon tunneling between the two layers. In the large-U limit, the effective interlayer exchange includes only Ising and Dzyaloshinskii-Moriya exchange interactions that do not split the two degenerate single-layer modes. This anisotropy is also responsible for the large spin gap revealed by RIXS measurements [24]. No excitation peak is found below the particle-hole continuum in the in-phase longitudinal response. To match the RIXS peak at ( π 2 , π 2 ), we added the intralayer next nearest neighbor hopping 0.012 eV ≈ t/10, which only changes the eigenenergy in the above argument. The overall dispersion curve for U = 0.33 eV, shown in Figs. 2(a) and 2(b), is consistent with the RIXS measurements [24][25][26][27]. The resultant charge gap (≈ 130 meV) is also consistent with the experimental observation [23]. The structure of the continuum reflects the underlying electron bands. Unfortunately, the resolution of the reported RIXS data [24][25][26][27] is not enough to extract structures in the continuum. It would be of interest to compare the structure and the onset of the continuum to our calculation.
The most salient feature of the results shown in Fig. 2 is the sharp OL mode near q ⊥ = (0, 0) and (π, π). Interestingly, the sharp OL mode appears only at restricted wave vectors because the particle-hole continuum exists in the same energy scale. The OL mode becomes gapless at U = U c , while the OT mode remains gapped, as shown in Figs. 2(c) and 2(d); the exciton peaks in the band insulator for U < U c remain sharp, as shown in Figs. 2(e) and 2(f). The U dependence of the energy of this OL mode has important consequences for its stability. While the kinematic constraints do not allow for its decay into two transverse modes for the parameters of Sr 3 Ir 2 O 7 , the fact that ω L (ω T ) is of order U (t 2 /U ) in the large-U/t limit implies that the decay becomes kinematically and (e) S xx (q, ω) = S yy (q, ω) and (f) S zz (q, ω) for U = 0.23 eV. The hopping parameters were chosen to reproduce the RIXS data of the bilayer iridate Sr3Ir2O7 [24][25][26][27]: t = 0.11 eV, tz = 0.09 eV, α = 1.45, and the intralayer next nearest neighbor hopping 0.012 eV. The broadening factor is η = 10 −4 eV. Intensity larger than the maximum value in the scale (color) bar is plotted in the same color. allowed above a certain value of U . In other words, the sharp OL mode that we predict for Sr 3 Ir 2 O 7 is protected by the large spin gap generated by the easy-axis anisotropy and by the relatively small value of U/t. Thus, the sharpness of the OL peak and the fact that ω L and ω T are comparable energy scales are clear indicators of the proximity of Sr 3 Ir 2 O 7 to the QCP at U = U c . Because the real material has a small inter-bilayer hopping, the dimension of the effective theory that describes this QCP, D = 3 + 1, coincides with the upper critical dimension (Gaussian fixed point). This means that the mean field theory adopted here is correct up to logarithmic corrections.
Finally, we note that higher order processes through the J eff = 3/2 orbitals induce further neighbor hopping terms that can make the noninteracting system semimetallic without the FS nesting. In this case, a first-order metal-to-insulator transition occurs at a finite value of U [32], implying that there are two alternative scenarios for suppressing the AFM order via reduction of the coupling strength (U/t) in bilayer iridates. For the case of Sr 3 Ir 2 O 7 , we predict that, if the material transitions into a band insulator, a QCP must exist at U c ≈ 0.28 eV, which is only 15% smaller than the estimated value (U = 0.33 eV) at ambient pressure. The coupling strength U/t can be reduced by applying high pressure [33]. In this scenario, the LM becomes soft at U = U c . In contrast, if the transition is of first order, the material becomes metallic without the softening of the LM. While absent in a metallic phase, the characteristic S z = 0 exciton peak appears in the bilayer spin-orbit band insulator for U < U c , as shown in Fig. 2(f).

VI. SUMMARY AND DISCUSSION
We have reported the existence of a sharp LM in antiferromagnetic bilayer Mott insulators with large spin-orbit coupling. Whenever the noninteracting system is a band insulator, S z = 0, ±1 excitons emerge from the particle-hole continuum for an infinitesimally small value of U . For materials with easy-axis anisotropy, such as the bilayer iridate Sr 3 Ir 2 O 7 , the Ising-like AFM ordering results from the condensation of S z = 0 excitons at a critical value U = U c . If the noninteracting system is a metal with FS nesting, the formation and condensation of the particle-hole pairs that produce the magnetic moments occur simultaneously at U c = 0 + . Similarly to the case of single-layer Mott insulators, such as Sr 2 IrO 4 , the LM is absent in the resulting AFM state.
We emphasize that the LM that we discussed in this work exists as a sharp mode only close enough to the QCP between the AFM phase and the paramagnetic band insulator. A similar QCP still exists in the large-U limit, where it divides the AFM state from a quantum paramagnetic spin state, whose mean field description is a product state of rung singlets [25], which is adiabatically connected with the band insulator. Note that such a mean field state can only be captured by expanding the variational space of product single-site sates that is assumed by the conventional AFM mean field approach used in this work. However, our estimate of the hopping amplitudes for Sr 3 Ir 2 O 7 , consistent with the effective model obtained from the three-orbital model [32], leads to J z /J = (t z /t) 2 ≈ 0.67 in the strong coupling limit. This value of the exchange coupling ratio is significantly smaller than the critical value (J z /J) c . The bond-operator mean field theory gives (J z /J) c = 4 for the isotropic and the easy-axis cases [16,25], while the exact value for the SU(2) invariant case is very close to (J z /J) c = 2.5221 according to quantum Monte Carlo simulations [34]. Therefore, as we discussed in the previous section, the sharp LM mode does not survive in the large-U limit of the Hubbard model of Sr 3 Ir 2 O 7 (the LM loses its sharpness far enough from the QCP because the kinematic constraints allow it to decay into pairs of transverse modes). Moreover, as shown in Fig. 2, the sharp LM appears only at restricted wave vectors because the particle-hole continuum exists in the same energy scale. This characteristic feature of Sr 3 Ir 2 O 7 indicates that this material is an excitonic insulator [28][29][30] near the critical point U = U c . In other words, the proximity of Sr 3 Ir 2 O 7 to quantum criticality is caused by strong charge fluctuations, which are absent in the pure spin model.
Finally, it is worth mentioning that the longitudinal mode can also exist in bilayer Mott insulators with easy-plane anisotropy. The corresponding XY -AFM state results from the condensation of the S z = ±1 excitons, while the S z = 0 exciton remains gapped at U = U c and becomes the LM for U > U c . The main difference is that this LM, also known as "Higgs mode," is critically damped in (3+1)D because it is allowed to decay into two transverse magnons (Goldstone modes) of the AFM state [11-13, 15, 35, 36]. It is of great interest to seek the sharp LM and the Higgs mode in other physical systems. Recently, fermionic systems described by bilayer Hubbard models have been realized in cold atoms [37], implying that the exciton condensation that we have discussed here can also be realized in these systems.
In view of the U (1) invariance of H and of the ground state, the transverse and longitudinal spin fluctuations are decoupled from each other. The same is true for the charge fluctuations. Consequently, we rewrite the (bare) interaction vertex in three different forms which account for the longitudinal and transverse spin fluctuations, respectively, where the interaction vertex takes three equivalent forms, σ 0 is the 2 × 2 identity matrix, and σ ± = 1 2 (σ x ± iσ y ). We are ready to compute the dynamic spin-charge susceptibility within the RPA. Figures 3(a) and 3(b) show the fermion propagator and the bare vertex, respectively. The free fermion Green's function is represented by a solid line, where P nσ (q) = |X nσ (q) X nσ (q)| is the projector on the eigenstate |X nσ (q) of the mean field Hamiltonian. The (bare) dynamic spin-charge susceptibility of the noninteracting mean field Hamiltonian is given by where γ 1 , γ 2 = A, B are sublattice indices and µ, ν = 0, ±, z are charge-spin indices. Here we have introduced the polarization factor F (γ1,µ);(γ2,ν) kσ;lσ where γ takes the values 1 (−1) for γ = A(B), and τ 0 and τ z are the 2 × 2 identity and the Pauli matrix of the sublattice space, respectively. The diagram of χ 0 is shown in Fig. 3(c). Because of the U (1) spin rotation symmetry about the z axis, F (γ1,µ);(γ2,ν) kσ;lσ (k; q) has the following structure: The same structure holds for χ (γ1,µ);(γ2,ν) 0 , confirming that the transverse spin fluctuations, longitudinal spin fluctuations, and charge fluctuations are decoupled from each other. In the following, we focus on the spin channel, which is of main interest for this work. Figure 3(d) represents the magnetic susceptibility at the RPA level. The results are for the transverse channel and for the longitudinal channel. Note that χ +− RPA , χ −+ RPA , and χ zz RPA are 2 × 2 matrices in the sublattice space, while χ +− 0 , χ −+ 0 , and χ zz 0 are the corresponding bare magnetic susceptibilities, respectively.

Appendix B: Exciton wave function
We here investigate the wave function of the exciton formed by multiple particle-hole pairs with S z = 0. In the large-U limit, it is mainly composed of a single particle-hole pair because the energy cost of creating a particle-hole is of order ∼ U . The wave function of the pair is anticipated to be a tightly bound state because the mean field bandwidth of an electron/hole ∼ t 2 /U is comparable to the binding energy ζt 2 z /U , where ζ 0.88 obtained by solving the pole equation for the relevant set of the hopping parameters.
Here iG (2) is the noninteracting particle-hole Green's function where P s,σ (k) =|X s,σ (k) X s,σ (k)| is the projector to the (s, σ) eigenstate at k. In the RPA (see Fig. 4), the T matrix is given by Note that the bare Hubbard interaction is written in terms of the longitudinal component of the spin operator, H int = r [(1/4)(n r ) 2 − (S z r ) 2 ] with n r = σ c † σ (r)c σ (r) and S z r = (1/2) σ σc σ (r) † c σ (r). Under the RPA, the T matrix that describes the renormalized interaction between electrons includes contributions from longitudinal spin fluctuations.
For each combination of sublattices, it reads , where ζ 0.88 determines the exciton binding energy E b ≡ ω ph − ω L ζt 2 z /U . The amplitude ψ ↓↓ (r 1 ∈ A, r 2 ∈ B) determines the spectral weight of the configuration with one hole with spin ↓ at r 1 ∈ A and one electron with spin ↓ at r 2 ∈ B. Since the ordered moment is −M on sublattice A and M on sublattice B, the exciton is created by moving either ↓ spin from sublattice A to B or a ↑ spin from sublattice B to A. On the condition that a hole with spin ↓ is pinned at r 1 = 0 (sublattice A and layer 1), the distribution of the spectral weight, i.e., |ψ ↓↓ (0, r 2 )| 2 , is shown in Fig. 5. In the large-U limit, a hole and an electron occupy different sublattices [ Fig. 5(h)]. The bound state with the particle-hole pair occupying a vertical bond has the largest spectral weight and the probability of finding the particle and the hole on different layers is higher than the probability of finding them on the same layer. As a result of the strong binding energy relative to the bandwidth of the single electron or hole (∼ E b ), the size of the bound state is comparable to one lattice constant. Figure 5 shows the evolution of the real space distribution of the particle-hole pair as a function of U . The ordered magnetic moment, M = | G|S z r |G | decreases upon reducing U/t within the ordered phase [Figs. 5(e)-5(g)] and the probability of finding the particle and the hole on the same site increases. Thus, S z r |G provides a reasonably good approximation of the exciton eigenstate. At U = U c , we find ψ A↓;B↓ (k) = ψ * B↑;A↑ (k) = 0 [ Fig. 5(d)]. The magnetically ordered state emerges as a superposition of the nonmagnetic ground state of the band insulator and states containing multiple coherent excitons. The choice of the phase factor carried by each condensing exciton, equal to 0 or π, corresponds to the Z 2 time-reversal symmetry breaking (M > 0 or < 0 on either sublattice). Note that the magnetic moments on the two sublattices must be opposite because of the relative "−" sign between ψ σ;σ (r 1 ∈ A, r 2 ∈ A) and ψ σ;σ (r 1 ∈ B, r 2 ∈ B), i. e., the system develops Néel magnetic ordering for U > U c .
The time-reversal symmetry is restored as U decreases further below U c (in the band insulator), and the wave function of the exciton becomes extended [see Figs. 5(a)-5(c)], as expected from the reduction of the binding energy. For zero binding energy (ζ = 0), the integral 1 Nu k ψ A↓;B↓ (k)e ik·(r2−r1) that determines the real space wave function ψ ↓↓ (r 1 ∈ A, r 2 ∈ B) is dominated by the singular points k 0 that minimize the particle-hole excitation energy 2 δ 2 + b 2 k . In the simplest case where there is a unique singular point k 0 , ψ ↓↓ (r 1 ∈ A, r 2 ∈ B) ∝ 1 Nu e ik0·(r2−r1) takes the form of a plane wave.