Dynamics of a single exciton in strongly correlated bilayers

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 last decade, increasingly focused its attention on 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 two-dimensional electron gases (2DEGs) under the application of a perpendicular magnetic field [3]. Since then, many other candidate materials have been suggested that should support interlayer exciton condensation in the absence of magnetic fields, such as graphene [4–7] or topological insulators [8]. One class of candidate materials has not been considered so far, namely Mott insulators [9]. The strong interactions between electrons render 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₂ 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 a strongly correlated system such as exciton condensation, understanding the dynamics of a single exciton would be the first step.

Heterostructures of p- and n-doped cuprates can be typically described by a strongly correlated model, namely the bilayer t–J model, which is extended from two single-band t–J models for each layer with coupling terms between the layers as follows:

$$\begin{equation} H_{bt-J}=H_t+H_J+H_V, \end{equation} \tag{ 1 }$$

where H_t is the hopping of electrons in each layer,

$$\begin{equation} H_{t}=-t_{\mathrm{e}}\sum_{\langle ij \rangle\sigma,l}c_{il\sigma}^\dagger c_{jl\sigma}+\mathrm{h.c.} \end{equation} \tag{ 2 }$$

and H_J is the bilayer Heisenberg model describing the undoped Mott insulating state

$$\begin{equation} H_J=J\sum_{\langle ij \rangle,l}\mathbf{s}_{il}\cdot\mathbf{s}_{jl}+J_{\perp}\sum_{i}\mathbf{s}_{i1}\cdot\mathbf{s}_{i2}. \end{equation} \tag{ 3 }$$

Here c_ilσ and s_il denote 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 a double-occupied site (doublon) in the same rung, described by

$$\begin{equation} H_V= V\sum_{i}n_{i1}n_{i2}, \end{equation} \tag{ 4 }$$

which is the force required to form an exciton in the same rung. Without loss of generality, we assume that layer '1' contains the excess electrons with the constraint $\sum _{\sigma }c_{i1\sigma }^\dagger c_{i1\sigma }\geqslant 1$ and layer '2' has the constraint $\sum _{\sigma }c_{i2\sigma }^\dagger c_{i2\sigma }\leqslant 1$. If one considers doped systems, this amounts to n-type doping in layer '1' and p-type doping in layer '2'.

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 the doublon on the same rung, as is shown in figure 1.

Zoom In Zoom Out Reset image size

Understanding of the bilayer Heisenberg model will be an important step towards analyzing the dynamics of a single exciton. The ground state and excitations of the bilayer Heisenberg Hamiltonian have been studied quite extensively using quantum Monte Carlo (QMC) methods [11, 12], dimer expansions [13–15] and the closely related bond operator theory [16, 17], the nonlinear sigma model [18, 19] and spin-wave theory [20–23]. All results indicate an 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.

Zoom In Zoom Out Reset image size

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 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 2. 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 3. Based on these two developments, we can compute the exciton spectral function using the self-consistent Born approximation (SCBA) in section 4. Finally, we connect the exciton spectral function to measurable quantities in section 5.

Zoom In Zoom Out Reset image size

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 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 of 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 χ₀(ω) in the absence of the Coulomb force using the ladder diagram approximation (figure 4),

$$\begin{equation} \chi(\omega) = \frac{\chi_0(\omega)}{1 - V \chi_0(\omega)}. \end{equation} \noindent \tag{ 5 }$$

Since the undoped state is a Mott insulator, there is a gap in the imaginary part of the bare susceptibility χ₀''. 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 χ₀' = 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 5, we will henceforth focus our attention on the strong coupling limit.

Zoom In Zoom Out Reset image size

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 [31], in a manner similar to the derivation of the t–J model from the Hubbard model [32–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

$$\begin{equation} E_{\skew3\widetilde{N}} = V ( N - N_0 + \skew3\widetilde{N} ) = E_0 + V \skew3\widetilde{N}, \end{equation} \noindent \tag{ 6 }$$

where N is the total number of sites, N₀ is the number of dopants per layer and $\widetilde {N}$ is the number of double occupied sites that do not lie above a vacant site. It is clear that the ground state of H_V is given by the state where all double occupied and vacant sites lie above each other, as depicted in figure 1. As mentioned before, an exciton consists of a double-occupied site and a vacant site bound on top of each other. Consequently, the ground state of H_V is the state where all dopants are bound to excitons.

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 zeroth order, hopping of electrons is forbidden since that would break up an exciton state. Therefore the zeroth order Hamiltonian only contains Heisenberg terms H⁽⁰⁾ = H_J.

Zoom In Zoom Out Reset image size

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

$$\begin{equation} -\frac{1}{2V} P_{\mathrm{e}} \left( H_t \right) (1 - P_{\mathrm{e}}) \left( H_t\right) P_{\mathrm{e}}, \end{equation} \noindent \tag{ 7 }$$

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,

$$\begin{equation} E^\dagger_i = c_{i1 \uparrow}^\dagger c_{i1 \downarrow}^\dagger ( 1 - \rho_{i2} ), \end{equation} \noindent \tag{ 8 }$$

where $\rho _{i2} = \sum _\sigma c^\dagger _{i2\sigma } c_{i2\sigma }$ is the density operator in the p-type layer, the exciton hopping process can be formulated as

$$\begin{equation} H_{t,\mathrm{ex}} = - \frac{t_{\mathrm{e}}t_{\mathrm{h}}}{V} \sum_{\langle ij\rangle \sigma \sigma'} E^\dagger_j \left[ c^\dagger_{i1 \sigma'} c^\dagger_{i 2 \sigma} c_{j 2 \sigma} c_{j 1 \sigma'} \right] E_i. \end{equation} \noindent \tag{ 9 }$$

Note that in this Hamiltonian, no breakup 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,

$$\begin{equation} t = \frac{t_{\mathrm{e}} t_{\mathrm{h}}}{V}. \end{equation} \noindent \tag{ 10 }$$

Here t_e is the hopping energy for a single electron, t_h is 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 are equal to 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

$$\begin{equation} H = H_{t,\mathrm{ex}} +H_J. \end{equation} \noindent \tag{ 11 }$$

We will refer to this model as the exciton t–J model.

2.1. 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 singlet–triplet 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 a basis in terms of five hard-core bosons as one interlayer exciton state |E〉_i and four different spin states. In the singlet–triplet basis, which is valid for both the doped and the undoped case, we can introduce the four hard-core bosons as one singlet state and three triplet states:

$$\begin{equation} | 0 \; 0 \rangle_i= {1\over\sqrt2}\left(c^\dagger_{i1\uparrow}c^\dagger_{i2\downarrow}-c^\dagger_{i1\downarrow}c^\dagger_{i2\uparrow}\right) |0\rangle, \end{equation} \tag{ 12a }$$

$$\begin{equation} | 1 \; 0\rangle_i= {1\over\sqrt2}\left(c^\dagger_{i1\uparrow}c^\dagger_{i2\downarrow}+c^\dagger_{i1\downarrow}c^\dagger_{i2\uparrow}\right)|0\rangle, \end{equation} \tag{ 12b }$$

$$\begin{equation} | 1 \; 1 \rangle_i=c^\dagger_{i1\uparrow}c^\dagger_{i2\uparrow} |0\rangle, \end{equation} \tag{ 12c }$$

$$\begin{equation} | 1 \; -1 \rangle_i=c^\dagger_{i1\uparrow}c^\dagger_{i2\uparrow}|0\rangle. \end{equation}\noindent \tag{ 12d }$$

Then the hopping term (9) can be re-expressed as

$$\begin{equation} H_{t,\mathrm{ex}} = - t \sum_{\langle ij\rangle } | E_j \rangle \left( | 0 \; 0 \rangle_i \langle 0 \; 0 |_j + \sum_m |1 \; m \rangle_i \langle 1 \; m |_j \right) \langle E_i |. \end{equation} \tag{ 13 }$$

We can introduce the total spin operator

$$\begin{equation} \mathbf{S}_i = \mathbf{s}_{i1} + \mathbf{s}_{i2} \end{equation} \tag{ 14 }$$

and the spin difference operator

$$\begin{equation} \widetilde{\mathbf{S}}_i = \mathbf{s}_{i1} - \mathbf{s}_{i2}. \end{equation} \tag{ 15 }$$

Explicitly, in terms of singlet and triplet rung states for $S=\frac {1}{2}$, this reads

$$\begin{equation} S^z_i = | 1 \; 1 \rangle \langle 1 \; 1 | - | 1 \; -1 \rangle \langle 1 \; -1 |, \end{equation} \tag{ 16a }$$

$$\begin{equation} S^+_i = \sqrt{2} \left( | 1 \; 1 \rangle \langle 1 \; 0 | + | 1 \; 0 \rangle \langle 1 \; -1 | \right), \end{equation} \tag{ 16b }$$

$$\begin{equation} \widetilde{S}^z_i = - | 0 \; 0 \rangle \langle 1 \; 0 | - | 1 \; 0 \rangle \langle 0 \; 0 |, \end{equation} \tag{ 16c }$$

$$\begin{equation} \widetilde{S}^+_i = \sqrt{2} \left( | 1 \; 1 \rangle \langle 0 \; 0 | - | 0 \; 0 \rangle \langle 1 \; -1 | \right). \end{equation}\noindent \tag{ 16d }$$

In general, we see that the operator S_i conserves the total on-site spin, while $\widetilde {\mathbf {S}}_i$ 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

$$\begin{equation} H_J = \frac{J}{2} \sum_{\langle ij\rangle} \left( \mathbf{S}_i \cdot \mathbf{S}_j + \widetilde{\mathbf{S}}_i \cdot \widetilde{\mathbf{S}}_j \right) + \frac{J_\perp}{4} \sum_i \left( \mathbf{S}^2_i - \widetilde{\mathbf{S}}^2_i \right). \end{equation}\noindent \tag{ 17 }$$

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.

2.2. Sign problem

Note 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 difficulty 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, we define a spin basis state with a built-in Marshall sign transformation of the form (compare with [37])

$$\begin{equation} | \phi \rangle = (-1)^{N_{A\mathrm{n}}^\downarrow + N_{B\mathrm{p}}^\downarrow} \left| \cdots \begin{array}{ccc} \downarrow & \uparrow \downarrow & \uparrow \\ \downarrow & 0 & \downarrow \end{array} \cdots \right>, \end{equation} \tag{ 18 }$$

where N^↓_An is the number of down spins on the A sublattice in the n-layer and similarly 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 an 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². However, as will be further elaborated upon in section 4, 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.

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, and then develop the linear order for the spin-wave theory from the mean field ground state. The method we present here is similar to that presented in [26].

3.1. Mean field ground state

The singlet–triplet basis (17) of the bilayer Heisenberg model is convenient for mean field theory. Mean field theory tells us that for a 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 $\widetilde{S}_i^z$. By setting $\langle \widetilde{S}_i^z \rangle = (-1)^i \widetilde{m}$ we obtain a mean field Hamiltonian

$$\begin{equation} H^{\mathrm{MF}}_J = \sum_{i} \left[ \frac{1}{4} J z \widetilde{m}^2 + \frac{J_\perp}{4} \left( S^2_i - \widetilde{S}^2_i \right) - \frac{1}{2} J z \widetilde{m} (-1)^i \widetilde{S}^z_i \right], \end{equation}\noindent \tag{ 19 }$$

which has an ordered–disordered transition point at

$$\begin{equation} \alpha_{\mathrm{c}} \equiv \left( \frac{J_\perp}{Jz} \right)_{\mathrm{c}} = \frac{4}{3} S (S+1), \end{equation} \tag{ 20 }$$

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 [39–41] 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. Firstly, the mean field ground state is no longer a Néel state for a finite α. Secondly, where Holstein–Primakoff describes one and Schwinger describes two on-site 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 (19), we can 'reach' all states in the local Hilbert space with three types of excitations: a longitudinal e^† that keeps the magnetic number m constant and a transversal b^†_± that changes 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 (19) and these three excitations as the starting point for the linear spin-wave theory.

Finally, we must mention an 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 (20) predicts α_c = 1 for $S = \frac {1}{2}$! Nevertheless, as we show in appendix C concerning $S=\frac {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 $\left ( \frac {J_\perp }{Jz} \right )_{\mathrm {c}}$ calculated above. Hence, because of correctly taking into account the ground state energy shifts, one finds that the accurate critical value for α is consistent with numerical calculations.

3.2. Spin-wave theory

We will now construct explicitly the spin-wave theory described above for $S=\frac {1}{2}$. First, one needs to find the ground state following equation (19). In the $S=\frac {1}{2}$ case, this amounts to 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,

$$\begin{equation} | G \rangle_i = \eta_i\cos \chi | 0 \; 0 \rangle_i - \sin \chi |1 \; 0 \rangle_i, \end{equation}\noindent \tag{ 21 }$$

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)ⁱ 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

$$\begin{equation} e^\dagger_{i} = \left( \eta_{i}\sin \chi | 0 \; 0 \rangle_i + \cos \chi | 1 \; 0 \rangle_i \right) \langle G |_i , \end{equation} \tag{ 22a }$$

$$\begin{equation} b^\dagger_{i+} = | 1 \; 1 \rangle_i \langle G | _i , \end{equation} \tag{ 22b }$$

$$\begin{equation} b^\dagger_{i-} = | 1 \; -1 \rangle_i \langle G | _i. \end{equation} \noindent \tag{ 22c }$$

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 of the Ising limit (λ = 0) with the Heisenberg limit (λ = 1),

$$\begin{equation} \mathbf{S}_1 \cdot \mathbf{S}_2 = S^z_1 S^z_2 + \textstyle\frac{1}{2} \lambda ( S^+_1 S^-_2 + S^-_1 S^+_2 ). \end{equation} \tag{ 23 }$$

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 χ,

$$\begin{equation} ( \cos 2 \chi - \alpha \lambda ) \sin 2 \chi = 0, \end{equation} \tag{ 24 }$$

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 the 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

$$\begin{equation} e^\dagger_{i \sigma} = \sqrt{ \frac{2}{N} } \sum_{k} e^\dagger_{k\sigma} \mathrm{e}^{\mathrm{i} k \cdot r_i}, \end{equation} \tag{ 25 }$$

where the sum over k runs over the N/2 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

$$\begin{equation} e^\dagger_{k,p} = \frac{1}{\sqrt{2}} ( e^\dagger_{kA} + p e^\dagger_{kB} ), \end{equation} \tag{ 26 }$$

where p = ± stands 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,

$$\begin{equation} e^\dagger_{k,p} = \cosh \varphi_{k,p} \zeta^\dagger_{k,p} + \sinh \varphi_{k,p} \zeta_{-k,p}, \end{equation} \tag{ 27a }$$

$$\begin{equation} b^\dagger_{k,p,+} = \cosh \theta_{k,p} \alpha^\dagger_{k,p} + \sinh \theta_{k,p} \beta_{-k,p}, \end{equation} \tag{ 27b }$$

$$\begin{equation} b^\dagger_{k,p,-} = \cosh \theta_{k,p} \beta^\dagger_{k,p} + \sinh \theta_{k,p} \alpha_{-k,p}. \end{equation}\noindent \tag{ 27c }$$

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 waves). In doing so, we introduced the 'ideal' spin-wave approximation in which we assume that the spin-wave operators obey bosonic commutation relations [41]. This assumption is exact in the large S limit. For $S=\frac {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

$$\begin{equation} \tanh 2\varphi_{k,p} = \frac{ -p \frac{1}{2} \cos^2 2 \chi \gamma_k} {\sin^2 2\chi + \lambda \alpha \cos 2\chi -p {1\over2} \cos^2 2 \chi \gamma_k} , \end{equation} \tag{ 28 }$$

$$\begin{equation} \tanh 2\theta_{k,p} = \frac{p \lambda \gamma_k} {\sin^2 2\chi+ (1+\lambda) \alpha \cos^2 \chi -p \lambda \cos 2\chi\gamma_k}. \end{equation} \tag{ 29 }$$

The factor γ_k encodes for the lattice structure, and it equals, for a square lattice,

$$\begin{equation} \gamma_k = \frac{1}{z} \sum_{\delta} \mathrm{e}^{\mathrm{i} k \cdot \delta} = \frac{1}{2} \left( \cos k_x + \cos k_y \right), \end{equation} \tag{ 30 }$$

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

$$\begin{equation} \tanh 2\varphi_{k,p} = \frac{-p \alpha^2 \gamma_k}{ 2 - p \alpha^2 \gamma_k } , \end{equation} \noindent \tag{ 31 }$$

$$\begin{equation} \tanh 2\theta_{k,p} = \frac{p \gamma_k}{ 1+\alpha - p \alpha \gamma_k }. \end{equation} \tag{ 32 }$$

We can distinguish between the longitudinal and transversal spin excitations, with their dispersions given by

$$\begin{equation} \epsilon^{\mathrm{L}}_{k,p} = Jz \sqrt{1 - p \alpha^2 \gamma_k }, \end{equation} \noindent \tag{ 33 }$$

$$\begin{equation} \epsilon^{\mathrm{T}}_{k,p} = \frac{1}{2} Jz \sqrt{ ( 1 + \alpha ( 1 - p \gamma_k) )^2 - \gamma_k^2 }. \end{equation} \tag{ 34 }$$

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 k-dependent 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

$$\begin{equation} \tanh 2\varphi_{k,p} = - \tanh 2\theta_{k,p} = \frac{- p \gamma_k}{ 2 \alpha - p \gamma_k }, \end{equation} \tag{ 35 }$$

and the dispersion of the triplet spin waves is

$$\begin{equation} \epsilon_{k,p} = Jz \sqrt{\alpha ( \alpha - p \gamma_k) }. \end{equation} \tag{ 36 }$$

These dispersions (see figure 6) 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].

Zoom In Zoom Out Reset image size

The above derivation adds to earlier studies of the bilayer Heisenberg model in that we have now found explicit expressions of how the spin waves are related to local spin flips, equations (27a)–(29). 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 (figure 6).

As was pointed out in section 2, 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 a 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 either by exciting an interlayer charge-transfer exciton in the undoped bilayer or by infinitesimal small chemical doping of layered structures.

From a theoretical viewpoint, the spin wave theory we derived in section 3 can be used to construct the effective theory of the single exciton and apply the SCBA. 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 a single exciton as

$$\begin{equation} \fl H_{t,\mathrm{ex}} = t\sum_{\langle ij \rangle}E_j^\dagger E_i \left[\cos2\chi(1-e_i^\dagger e_j) +\sin2\chi(e_i^\dagger+e_j) -\sum_\sigma b_{i\sigma}^\dagger b_{j\sigma} \right] + \mathrm{h.c.} \end{equation} \tag{ 37 }$$

The dynamics of a single exciton are contained in the dressed Greens function, formally written as

$$\begin{equation} G^p (k, \omega) = \left\langle \psi_0 | E_{k,p} \frac{1}{\omega - H + \mathrm{i} \epsilon} E_{k,p}^\dagger | \psi_0 \right\rangle, \end{equation} \tag{ 38 }$$

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 (26). The |ψ₀〉 denotes the ground state that arises from the spin-wave approximation [24]; that is, it is defined by the conditions

$$\begin{equation} \zeta_{k,p} | \psi_0 \rangle = \alpha_{k,p} | \psi_0 \rangle = \beta_{k,p} | \psi_0 \rangle = 0 \end{equation} \tag{ 39 }$$

for all k,p. Note that |ψ₀〉 is not equal to the mean field ground state |G〉 defined in equation (21).

Now the Greens function cannot be solved exactly and one needs to write out a diagrammatic expansion in the parameter t. 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

$$\begin{equation} G^p(k, \omega) = \frac{1}{\omega - \epsilon_0^p (k) - \Sigma^p (k, \omega) + \mathrm{i}\epsilon}, \end{equation} \tag{ 40 }$$

where ^p₀(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

$$\begin{equation} \epsilon_{0}^p (k) = p \; zt \cos 2 \chi \; \gamma_k, \end{equation} \tag{ 41 }$$

where cos 2χ is equal to αλ in the antiferromagnetic phase and equal to 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 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: (i) that one can neglect vertex corrections and (ii) one uses only the bare spin-wave propagators. The first assumption is motivated by an extension of Migdal's theorem³, the second by linear spin-wave approximation. Consequently, all the remaining diagrams are of the 'rainbow' type which can be summed over using a self-consistent equation. The assumption that the vertex corrections are irrelevant allows us to completely resume Feynman diagrams up to all orders in t. The SCBA is therefore not a perturbation series expansion and consequently t does not necessarily have to be a small parameter.

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. Usual Feynman rules dictate that we need to integrate over all intermediate frequencies of the virtual spin waves. However, under the linear spin-wave approximation the spin-wave propagator is i/(ω' − (k) + i) which amounts to a Dirac delta function in the frequency domain integration [28]. For example, the first diagram of figure 7 is reduced as follows:

$$\begin{eqnarray} &&\fl \frac{1}{N} \sum_{q,p} \int^\infty_{-\infty} \frac{\mathrm{d} \omega'}{\pi} M_{k,q}^2 G^p (k-q, \omega-\omega') \left[ \frac{\mathrm{i}}{\omega' - \epsilon^{\mathrm{L}}_{k,p} + \mathrm{i}\epsilon} \right] = \frac{1}{N} \sum_{q,p} M_{k,q}^2 G^p (k-q, \omega-\epsilon^{\mathrm{L}}_{q,p}),\nonumber\\ \end{eqnarray} \noindent \tag{ 42 }$$

where M_k,q is the vertex contribution and G^p(k,ω) is the exciton propagator. Emission (or absorption) of a spin wave by an exciton can thus be incorporated by changing the momentum and energy of the exciton propagator. Analytically, we write for the in-phase exciton self-energy,

$$\begin{eqnarray} &&\fl \Sigma^+ (k, \omega) = \frac{z^2 t^2}{N} \sin^2 2\chi \sum_{q,p} \left( \gamma_{k-q} \cosh \varphi_{q,p} + p \gamma_k \sinh \varphi_{q,p} \right)^2 G^p \left(k-q, \omega-\epsilon^{\mathrm{L}}_{q,p} \right) \nonumber \\ &&+ \frac{z^2 t^2}{N^2} \cos^2 2 \chi \sum_{q,q'} \sum_{\pm, p} ( \gamma_{k+q'} \cosh \varphi_{q,p} \sinh \varphi_{q',\pm p}\pm \gamma_{k+q} \cosh \varphi_{q',\pm p} \sinh \varphi_{q,p} )^2 \nonumber \\ &&\times G^\pm \left(k-q-q', \omega - \epsilon^{\mathrm{L}}_{q,p} -\epsilon^{\mathrm{L}}_{q',\pm p} \right) + \frac{z^2 t^2}{N^2} \sum_{q,q'} \sum_{\pm, p} ( \gamma_{k-q} \cosh \theta_{q,p} \sinh \theta_{q', \pm p} \nonumber \\ &&\pm \gamma_{k-q'} \cosh \theta_{q',\pm p} \sinh \theta_{q, p} )^2 G^\pm \left(k-q-q', \omega - \epsilon^{\mathrm{T}}_{q,p} -\epsilon^{\mathrm{T}}_{q',\pm p} \right), \end{eqnarray} \noindent \tag{ 43 }$$

which depends on the exciton propagator and the Bogolyubov angles derived in the previous section. A similar formula as (43) applies as Σ⁻. However, it is easily verified that

$$\begin{equation} \Sigma^- (k, \omega) = \Sigma^+ (k+(\pi,\pi), \omega) \end{equation}\noindent \tag{ 44 }$$

since γ_k+(π,π) = −γ_k. In general the SCBA (43) cannot be solved analytically, and hence we have obtained the exciton spectral function

$$\begin{equation} A (k, \omega) = - \frac{1}{\pi} \mathrm{Im} \left[ G (k, \omega) \right] \end{equation}\noindent \tag{ 45 }$$

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–11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 \sqrt {\alpha (\alpha -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 quasi-particle peak is diminished and spectral weight is transferred to a polaronic bump at a higher energy than the quasi-particle 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. As α further decreases to values α < 1, the 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 footprint in the exciton spectrum.

The most striking thing happens then at α = 0, when the two layers are effectively decoupled, and we can expect a similar behavior for an interlayer exciton as for a hole or electron in a single layer. Indeed conforming 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 the magnetic background. Therefore, we expect a localization of the exciton, which is reflected in spectral data by an almost nondispersive quasi-particle peak. This peak has a bandwidth proportional to t²/J and carries most of the spectral weight, $1 - \mathcal {O} (t^2/J^2)$. The remaining spectral weight is carried by a second peak, at an energy Jz above the main peak.

A more complex behavior at α = 0 arises in the anti-adiabatic limit t ≫ J, where the kinetic energy of the exciton is large compared with 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 quasi-particle peak remains with a bandwidth of order J. Nonetheless, the majority of spectral weight is carried by the incoherent many-spin wave part.

However, 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 aforementioned cases yields that the bandwidth of the quasi-particle peak will reach its maximum value at J ≈ t. Similar extrapolations suggest that about half of the spectral weight will be carried by the quasi-particle peak. However, inclusion of a finite value of α is not so trivial on an analytical level. Numerical results are therefore needed, and an overview of spectral functions for different ratios of t/J and small values of α is given in figure 10.

4.1. Development of Ising-like confinement

Upon the inclusion of a small nonzero interlayer coupling α, a ladder spectrum seems to appear, reminiscent of the spectrum of a single hole in an Ising antiferromagnet. Physically, this can be understood as follows. In the α = 0 limit, the magnetic interactions are dominated by transverse excitations, which are just single layer spin waves. For any finite α > 0, the (interlayer) longitudinal spin waves become increasingly relevant. To understand their effect on the exciton spectral function, consider the SCBA equation (43), neglect the diagrams involving transversal spin waves and expand the self-energy up to first order in α. Only the single spin-wave diagram contributes and it is equal to

$$\begin{equation} \Sigma^+ (k, \omega) = \frac{z^2t^2}{N} \sum_{q,\pm} \gamma_{k-q}^2 G^{\pm} (k-q, \omega-Jz) \end{equation} \tag{ 46 }$$

from which we deduce, observing that Σ⁻ = Σ⁺ and shifting the momentum summation, that the self-energy must be momentum independent and given by the self-consistent equation

$$\begin{equation} \Sigma (\omega) = \frac{\frac{1}{2} z^2 t^2}{\omega - Jz - \Sigma (\omega - Jz)}. \end{equation} \tag{ 47 }$$

This self-energy is exactly the same as the self-energy of a single dopant moving through an Ising antiferromagnet [29]. In fact, in any system where a moving particle automatically excites a gapped and flat mode, the self-consistent equation (47) applies.

As described in [29], a hole in an Ising antiferromagnet is effectively confined by the surrounding magnetic texture. Each hop away from its initial point increases the energy, thus creating a linear potential well for the hole. In such a linear confinement potential, a ladder spectrum appears where 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 [29].

The Ising-like features in the exciton spectral function are explicitly visible in the numerically computed dispersions shown in our figure 10 and in figure 2 of [10]. We indeed conclude that the visibility of the ladder spectrum is actually enhanced in the bilayer case presented here relative to the hole in the single layer due to the nondispersive interlayer spin excitations.

Of course, the exciton ladder spectrum in figure 10 is not exactly sharp. From 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 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 quasi-particle peak broadening.

The formation of bound exciton states can be experimentally verified in indirect measurements of the dielectric function or any other charge-excitation measurements. One particular example of the former is electron energy loss spectroscopy (EELS) which earlier showed clear signatures of the in-plane charge transfer excitons in cuprates [44, 45]. The EELS cross-section is directly related to the dielectric function [46] via the dynamic structure factor S(q,ω),

$$\begin{equation} \mathrm{d}\sigma \propto \frac{1}{q^4} S(q, \omega) \propto \frac{1}{q^2} \mathrm{Im} \left[ \frac{-1}{\epsilon(q, \omega)} \right], \end{equation}\noindent \tag{ 48 }$$

where the dynamic structure factor is equal to

$$\begin{eqnarray} &&\fl S(q, \omega) = \frac{1}{N} \int \frac{\mathrm{d}t}{2\pi}\,\mathrm{e}^{- \epsilon |t|} \sum_\lambda \langle \psi_0 | \sum_i \mathrm{e}^{-\mathrm{i} q \cdot r_i} \mathrm{e}^{ \mathrm{i} (\omega - H)t} | \lambda \rangle \langle \lambda | \sum_j \mathrm{e}^{\mathrm{i} q \cdot r_j} | \psi_0 \rangle, \end{eqnarray} \tag{ 49 }$$

where the sum λ runs over all intermediate states, |ψ₀〉 is the ground state wavefunction. We use the dipole expansion such that

$$\begin{equation} \mathrm{e}^{\mathrm{i} q r_i} = 1 + \mathrm{i} \vec{q} \cdot \vec{r}_i + \cdots, \end{equation} \noindent \tag{ 50 }$$

where the electron position operator can be expanded in terms of the possible electron wave functions in the tight binding approximation,

$$\begin{equation} \sum_i \vec{r}_i = \sum_{ij \sigma} c^\dagger_{i \sigma} c_{j \sigma} \langle \phi_{i} | \vec{r} | \phi_{j} \rangle, \end{equation} \noindent \tag{ 51 }$$

where |ϕ_i〉 are the Wannier wave functions of the electron on site i. The z component of $\langle \phi _{i} | \vec {r} | \phi _{j} \rangle$ is proportional to the interlayer hopping energy t_⊥, which in turn is equal to the creation operator of an exciton,

$$\begin{equation} r^z \propto t_\perp \sum_{i \sigma} c^\dagger_{i n \sigma} c_{i p \sigma} + \mathrm{h.c.} \end{equation} \tag{ 52 }$$

$$\begin{equation} \hspace{1pc}\propto t_\perp \sum_{i} \left( E^\dagger_i + E_i \right). \end{equation} \tag{ 53 }$$

We recognize the Fourier transform of the k = 0 excitonic state, so that we find that

$$\begin{eqnarray} S(q^z, \omega) &\propto& (q^z t_\perp)^2 \int \frac{\mathrm{d}t}{2\pi} \mathrm{e}^{- \epsilon |t|} \sum_\lambda \langle \psi_0 | E_{k=0} \; \mathrm{e}^{ \mathrm{i} (\omega - H)t} | \lambda \rangle \langle \lambda | \; E_{k=0}^\dagger | \psi_0 \rangle. \end{eqnarray}\noindent \tag{ 54 }$$

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

$$\begin{eqnarray} S(q^z, \omega) &\propto& (q^z t_\perp)^2 \left\langle \psi_0 | E_{k=0} \; \left( \frac{\mathrm{i}}{\omega - H + \mathrm{i}\epsilon} -\frac{\mathrm{i}}{\omega - H - \mathrm{i}\epsilon} \right) \; E_{k=0}^\dagger | \psi_0 \right\rangle\nonumber\\\ms\ms &\propto& (q^z t_\perp)^2 A(k=0, \omega)\end{eqnarray} \tag{ 55 }$$

or in other words

$$\begin{equation} \textrm{Im} \left[ \epsilon^{-1} (q^z, \omega) \right] \sim (t_\perp)^2 A(k=0, \omega). \end{equation}\noindent \tag{ 56 }$$

Consequently, one expects the bound exciton states to show up in EELS measurements when probing the z-axis excitations. In addition to the bound exciton states, a broad electron–hole continuum will show up at high energies.

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 [47] by

$$\begin{equation} \epsilon^{-1} (q, \omega) = 1 - \mathrm{i} \frac{q^2}{\omega} V_{\mathrm{c}} (q) \sigma( q, \omega), \end{equation} \tag{ 57 }$$

where V_c(q) is the Fourier transform of the Coulomb potential $\frac {1}{\epsilon _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) [48] spectrum.

When comparing the dielectric function with the computed spectral functions in figures 8–11, bear in mind that the latter are shifted over the energy E₀ required for exciting an interlayer exciton. This energy is of the order of electronvolts. For example, along the ab-plane in cuprates charge-transfer excitons are observed in the range of 1–2 eV [49]. Since the energy required for a charge-transfer excitation is largely dependent on the on-site repulsion, we expect that the c-axis exciton will be visible at comparable energy scales.

How then would the exciton spectrum look for a realistic material, such as the bilayer cuprate YBa₂Cu₃O_7−δ (YBCO)? Following earlier neutron scattering experiments [9, 50], 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 as to what a realistic estimate of the exciton binding energy is. The planar excitons are known to be strongly bound [45] with binding energy of the order of 1–2 eV. Since the Coulomb repulsion scales as V ∼(r)⁻¹, we can relate the binding energy of the interlayer excitons to that of the planar excitons. The distance between the layers is about two times 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 2.

The spectral function corresponding to these parameters is shown in figure 11. 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 quasi-particle peak. To the best of our knowledge and to our surprise, the c-axis optical conductivity of YBCO has not been measured before in the desired regime with energies above 1 eV.⁴ Detection of this second ladder peak in future experiments would suggest that indeed the interlayer excitons in cuprates are frustrated by the spin texture.

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 $\alpha = \frac {J_\perp }{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−xCe_xCuO₄/La_2−xSr_xCuO₄. In YBCO or Bi₂Sr₂CaCu₂O_8+δ, the copper oxide layers come in pairs, which suggests the possibility of interlayer charge-transfer excitons. A spectrum of c-axis excitons in undoped YBCO is shown in figure 11.

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 were 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 that 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 from the interlayer Coulomb interaction, from which we consider only 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 that 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 this work.

Besides the interesting properties of the exciton formation process, we think that further research should be directed towards finite densities of excitons [51]. 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 Mott-insulating 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₂Ca₃Cu₄O₈F₂ where self-doping effects occur creating simultaneously p- and n-doped layers [52]. 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, the formation of finite exciton densities is still far from achieved, although recent developments on oxide interfaces indicate exciting potential (see, for example, [53]). Besides the closely coupled p- and n-doped conducting interface layers in these SrTiO₃–LaAlO₃–SrTiO₃ heterostructures, further candidates would be closely coupled p- and n-doped cuprates, such as YBa₂Cu₃O_7−x or La_2−xSr_xCuO₄ with Nd_2−xCe_xCuO₄. The feasibility of this has already been experimentally demonstrated, e.g., in [54], but the exact interface effects need to be investigated in more detail, both experimentally and theoretically [51, 55].

This research was supported by the Dutch NWO foundation through a VICI grant. The authors thank Hans Hilgenkamp, Jeroen van der Brink, Sergei Mukhin, Matthias Vojta and Dirk van der Marel for helpful discussions.

In this appendix, we will prove equation (20). The mean field Hamiltonian (19) depends on the antiferromagnetic order parameter $\widetilde {m}$. We must find the ground state energy of (19) as a function of $\widetilde {m}$ and then minimize with respect to $\widetilde {m}$, thus yielding the mean field value of the antiferromagnetic order parameter.

However, since we are only interested in the critical value α_c where $\widetilde {m}$ changes from nonzero to zero, we can proceed as follows. In the singlet phase ($\widetilde {m}=0$) the mean field Hamiltonian is reduced to

$$\begin{equation} H^{(0)} = J_\perp S_1 \cdot S_2, \end{equation} \noindent \tag{ A.1 }$$

which has as ground state the singlet |0 0〉 and as the first excited state the triplet |1 0〉 with energy difference E₁ − E₀ = J_⊥. We will treat the Hamiltonian terms that depend on $\widetilde {m}$ as a perturbation, and compute the ground state energy in second order perturbation theory for small $\widetilde {m}$. If the ground state energy decreases with nonzero $\widetilde {m}$, then there is an instability towards antiferromagnetism. The perturbation Hamiltonian is

$$\begin{equation} H^{(1)} = \textstyle\frac{1}{4} Jz \widetilde{m}^2 - \frac{1}{2} Jz \widetilde{m} (-1)^i \widetilde{S}^z \end{equation} \noindent \tag{ A.2 }$$

and the first- and second-order corrections to the ground state energy are

$$\begin{equation} E^{(1)}_0 + E^{(2)}_0 = \langle 0 \; 0 | H^{(1)} | 0 \; 0 \rangle + \sum_{s=1}^{2S} \frac{|\langle s \; 0| H^{(1)} | 0 \; 0 \rangle|^2} {E_0^{(0)} - E_s^{(0)}}. \end{equation} \noindent \tag{ A.3 }$$

Now H⁽¹⁾ contains one term that is just an identity operator, and the $\widetilde {S}^z$ operator can only change the total spin number s by a single unit. This means that the former expression yields

$$\begin{eqnarray} E^{(1)}_0 + E^{(2)}_0 &=& \frac{1}{4} Jz \widetilde{m}^2 - \frac{(Jz)^2 \widetilde{m}^2} {4 J_\perp} \; |\langle 1 \; 0| \widetilde{S}^z | 0 \; 0 \rangle|^2 \nonumber \\ &=& \frac{Jz \widetilde{m}^2}{4 \alpha} \left[ \alpha - |\langle 1 \; 0| \widetilde{S}^z | 0 \; 0 \rangle|^2 \right]. \end{eqnarray} \noindent \tag{ A.4 }$$

We see that whenever $\alpha > |\langle 1 \; 0| \widetilde {S}^z | 0 \; 0 \rangle |^2$, the ground state energy always increases when $\widetilde {m}$ is nonzero. Hence the critical value of α is given by

$$\begin{equation} \alpha_{\mathrm{c}} = |\langle 1 \; 0| \widetilde{S}^z | 0 \; 0 \rangle|^2. \end{equation} \noindent \tag{ A.5 }$$

The right-hand side can be evaluated explicitly using Clebsch–Gordan coefficients, since

$$\begin{eqnarray} \langle 1 \; 0| \widetilde{S}^z | 0 \; 0 \rangle &=& \sum_{m=-S}^S 2\,m \; C^{S S 1}_{m, -m , 0} \; C^{SS 0}_{m, -m, 0} \nonumber \\ &=& \frac{2}{2S+1} \sqrt{ \frac{3}{S (S+1)}} \sum_{m=-S}^S m^2 \nonumber \\ &=& \frac{2}{\sqrt{3}} \sqrt{S(S+1)} \end{eqnarray} \noindent \tag{ A.6 }$$

from which we indeed conclude that

$$\begin{equation} \alpha_{\mathrm{c}} = \textstyle\frac{4}{3} S (S+1). \end{equation} \tag{ A.7 }$$

The bilayer Heisenberg operators (total spin and spin difference) can be expressed, in terms of the local spin excitations e^† and b^†, by

$$\begin{equation} \fl S^z_{i \sigma} = b^\dagger_{+i\sigma} b_{+i\sigma} - b^\dagger_{-i\sigma} b_{-i\sigma}, \end{equation} \noindent \tag{ B.1 }$$

$$\begin{equation} \fl S^+_{i \sigma} = \sqrt{2} \left( - \sin \chi ( b^\dagger_{+i\sigma} + b_{-i\sigma} ) + \cos \chi ( b^\dagger_{+i\sigma} e_{i\sigma} + e^\dagger_{i\sigma} b_{-i\sigma} ) \right), \end{equation} \noindent \tag{ B.2 }$$

$$\begin{equation} \fl \widetilde{S}^z_{i\sigma} = (-1)^{\sigma_i} \left( \sin 2 \chi \left( 1 - \sum_\pm b^\dagger_{\pm i\sigma} b_{\pm i\sigma} - 2 e^\dagger_{i\sigma} e_{i\sigma} \right) - \cos 2 \chi ( e_{i\sigma}^\dagger + e_{i\sigma} ) \right), \end{equation} \noindent \tag{ B.3 }$$

$$\begin{equation} \fl \widetilde{S}^+_{i\sigma} = \sqrt{2} (-1)^{\sigma_i} \left( \cos \chi ( b^\dagger_{+i\sigma} - b_{-i\sigma} ) + \sin \chi ( b^\dagger_{+i\sigma} e_{i\sigma} - e^\dagger_{i \sigma} b_{-i\sigma} ) \right), \end{equation} \noindent \tag{ B.4 }$$

where σ represents the sign of the sublattice of site i. Consequently, the bilayer Heisenberg model in terms of these new operators reads (with $\alpha \equiv \frac {J_\perp }{Jz}$ and σ = A,B denotes the sublattice index)

$$\begin{eqnarray} &&\fl H = \frac{1}{4} Jz N (- \alpha - 2 \lambda \alpha \cos 2 \chi - \sin^2 2 \chi ) + \frac{1}{2} Jz \sum_{i} ( \cos 2 \chi - \alpha \lambda ) \sin 2 \chi (e_{i \sigma}^\dagger + e_{i \sigma} ) \nonumber \\ &&+ Jz \sum_{i} (\sin^2 2 \chi + \alpha \lambda \cos 2 \chi ) e^\dagger_{i\sigma} e_{i\sigma} - \frac{1}{2} J \sum_{i \in A, \delta} \cos^2 2 \chi ( e^\dagger_{iA} + e_{iA} )( e^\dagger_{i+\delta,B} + e_{i+\delta, B}) \nonumber \\ &&+ \frac{1}{2} J z \sum_{i \pm} ( \alpha \!+\! \sin^2 2 \chi \!+ \!\lambda \alpha \cos 2 \chi) b^\dagger_{\pm i \sigma} b_{\pm i \sigma} \!+\! \frac{1}{2} J \lambda \sum_{i \in A, \delta} ( b_{+iA}^\dagger b_{-,i+\delta,B}^\dagger \!+\! b_{-iA} b_{+,i+\delta,B} \nonumber \\ &&- \cos 2 \chi ( b^\dagger_{+iA} b_{+,i+\delta,B} \!+\! b_{-iA} b^\dagger_{-,i+\delta,B} ) \!+\! \mathrm{h.c.}) \!+\! \mathcal{O} ( b^\dagger b^\dagger e + e^\dagger b b ) + \mathcal{O} \left( [e^\dagger e + b^\dagger b]^2 \right).\nonumber\\ \end{eqnarray} \noindent \tag{ B.5 }$$

We explicitly neglect the interaction terms, which are cubic and quartic in the spin-wave operators. The above Hamiltonian contains a constant term (depends only on α, χ and λ) that describes the ground state energy competition between the singlet and antiferromagnetic phase, see appendix C. The term linear in spin operators gives us the self-consistent condition for χ. The quadratic terms will be diagonalized using the Fourier and Bogolyubov transformation as described in the main text.

The ideal spin-wave approximation introduces a shift in the ground state energy similar to that in the single layer Heisenberg model [24]. However, in the bilayer model there will be competition between the ordered phase (cos 2χ = αλ) and the disordered phase (χ = 0). Note that for αλ > 1 we automatically end up in the disordered phase.

For αλ < 1, the ground state energy of both phases is given by the expression

$$\begin{eqnarray} &&\fl E_0 = \frac{1}{4} J z N (- \alpha - 2 \alpha \lambda \cos 2 \chi - \sin^2 2 \chi) \nonumber \\ &&+ Jz \sum_k \left[ (\sin^2 2 \chi + \alpha \lambda \cos 2 \chi) \sinh^2 \varphi_k - \frac{1}{4}\cos^2 2\chi \gamma_k ( \cosh 2 \varphi_k - \sinh 2 \varphi_k ) \right] \nonumber \\ &&+ Jz \sum_k \left[\!(\alpha \!+\! \sin^2 2 \chi + \alpha \lambda \cos 2 \chi ) \sinh^2 \theta_k \!-\! \frac{1}{2} \lambda \gamma_k ( \cos 2 \chi \cosh 2 \theta_k + \sinh 2 \theta_k ) \right],\nonumber\\ \end{eqnarray} \noindent \tag{ C.1 }$$

where we have to fill in the right values of χ, θ_k and φ_k depending on the phase. As can be seen in figure C.1, the spin waves drive the system earlier into the singlet phase, namely at α_c ≈ 0.605. For smaller values of λ, this critical value increases, in proportion to λ⁻¹.

Zoom In Zoom Out Reset image size

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 $\alpha = \frac {J_\perp }{Jz}$ by the renormalized α* = α/α_c when computing the exciton spectral function.

We can rewrite the hopping Hamiltonian (13) from the singlet–triplet basis using the local spin excitation operators defined in equations (22a)–(22c),

$$\begin{equation} H_t = t \sum_{\langle ij \rangle} E_j^\dagger E_i \left( \cos 2 \chi ( 1- e^\dagger_i e_j ) + \sin 2 \chi (e_i^\dagger + e_j) - \sum_\sigma b^\dagger_{i \sigma} b_{j \sigma} \right), \end{equation} \tag{ D.1 }$$

where σ is the sum over spins ±1 and 〈ij〉 denotes nearest neighbor pairs. We need to rewrite this in terms of the longitudinal (ζ) and transversal (α and β) modes derived in the main text. Therefore we first split all operators into the ones that live on sublattice A and the ones that live on B,

$$\begin{eqnarray} &&\fl H_t = t \sum_{i \in A, \delta} E_{i+\delta,B}^\dagger E_{i,A} \left( \cos 2 \chi ( 1- e^\dagger_{i,A} e_{i+\delta,B} ) + \sin 2 \chi (e_{i,A}^\dagger + e_{i+\delta,B}) - \sum_\sigma b^\dagger_{i,A, \sigma} b_{i+\delta,B, \sigma} \right) +\mathrm{ h.c.}\nonumber\\ \end{eqnarray} \tag{ D.2 }$$

As described in the main text, the Fourier transform for the sublattice operators is $E^\dagger _{i, A} = \sqrt {\frac {2}{N}} \sum _{k} E^\dagger _{kA} \mathrm {e}^{\mathrm {i}kr_i}$ and we introduce the in-phase p = 1 and out-phase p = −1 exciton operators $E^\dagger _{k,p} = \frac {1}{\sqrt {2}} ( E^\dagger _{kA} + p E^\dagger _{kB} )$; similar expressions hold for the spinon operators. The hopping Hamiltonian can now be written as

$$\begin{equation} H_0 = zt \cos 2 \chi \sum_{k,p} p \; \gamma_k \; E^\dagger_{k,p} E_{k,p}, \end{equation} \tag{ D.3 }$$

$$\begin{equation} H_1 = \frac{zt}{\sqrt{N}} \sin 2 \chi \sum_{k,q} \sum_{p,p'} p E^\dagger_{k+q,p} E_{k,pp'} (\gamma_{k+q} e^\dagger_{-q,p'} + p' \gamma_k e_{q,p'} ), \end{equation} \tag{ D.4 }$$

$$\begin{equation} H_2^L = - \frac{zt}{N} \cos 2\chi \sum_{k,k',q} \sum_{p,p'} \sum_{\pm} p p' \gamma_{k-k'+q} E^\dagger_{k+q,p} E_{k,\pm p} e^\dagger_{k'-q,\pm p'} e_{k',p'}, \end{equation} \tag{ D.5 }$$

$$\begin{equation} H_2^{\mathrm{T}} = - \frac{zt}{N} \sum_{k,k',q} \sum_{\sigma} \sum_{p,p'} \sum_{\pm} p p' \gamma_{k-k'+q} E^\dagger_{k+q,p} E_{k,\pm p} b^\dagger_{k'-q,\pm p',\sigma} b_{k',p',\sigma}. \end{equation} \tag{ D.6 }$$

Note that this Hamiltonian contains four different types of processes. The first line H₀ contains a free part of the exciton motion. The bandwidth of the free exciton dispersion increases linearly in α in the antiferromagnetic phase until it saturates at 2zt in the disordered phase. The next term H₁ describes the creation and the annihilation of a single longitudinal mode due to exciton motion. This term is only present in the antiferromagnetic phase and is comparable to the hole–spin vertex in the single layer t − J model. Finally, there are two H₂ interactions where an exciton scatters off a transversal (H^T₂) or longitudinal (H^L₂) mode. These processes can also be changed into the creation or annihilation of a pair of spin modes. All processes can be characterized by a conservation of total phase index p and conservation of total momentum.

The remaining step is to write out the interaction vertices explicitly in terms of the Bogolyubov transformed spin waves. The single-magnon process is equal to

$$\begin{eqnarray} &&\fl H_1 = \frac{zt\sin2\chi}{\sqrt N}\sum_{k_1 \cdots k_3} \sum_{p_1 \cdots p_3} \delta^{(2)} (k_1 - k_2 - k_3) \delta\left( \prod_{i=1}^3 p_i - 1 \right) \; E_{k_1p_1}^\dagger E_{k_2p_2}\zeta_{k_3p_3} \nonumber \\ &&\times p_1 (p_3 \gamma_{k_2} \cosh \varphi_{k_3p_3} + \gamma_{k_1} \sinh \varphi_{k_3 p_3} ) +\mathrm{h.c.} \end{eqnarray} \tag{ D.7 }$$

The process that involves two longitudinal spin waves is given by

$$\begin{eqnarray} &&\fl H_2^{\mathrm{L}} = - \frac{zt \cos 2\chi}{N} \sum_{k_1 \cdots k_4} \sum_{p_1 \cdots p_4} \delta^{(2)} (k_1 - k_2 + k_3 - k_4) \delta\left( \prod_{i=1}^4 p_i - 1 \right) \; E_{k_1p_1}^\dagger E_{k_2p_2}\zeta^\dagger_{k_3p_3}\zeta_{k_4p_4} \nonumber \\&&\ms \times p_1 \left( p_4 \gamma_{k_1-k_4} \cosh \varphi_{k_3p_3} \cosh \varphi_{k_4p_4} + p_3 \gamma_{k_1-k_3} \sinh \varphi_{k_3p_3} \sinh \varphi_{k_4p_4} \right) \nonumber \\&&\ms - \frac{zt \cos 2\chi}{N} \sum_{k_1 \cdots k_4} \sum_{p_1 \cdots p_4} \delta^{(2)} (k_1 - k_2 + k_3 + k_4) \delta\left( \prod_{i=1}^4 p_i - 1 \right) \; E_{k_1p_1}^\dagger E_{k_2p_2}\zeta^\dagger_{k_3p_3}\zeta^\dagger_{k_4p_4} \nonumber \\&&\ms \times p_1 \left( p_4 \gamma_{k_1+k_4} \cosh \varphi_{k_3p_3} \sinh \varphi_{k_4p_4} + p_3 \gamma_{k_1+k_3} \sinh \varphi_{k_3p_3} \cosh \varphi_{k_4p_4} \right) + \mathrm{h.c.} \end{eqnarray} \noindent \tag{ D.8 }$$

Finally, we can also write out the Hamiltonian for the interaction vertex with the transversal spin waves. We can write this Hamiltonian term explicitly using phase and momentum conservation,

$$\begin{eqnarray} &&\fl H^{\mathrm{T}}_2 = - \frac{zt}{N} \sum_{k_1 \cdots k_4}\sum_{p_1 \cdots p_4} \delta^{(2)}( k_1 -k_2+k_3-k_4) \delta\left( \prod_{i=1}^4 p_i - 1\right ) \; \; E^\dagger_{k_1 p_1} E_{k_2 p_2} \left( \alpha^\dagger_{k_3 p_3} \alpha_{k_4 p_4} + \beta^\dagger_{k_3 p_3} \beta_{k_4 p_4} \right) \nonumber \\ &&\times p_1 \left( p_4 \gamma_{k_1 - k_4} \cosh \theta_{k_3 p_3} \cosh \theta_{k_4 p_4} + p_3 \gamma_{k_1 - k_3} \sinh \theta_{k_4 p_4} \sinh \theta_{k_3 p_3} \right) \nonumber \\ &&- \frac{zt}{N} \sum_{k_1 \ldots k_4}\sum_{p_1 \ldots p_4} \delta^{(2)}( k_1 -k_2+k_3+k_4)\delta\left( \prod_{i=1}^4 p_i - 1 \right) \; \; E^\dagger_{k_1 p_1} E_{k_2 p_2} \alpha^\dagger_{k_3 p_3} \beta^\dagger_{k_4 p_4} \nonumber \\ &&\times p_1 \left( p_4 \gamma_{k_1 + k_4} \cosh \theta_{k_3 p_3} \sinh \theta_{k_4 p_4} + p_3 \gamma_{k_1 + k_3} \cosh \theta_{k_4 p_4} \sinh \theta_{k_3 p_3} \right) + \mathrm{h.c.} \end{eqnarray} \tag{ D.9 }$$

These expressions are used to transform the Feynman diagrammatic representation of the SCBA of figure 7 into the explicit formula (43).