Dynamics of a single exciton in strongly correlated bilayers

We formulated an effective theory for a single interlayer exciton in a bilayer quantum antiferromagnet, in the limit that the holon and doublon are strongly bound onto one interlayer rung by the Coulomb force. Upon using a rung linear spin wave approximation of the bilayer Heisenberg model, we calculated the spectral function of the exciton for a wide range of the interlayer Heisenberg coupling \alpha=J_{\perp}/Jz. In the disordered phase at large \alpha, a coherent quasiparticle peak appears representing free motion of the exciton in a spin singlet background. In the N\'{e}el phase, which applies to more realistic model parameters, a ladder spectrum arises due to Ising confinement of the exciton. The exciton spectrum is visible in measurements of the dielectric function, such as c-axis optical conductivity measurements.


I. INTRODUCTION
An exciton is the bound state of an electron and a hole, and considering their bosonic character the question immediately arises whether they can condense into an exciton Bose condensate 1 . The quest for such exciton superfluidity has, over the past decade, increasingly focussed its attention to layered structures where one layer contains holes and the other layer contains electrons 2 . The Coulomb attraction between the electrons and holes then allows for the formation of so-called interlayer excitons. In 2004 a condensate of interlayer excitons was successfully created in a heterostructure of two 2DEGs under the application of a perpendicular magnetic field 3 . Since then many other candidate materials were suggested that should support interlayer exciton condensation in the absence of magnetic fields, such as graphene [4][5][6][7] or topological insulators 8 . One class of candidate materials has not been considered yet, namely the Mott insulators 9 . The strong interactions between electrons make these materials currently one of the most fascinating and the least understood solid state compounds. When making heterostructures of p and n-doped quasi-two-dimensional CuO 2 layers one expects the formation of interlayer excitons, and these excitons will interact strongly with magnetic excitations, possibly leading to unexpected dynamics. To explore all these unexpected dynamics of the excitons in the strongly correlated system such as the exciton condensation, understanding the dynamics of single exciton will be the first step.
Heterostructures of p and n-doped cuprates can be typically described by a strongly correlated model: the bilayer t−J model, which is extended from two single-band t − J model for each layer with coupling terms between the layers as following: where H t is the hopping of electrons in each layer and H J is the bilayer Heisenberg model describing the undoped Mott insulating state Here c ilσ and s il denotes the electron and spin operators respectively on site i in layer l = 1, 2. The Heisenberg H J is antiferromagnetic with J > 0 and J ⊥ > 0. The last term H V in (1) is the Coulomb attraction between a vacant site (holon) and double-occupied site (doublon) in the same rung, described by which is the force required to form an exciton in the same rung. Without loss of generality, we assume that layer '1' has n-type doping with the constraint σ c † i1σ c i1σ ≥ 1 and layer '2' has p-type with the constraint σ c † i1σ c i1σ ≤ 1.
In this paper we will present a theoretical framework describing the dynamical properties of a single exciton in a strongly correlated bilayer described by (1), following our previous shorter publication on this topic 10 . The binding of the holon and the doublon is determined by the interlayer Coulomb repulsion and we will focus on the strong coupling limit (V > t). This implies that the exciton is formed by the holon and doublon on the same rung, as is shown in Figure 1.
Understanding of the bilayer Heisenberg model will be an important step towards analysing the dynamics of a single exciton. The   methods 11,12 , dimer expansions [13][14][15] and the closely related bond operator theory 16,17 , the nonlinear sigma model 18,19 and spin wave theory [20][21][22][23] . All results indicate a O(3) universality class quantum phase transition at a critical value of J ⊥ /J from an antiferromagnetically ordered to a disordered state, see Figure 2. A naive mean field picture of the antiferromagnetic ground state is provided by the Néel state, in which each of the sublattices are occupied by either spin up or spin down electrons as shown in Figure 1. However, the exact ground state is scrambled up by spin flip interactions reducing the Néel order parameter to about 60% of its mean field value 24 . A finite interlayer coupling J ⊥ influences the antiferromagnetic order. In the limit of infinite J ⊥ , the electrons on each interlayer rung tend to form singlets destroying the antiferromagnetic order.
Standard spin wave theories however cannot account for the critical value of J ⊥ /J ∼ 2.5 found in QMC and series expansion studies. This discrepancy between numerical results and the spin wave theory has a physical Moved exciton Spin mismatch FIG. 3: Exciton motion in a naive real space picture. In a perfect Néel state, the motion of an exciton (with respect to the situation in Figure 1) causes a mismatch in the spin ordering. The kinetic energy gained by moving the exciton is proportional to the energies of the doublon te and holon t h divided by the exciton binding energy V .
origin. Chubukov and Morr 25 pointed out that standard spin wave theories do not take into account the longitudinal (that is, the interlayer) spin modes. By taking into account those longitudinal spin waves one can derive analytically the right phase diagram 26 . Another correct method is to introduce an auxiliary interaction which takes care of the hard-core constraint on the spin modes 27 .
If one wants to study the doped bilayer antiferromagnet however, one needs explicit expressions of how a moving dopant (be it a hole, electron or exciton) interacts with the spin excitations. Even though the Néel state is just an approximation to the antiferromagnetic ground state, it provides an intuitive explanation of the major role spins play in the dynamics of any dopant. As can be seen in Figure 3, a moving exciton causes a mismatch in the previously perfect Néel state. Consequently, the motion of an exciton is greatly hindered and a full understanding of possible spin wave interactions is needed to describe the exciton dynamics. This is of course similar to the motion of a single hole in a single Mott insulator layer 28,29 . Vojta and Becker 30 have computed the spectral function of a single hole in the Heisenberg bilayer. A rung linear spin wave approximation 26 is needed to obtain the expressions for the spin waves in terms of single site spin operators. Summarizing, we will formulate first an effective exciton t − J model from the bilayer t − J model in the limit of strong Coulomb attraction in section II. In order to find the interaction coefficients between excitons and spin excitations we will construct a spin wave theory of the bilayer Heisenberg model in section III. Based on these two developments, we can compute the exciton spectral function using the self-consistent Born approximation in section IV. Finally, we connect the exciton spectral function to measurable quantities in section V.

II. THE BILAYER EXCITON t − J MODEL
The bilayer t − J model (1) describes generally the p/n-doped bilayer antiferromagnet. The behavior of a bound exciton however depends on the magnitude of the FIG. 4: In weak coupling the spectrum of an exciton is obtained by the ladder diagram approximation from the spectrum of the single doped hole. The χ 0 and χ 0 are respectively the imaginary and real part of the bare exciton susceptibility. The χ and χ are the imaginary and real part of the full exciton susceptibility obtained in the ladder diagram approximation (5). Besides the continuous particle-hole spectrum above the gap, there can only be a single exciton peak determined by V χ 0 = 1 in the weak coupling limit.
Coulomb force V in H V , equation (4). If the Coulomb repulsion is relatively weak, the motion of holons and doublons will be relatively independent with each other and the H V can be treated as a perturbation on H t + H J . The full exciton-susceptibility χ(ω) can be obtained from the bare susceptibility χ 0 (ω) in the absence of the Coulomb force using the ladder diagram approximation, Since the undoped state is a Mott insulator, there is a gap in the imaginary part of the bare susceptibility χ 0 . Above this gap there is an onset of the particle-hole continuum. In the ladder diagram approximation, there can only be a single delta function peak in the full susceptibility at V χ 0 = 1 signaling the formation of an exciton. We conclude that in the weak coupling limit no special exciton features other than a single delta function peak can appear in the gap. Following our expectation that realistic materials are in fact in the strong coupling limit, as explained in section V, we will henceforth focus our attention to the strong coupling limit. In the strongly coupling limit (V t), the hopping term H t can be treated as a perturbation on the unperturbed H V using the perturbation method developed by Kato 34 . In this method, one considers first an exact solvable part of the Hamiltonian, in this case the interlayer Coulomb interaction H V . It has the eigenvalues where N is the total number of sites, N 0 is the number of dopants per layer and N is the number of double occupied sites that do not lie above a vacant site. It is clear that The essence of Kato's perturbation method is that we now forbid all states with higher H V eigenvalues. In our model, this implies that we forbid states such as the one depicted in Figure 5 where the double occupied site is not on top of the vacant site. In zeroeth order, hopping of electrons is forbidden since that would break up an exciton state. Therefore the zeroeth order Hamiltonian only contains Heisenberg terms H (0) = H J .
In second order processes are allowed that virtually break up excitons, but end up with only bound excitons. The corresponding effective Hamiltonian is given by where P e is the operator that projects out states with unbound dopants. As can be verified from Figure 5 this process allows the hopping of excitons by virtually breaking the dopants apart. If we define the exciton operator in terms of electron creation operators the exciton hopping process can be formulated as Note that in this Hamiltonian, no break-up of the exciton is required. The virtual process as described before only enabled us to relate the single layer kinetic energies to the bilayer exciton kinetic energy, Here t e is the hopping energy for a single electron, t h the hopping energy for a single hole and t is the hopping energy for a bound exciton. In addition to this hopping process there are also second order processes that equal a shift in chemical potential of the excitons. In the limit that we are interested in, that of a single exciton, we neglect chemical potential terms.
In conclusion, we formulated a model for the strong coupling limit of H V that describes the motion of bound excitons in a Mott insulator double layer. The corresponding Hamiltonian is We will refer to this model as the exciton t − J model.

A. The singlet-triplet basis
The hopping term (9) represents an exciton E i on site i swapping places with the spin background c jpσ c jnσ on site j. This Hamiltonian is in the electron Fock state representation with the background determined by the bilayer Heisenberg model (3). Historically the spin singlettriplet basis turned out to be convenient in treating the bilayer Heisenberg model, and consequently we will apply this representation also to the hopping term (9).
Unlike the fermionic holes in the single layer case, the exciton is composed of a fermionic doublon and holon in the same rung, and hence is a bosonic particle. The local Hilbert space on each interlayer rung is five dimensional with basis in terms of five hard-core bosons as one interlayer exciton state |E i and four different spin states. In the single-triplet basis, which is valid for both the doped and undoped case, we can introduce the four hard coreboson as one singlet state and three triplet states: Then the hopping term 9 can be reexpressed as: (16) We can introduce the the total spin operator and the spin difference operator Explicitly in terms of singlet and triplet rung states for In general, we see that the operator S i conserves the total onsite spin, while S always changes the total spin number s by a unit. The z-components of the spin operators do not change the magnetic number m, while the ±-components of the spin operators change the magnetic number by a unit. The bilayer Heisenberg model is now written as In conclusion, we formulated the exciton t-J model in the singlet-triplet basis which will be a starting point to solve the dynamics of the single exciton.

B. Sign problem
Notice also that the Hilbert space no longer contains fermionic degrees of freedom. The question is whether the disappearance of the fermionic structure also leads to the disappearance of the fermionic sign structure, which causes so much difficulties in the single layer t − J model 35 .
The sign structure can be investigated as follows. Remember that at half-filling the fermionic signs in the standard t − J model on a bipartite lattice can be removed by a Marshall sign transformation 36 . Upon doping, signs reappear whenever a hole is exchanged with (for example) a down spin. Which matrix elements of the Hamiltonian become positive (and thus create a minus sign in the path integral loop expansion) depends on the specific basis and on the specific Marshall sign transformation.
For the double layer exciton model, define a spin basis state with a built-in Marshall sign transformation of the form (compare to Ref. 37 ) where N ↓ An is the number of down spins on the A sublattice in the n-layer and similary we define N ↓ Bp . With these basis states the Heisenberg terms are sign-free and the only positive matrix elements come from the exchange of an exciton with a m = ±1 triplet.
We conclude that, even though the model is purely bosonic, the exciton t − J model is not sign-free and it is not possible to remove this sign structure using a Marshall or similar transformation. 52 However, as will be further elaborated upon in section IV, for both the antiferromagnetic and singlet ground states these signs do cancel out. Therefore for such ordered bilayers the problem of exciton motion turns out to be effectively bosonic.

III. UNDOPED CASE: THE BILAYER HEISENBERG MODEL
Before considering the dynamics of the exciton and expressing the interaction between the exciton and the spin background, we need to derive a spin wave theory for the bilayer Heisenberg model. Similar to the traditional Holstein-Primakoff spin-wave theory, we need a classical reference state, i.e. the mean field ground state of the bilayer Heisenberg model, then develop the linear order for the spin wave theory from the mean field ground state. The method we present here is similar to the one presented in 26 .

A. Mean field ground state
The singlet-triplet basis (23) of the bilayer Heisenberg model is convenient for mean field theory. Mean field theory tells us that for large ratio J ⊥ /J the ground state is the singlet configuration |0 0 . For small J ⊥ /J, we expect antiferromagnetic ordering, which amounts to staggered condensation of S z . By setting S z = (−1) i m we obtain a mean field Hamiltonian where S is the magnitude of spin of the spin operator on each site. A proof of this result can be found in appendix A.
The basic idea of a spin wave theory 38-40 is to start from this semiclassical (mean field) ground state and describe the local excitations with respect to this ground state. One can immediately infer why the Holstein-Primakoff or Schwinger approach to spin wave theories fails for the bilayer Heisenberg model. First, the mean field ground state is no longer a Néel state for finite α. Secondly, where Holstein-Primakoff describes one and Schwinger describes two onsite spin excitations, the bilayer Heisenberg has in fact three types of excitations. This has been pointed out by Chubukov and Morr 25 , who called the 'third' excitation the longitudinal mode.
Here we want to point out that due to the local Hilbert space and the mean field ground state as described by (25) we can 'reach' all states in the local Hilbert space with three types of excitations: a longitudinal e † which keeps the magnetic number m constant, and transversal b † ± who change the magnetic number m by either ±1.
In the limit of large S these excitations tend to become purely bosonic. We will take the mean field ground state of (25) and these three excitations as the starting point for the linear spin wave theory.
Finally, we must mention the obvious flaw in the above reasoning. Where we criticized earlier spin wave theories because they predicted the wrong critical value of J ⊥ /Jz, we now apparently adopt such a 'wrong' theory since (26) predicts α c = 1 for S = 1 2 ! Nevertheless, as we show in appendix C concerning S = 1 2 , the presence of spin waves changes the ground state energy which makes the disordered state more favorable even below the mean field critical J ⊥ Jz c calculated above. Hence, due to correctly taken the ground state energy shifts into account, one finds the accurate critical value for α consistent with numerical calculations.

B. Spin wave theory
We will now construct explicitly the spin wave theory described above for S = 1 2 . First, one needs to find the ground state following equation (25). In the S = 1 2 case, this amounts to a competition between the singlet state |s = 0, m = 0 and the triplet |s = 1, m = 0 . The mean field ground state on each rung is given by a linear superposition of those two, which interpolates between the Néel state (χ = π/4) and the singlet state (χ = 0). The onset of antiferromagnetic order can thus be viewed as the condensation of the triplet state in a singlet background. 26 With η i = (−1) i alternating we have introduced a sign change between the two sublattices A and B. The angle χ will be determined later by self-consistency conditions. The three operators that describe excitations with respect to the ground state are The e-operators will later turn out to represent the longitudinal spin waves, whereas the b-operators represent the two possible transversal spin waves. The bilayer Heisenberg model can be rewritten in terms of these operators. For completeness we include the parameter λ that enables a comparison with the Ising limit (λ = 0) with the Heisenberg limit (λ = 1), Given this, we can explicitly write down the spin operators in terms of the new e and b operators, as is done in appendix B.
From the requirement that the Hamiltonian does not contain terms linear in spin wave operators we obtain the self-consistent mean field condition for the ground state angle χ, (cos 2χ − αλ) sin 2χ = 0 (32) which has two possible solutions. Either χ = 0, which corresponds to a singlet ground state configuration, the disordered phase. If cos 2χ = αλ, there exists an antiferromagnetic ordered phase. These are indeed the two phases represented in Figure 2. Which of the two solutions ought to be chosen, depends on the ground state energy competition. In appendix C we compare the ground state energy of both phases, from which we can deduce that the critical point lies at α c ≈ 0.6, consistent with numerical literature 11,12 . The dispersion of the spin wave excitations can be found when if we consider only the quadratic terms in the Hamiltonian. This is called the 'linear' spin wave approximation, and it amounts to neglecting the cubic and quartic interaction terms. First take a Fourier transform of the spin wave operators where the sum over k runs over the 2/N momentum points in the domain [−π, π] × [−π, π] and σ = A, B represents the sublattice index. A similar definition is used for the b-operators. Upon Fourier transformation, we can decouple the spin waves from the two sublattices A and B by introducing where p = ± stand for the phase of the spin mode. Modes with p = −1 are out-of-phase and have the same dispersion as the in-phase p = 1 modes but shifted over the antiferromagnetic wavevector Q = (π, π). Again similar considerations hold for the b operators. Next we perform the Bogolyubov transformation on the magnetic excitations, The corresponding transformation angles are set by the requirement that the Hamiltonian becomes diagonal in the new operators ζ (the longitudinal spin wave) and α, β (the transversal spin wave). In doing so, we introduced the 'ideal' spin wave approximation in which we assume that the spin wave operators obey bosonic commutation relations 40 . This assumption is exact in the large S limit. For S = 1 2 this approximation turns out to work extremely well 24 , since the corrections to the bosonic commutation relations are expressed as higher order spin-wave interactions. The Bogolyubov angles are therefore given by tanh 2ϕ k,p = −p 1 2 cos 2 2χγ k sin 2 2χ + λα cos 2χ − p 1 2 cos 2 2χγ k , (38) tanh 2θ k,p = pλγ k sin 2 2χ + (1 + λ)α cos 2 χ − pλ cos 2χγ k .
The factor γ k encodes for the lattice structure, and it equals for a square lattice where the sum runs over all nearest neighbor lattice sites δ. The Bogolyuobov angles still depend on χ, which characterizes the ground state. In the antiferromagnetic phase cos 2χ = λα and for the Heisenberg limit λ = 1 these angles reduce to We can distinguish between the longitudinal and transversal spin excitations, with their dispersions given by The longitudinal spin wave is gapped and in the limit where the layers are decoupled (α = 0) completely nondispersive, while the transversal spin wave is always linear for small momentum k. This type of spectrum is similar to a phonon spectrum, which contains a linear kdependent acoustic mode and a gapped flat optical mode. This correspondence between spin waves and phonons enables us to use techniques from electron-phonon interaction studies for the exciton-spin wave interactions.
On the other hand, in the singlet phase (α > 1) one has trivially three identical triplet spin excitations. The Bogolyubov angles are given by and the dispersion of the triplet spin waves is These dispersions correspond to earlier numerical and series expansions results 13,14,25,27 . In fact, these results are exactly equal to the dispersions obtained from the nonlinear sigma model 18 .
The above derivation adds to earlier studies of the bilayer Heisenberg model in that we now found explicit expressions of how the spin waves are related to local spin In the antiferromagnetic phase (first three pictures) there is a clear distinction between the longitudinal spin waves (in green) and the transversal spin waves (in blue and red). The first is gapped, whilst the latter is zero at either k = (0, 0) or (π, π) with a linear energymomentum dependence. In the singlet phase, all spin waves are gapped triplet excitations.
flips, equations (35)- (39). This microscopic understanding of the magnetic excitations of the system enables us in the next section to derive exactly how magnetic interactions influence the dynamics of excitons.

IV. A SINGLE EXCITON IN A CORRELATED BILAYER
As was pointed out in section II, the exciton t − J model is still troubled by the sign problem even though it is purely bosonic. The sign-problem makes it difficult to say anything conclusive for systems with a finite density of excitons. Doping the single layer t − J model leads to similar loss of theoretical control, and is the consequence of the fact that the magnetic ground state changes rapidly with doping. However, we can derive the dynamics of a single exciton in the undoped bilayer. In the thermodynamic limit a single exciton will not change the ground state. Following the exciton hopping Hamiltonian (9) we can express the dynamics of the exciton upon interaction with the spin wave modes. A single exciton can be physically realized by either exciting a interlayer charge-transfer exciton in the undoped bilayer, or by infinitesimal small chemical doping of layered structures.
From a theoretical side, the spin wave we derived in the last section III can be used to constructed the effective theory of the single exciton and apply the self-consistent Born approximation. Similar to the single layer case 28 , we consider the mean field state |G as the vacuum state and it is straightforward to derive the effective theory for single exciton as: The dynamics of a single exciton are contained in the dressed Greens function, formally written as where E † k,p is the Fourier transformed exciton creation operator, and p indicates the same phase index as used for the spin waves in equation (34). The |ψ 0 denotes the ground state after the spin wave theory has been applied. The Greens function cannot be solved exactly and one needs to write out a diagrammatic expansion in the parameter t/J. For this purpose, we have derived the corresponding Feynman rules of the exciton t − J model in appendix D.
Using Dyson's equation one can rephrase the diagrammatic expansion in terms of the self-energy Σ p (k, ω) such that where p 0 (k) is the dispersion in the absence of spin excitations for the exciton with phase p. The self-energy can be computed by summing all one-particle irreducible Feynman diagrams. The degree to which exciton motion contains a free part grows with α, and indeed the free dispersion is p 0 (k) = p zt cos 2χ γ k where cos 2χ equals αλ in the antiferromagnetic phase and equals 1 in the singlet phase. As we noted before, the spin wave spectrum resembles a phonon spectrum. Hence we can compute the exciton self-energy using the Self-Consistent Born Approximation (SCBA) 28,29 , an approximation scheme developed for electron-phonon interactions but subsequently successfully applied to the single layer t − J model.
The SCBA is based on two assumptions: 1) that one can neglect vertex corrections and 2) one uses only the bare spin wave propagators. The first assumption is motivated by an extension of Migdal's theorem 53 , the second by the linear spin wave approximation. Consequently, all remaining diagrams are of the 'rainbow' type which can be summed over using a self-consistent equation. For the exciton t − J model, the SCBA amounts to computing the self-energy for the in-phase exciton, as shown diagrammatically in Figure 7. Analytically we write for the in-phase exciton self-energy, which depends on the exciton propagator and the Bogolyubov angles derived in the previous section. A similar formula to (51) applies to Σ − . However, it is easily verified that Σ − (k, ω) = Σ + (k + (π, π), ω) since γ k+(π,π) = −γ k . In general the SCBA (51) cannot be solved analytically, and hence we have obtained the exciton spectral function using an iterative procedure with Monte Carlo integration over the spin wave momenta discretized on a 32 × 32 momentum grid. We start with Σ = 0 and after approximately 20 iterations the spectral function converged. The results for typical values of α, J and t are shown in Figures 8 to 13. We start from the situation with α > 1 where the magnetic background is a disorder phase with all spin singlet configuration in the same rung. In this case, the free dispersion of the exciton with bandwidth proportional to t survived because all the magnetic triplet excitations are gapped, with an energy of Jz α(α − 1). For t < J, the exciton-magnetic interactions will barely change the free dispersion while for t > J such exciton-magnetic interactions can still occur, leading to a small 'spin polaron' effect where the exciton quasiparticle (QP) peak is diminished and spectral weight is transferred to a polaronic bump at a higher energy than the quasiparticle peak. For most values of t/J this effect is however negligible already for α just above the critical point. The exciton spectral function for t = J and α = 1.4 can be seen in Figure 8.
As α decreases towards the quantum critical point at α = 1, the gap of the triplet excitations also decreases. The effect of the exciton-magnetic interactions become more significant, which leads to an increasing transfer of spectral weight from the free coherent peak to the incoherent parts. When α hits the quantum critical point the gap of all spin excitations vanishes. There the motion of the exciton is strongly scattered by the spin excitations which completely destroy the coherent peak and leads to an incoherent critical hump in the spectrum as shown in Figure 9. For α further decreases to values α < 1, magnetic background becomes antiferromagnetically ordered with two gapless transverse modes and one gapped longitudinal mode. In this case, the motion of the exciton is still strongly scattered with the spin excitations leaving a footprints in the exciton spectrum.
Most striking thing happens then at α = 0, when the two layers are effectively decoupled and we would expect similar behavior for an interlayer exciton as for a hole or electron in a single layer. Indeed conform with the single hole in the t − J model 28,29 we find that a moving exciton causes spin frustration with an energy proportional to J. In the limit where J t the kinetic energy of the exciton is too small to be able to move through magnetic background. Therefore, we expect a localization of the exciton which is reflected in spectral data by an almost non-dispersive quasiparticle peak. This peak has a bandwidth proportional to t 2 /J and carries most of the spectral weight, 1 − O(t 2 /J 2 ). The remaining spectral weight is carried by a second peak, at an energy Jz above the main peak. Numerical results are shown in Figure 10.
More complex behavior at α = 0 arises in the antiadiabatic limit t J, where the kinetic energy of the exciton is large compared to the energy required to excite (and absorb) spin waves. Consequently, many spin waves are excited as the exciton moves and the exciton becomes 'overdressed' with multiple spin waves. At nonzero J however, a very small quasiparticle peak remains with a bandwidth of order J. Nonetheless the majority of spectral weight is carried in the incoherent many-spin wave part, as is shown in Figure 11.
Finally, realistic physical systems are expected to have a small nonzero value of α and an intermediate value  of t/J. What happens here? A simple extrapolation of the two α = 0 limits yields that the bandwidth of the quasiparticle peak will reach its maximum value at J ≈ t. Similar extrapolations suggest that about half of the spectral weight will be carried by the QP peak. However, inclusion of a finite value of α is not so trivial. Our quantitative computations in Figure 12 show a remarkable emergent Ising-like ladder spectrum, which hints at confinement of the exciton. Physically, this can be understood as follows. In the α = 0 limit, the magnetic interactions are dominated by the transverse excitations which are just single layer spin waves. The (interlayer) longitudinal spin waves become only relevant at finite α. Since these longitudinal spin waves are gapped and nondispersive, interactions with longitudinal spin waves lead to a ladder spectrum for the exciton. In first order in α, the interaction with transversal spin waves is reduced and the longitudinal interactions become dominant. The resulting Ising-like spectrum contains multiple peaks, and the energy distance between the two lowest peaks scales as t(J/t) 2/3 . The spectral weight carried by higher order peaks vanishes as t/J → 0. This is exactly the behavior seen in Figure 12. The ladder-like spectrum for small nonzero values of α is quite remarkable, and is a typical bilayer feature. To understand this notice that the interactions with the longitudinal spin wave are responsible, in first order in α, for the Ising spectrum. Take the SCBA equation (51), neglect the transversal vertex and expand the self-energy up to first order in α. Only the single spin wave diagram contributes and it equals from which we deduce, observing that Σ − = Σ + and shifting the momentum summation, that the self-energy must be momentum-independent and given by the selfconsistent equation This self-energy is exactly the same as the selfenergy of a single dopant moving through an Ising antiferromagnet 29 . In fact, in any system where a moving particle must excite a gapped and flat mode an Ising confinement spectrum will appear. Of course the exciton ladder spectrum in Figure 12 is not exactly sharp. By the above analysis, we can infer that the incoherent broadening of peaks is due to interactions with the transversal spin waves. Indeed, the transversal spin waves can be viewed as the equivalent of the single layer spin waves. Therefore for small α the effect of transversal spin waves is to reproduce the results for a single hole in the t − J model, which is quasiparticle peak broadening.

V. RELATION TO EXPERIMENT
The formation of bound exciton states can be experimentally verified in indirect measurements of the di-electric function or any other charge-excitation measurements. One particular example of the former is electron energy loss spectroscopy (EELS) which showed earlier clear signatures of the in-plane charge transfer excitons in cuprates 41,42 . The EELS cross-section is directly related to the dielectric function 43 via the dynamic structure factor S(q, ω), where the dynamic structure factor equals where the sum λ runs over all intermediate states, |ψ 0 is the . We use the dipole expansion such that where the electron position operator can be expanded in terms of the possible electron wave functions in the tight binding approximation, where |φ i are the Wannier wave functions of the electron on site i. The z component of φ i | r|φ j is proportional to the interlayer hopping energy t ⊥ , which in turn is equal to the the creation operator of an exciton, We recognize the Fourier transform of the k = 0 excitonic state, so that we find We have introduced the term e − |t| to ensure convergence of the integral so that we can integrate over t. We find that the dynamic structure factor is directly related to the exciton spectral function FIG. 13: Expected exciton spectral function for the c-axis charge-transfer exciton in YBCO bilayers. We used model parameters J = 0.125 eV, t = 0.1 eV and α = 0.04. The exciton quasiparticle peak has a dispersion with bandwidth t 2 /J, and the quasiparticle peak is the most pronounced at the line between (π, 0) and (0, π). Following at a distance of zt(J/t) 2/3 , a secondary peak develops as a sign of Ising confinement.
or in other words Consequently, one expects that the bound exciton states to show up in EELS measurements when probing the zaxis excitations. In addition to the bound exciton states, a broad electron-hole continuum will show up at high energies.
When comparing the dielectric function with the computed spectral functions in Figures 8-13, bear in mind that the latter are shifted over the energy E 0 required to excite an interlayer exciton. This energy is typically of the order of the interlayer Coulomb energy.
Another possible way to detect interlayer excitons is to use optical probes. The optical conductivity σ(q, ω) of a material is related to the dielectric function 44 by where V c (q) is the Fourier transform of the Coulomb potential 1 0 |r−r | . The real part of the c-axis optical conductivity is therefore proportional to the exciton spectral function. Similar considerations hold when one measures the Resonant Inelastic X-ray Scattering (RIXS) 45 spectrum.
How would then the exciton spectrum look like for a realistic material, such as the bilayer cuprate YBa 2 Cu 3 O 7−δ (YBCO)? Following earlier neutron scattering experiments 9,46 one can deduce that the effective J = 125 ± 5 meV and J ⊥ = 11 ± 2 meV, which corresponds to an effective value of α = 0.04α c where α c is the critical value of α. 25 . The question remains what a realistic estimate of the exciton binding energy is. The planar excitons are known to be strongly bound 42 with binding energy of the order of 1-2 eV. Since the Coulomb repulsion scales as V ∼ ( r) −1 , we can relate the binding energy of the interlayer excitons to that of the planar excitons. The distance between the layers is about twice the in-plane distance between nearest neighbor copper and oxygen atoms, but simultaneously we expect the dielectric constant c along the c-axis to be smaller than ab due to the anisotropy in the screening. Combining these two effects, we consider it a reasonable assumption that the interlayer exciton binding energy is comparable to the in-plane binding energy. The hopping energy for electrons is approximately t e = 0.4 eV which yields, together with a Coulomb repulsion estimate of V ∼ 1.5 eV, an effective exciton hopping energy of t ∼ 0.1 eV. Note that these estimates of V /t justify our use of the strong coupling limit in section II.
The spectral function corresponding to these parameters is shown in Figure 13. Since t ∼ J the ladder spectrum is strongly suppressed compared to the aforementioned anti-adiabatic limit. However, the Ising confinement still shows its signature in a small 'second ladder peak' at 0.4 eV energy above the exciton quasiparticle peak. To the best of our knowledge and to our surprise, the c-axis optical conductivity of YBCO has not been measured before. Detection of this second ladder peak in future experiments would suggest that indeed the interlayer excitons in cuprates are frustrated by the spin texture.

VI. CONCLUSION
Using a rung linear spin wave theory for the bilayer Heisenberg model we constructed a theory of strongly bound excitons in a strongly correlated bilayer system. Surprisingly, for small but finite α = J ⊥ Jz the exciton becomes confined in a fashion similar to Ising confinement. The resulting ladder spectrum should be visible in measurements of the dielectric function, such as EELS, RIXS or optical conductivity.
Possible candidate materials are for example heterostructures of n and p-type doped cuprates such as Nd 2−x Ce x CuO 4 /La 2−x Sr x CuO 4 .
In YBCO or Bi 2 Sr 2 CaCu 2 O 8+δ , the copperoxide layers come in pairs which suggests the possibility of interlayer chargetransfer excitons. A spectrum of c-axis excitons in undoped YBCO is shown in Figure 13.
Our model can be extended to different stacking structures. For example, in 214 compounds the sites in adjacent cuprate layers do not lie above each other, and we might need to include new interlayer magnetic interactions such as the Moriya-Dzyaloshinskii interaction. Different lattice structures can also be studied, of which the hexagonal lattice (as in graphene) is the most relevant.
One may wonder to what extent the used approximations are generally valid, such as the linear and ideal spin wave approximation. For the single layer Heisenberg model, it was shown that the next-to-leading order corrections where indeed significantly smaller 24 , justifying the use of both approximations in that case. Together with the fact that we were able to reproduce the known phase diagram and excitation spectrum, this suggests our approach for the bilayer Heisenberg model is justifiable. Nevertheless, an exact computation of the next-to-leading order corrections can quantify the errors of the used spin wave approximations.
Another approximation we used was the expansion in large V , the exciton coupling strength. This coupling originates in the interlayer Coulomb interaction, from which we only consider the on-site and nearest neighbor terms. Therefore our model cannot describe accurately the process of how excitons are formed out of separate doublons and holons. We think this is a very interesting open question, especially at finite temperatures. In addition, the formation process is also accompanied by an exciton annihilation process which we neglected in our current work.
Besides the interesting properties of the exciton formation process, we think that further research should be directed towards finite densities of excitons 47 . The dynamical spin-hole frustration effects that are well known in the context of doped Mott insulators occur in a strongly amplified form dealing with interlayer excitons in Mottinsulating bilayer systems. This gives further impetus to the pursuit to create such finite density correlated exciton systems in the laboratory. One can wonder whether such physics is already at work in the four-layer material Ba 2 Ca 3 Cu 4 O 8 F 2 where self-doping effects occur creating simultaneously p and n-doped layers 48 . Much effort has been devoted to create equilibrium finite exciton densities using conventional semiconductors 1 , while exciton condensation has been demonstrated in coupled semiconductor 2DEGs 2,3 . In strongly correlated heterostructures, however, formation of finite exciton densities is still far from achieved, although recent developments on oxide interfaces indicate exciting potential (see for example 49 ). Besides the closely coupled p-and n-doped conducting interface-layers in these SrTiO 3 -LaAlO 3 -SrTiO 3 heterostructures, further candidates would be closely coupled p-and n-doped cuprates, such as YBa 2 Cu 3 O 7−x or La 2−x Sr x CuO 4 with Nd 2−x Ce x CuO 4 . The feasibility of this has already been experimentally demonstrated, e.g. in 50 , but the exact interface effects need to be investigated in more detail, both experimentally as well as theoretically 47,51 . Appendix A: Large S limit bilayer Heisenberg model In this appendix we will prove equation (26). The mean field Hamiltonian (25) depends on the antiferromagnetic (AF) order parameter m. We must find the ground state energy of (25) as a function of m and then minimize with respect to m, thus yielding the mean field value of the AF order parameter.
However, since we are only interested in the critical value α c where m changes from nonzero to zero, we can proceed as follows. In the singlet phase ( m = 0) the mean field Hamiltonian is reduced to which has as ground state the singlet |0 0 and as first excited state the triplet |1 0 with energy difference E 1 − E 0 = J ⊥ . We will treat the Hamiltonian terms that depend on m as a perturbation, and compute the ground state energy in second order perturbation theory for small m. If the ground state energy decreases with nonzero m, then there is an instability towards antiferromagnetism. The perturbation Hamiltonian is and the first and second order corrections to the ground state energy are Now H (1) contains one term that is just an identity operator, and the S z operator can only change the total spin number s by one single unit. This means that the former expression yields We see that whenever α > | 1 0| S z |0 0 | 2 , the ground state energy always increases when m is nonzero. Hence the critical value of α is given by The right hand side can be evaluated explicitly using Clebsch-Gordan coefficients, since 1 0| S z |0 0 = the spin waves drive the system earlier into the singlet phase, namely at α c ≈ 0.605. For smaller values of λ this critical value increases, proportional to λ −1 .
The critical value α c = 0.605 for our spin wave theory closely resembles the numerical results of α c = 0.63. Since this ground state energy competition the system is driven into the disordered state for a different α than mean field theory suggests, we should replace bare values of α = J ⊥ Jz by the renormalized α * = α/α c when computing the exciton spectral function.
The remaining step is to write out the interaction vertices explicitly in terms of the Bogolyubov transformed spin waves. The single-magnon process equals ×p 1 (p 3 γ k2 cosh ϕ k3p3 + γ k1 sinh ϕ k3p3 ) + h.c.
These expressions are used to transform the Feynman diagrammatic representation of the SCBA of Figure 7 into the explicit formula (51).