{'fraction_non_alphanumeric': 0.025907539559416694, 'fraction_numerical': 0.007601613403661185, 'mean_word_length': 4.73820506595452, 'pattern_counts': {'":': 1, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exciton condensation in strongly correlated electron bilayers

We studied the possibility of exciton condensation in Mott insulating bilayers. In these strongly correlated systems an exciton is the bound state of a double occupied and empty site. In the strong coupling limit the exciton acts as a hard-core boson. Its physics are captured by the exciton t-J model, containing an effective XXZ model describing the exciton dynamics only. Using numerical simulations and analytical mean field theory we constructed the ground state phase diagram. Three homogeneous phases can be distinguished: the antiferromagnet, the exciton checkerboard crystal and the exciton superfluid. For most model parameters, however, we predict macroscopic phase separation between these phases. The exciton superfluid exists only for large exciton hopping energy. Additionally we studied the collective modes and susceptibilities of the three phases. In the superfluid phase we find the striking feature that the bandwidth of the spin-triplet excitations, potentially detectable by resonant inelastic x-ray scattering (RIXS), is proportional to the superfluid density. The superfluid phase mode is visible in the charge susceptibility, measurable by RIXS or electron energy loss spectroscopy (EELS).


I. INTRODUCTION
Strongly correlated electron systems exhibit the highest attained superconducting transition temperatures currently known, and a rich variety of complex electronic phases 1,2 . Many compounds among this family of Mott insulators, such as the cuprates, are quasi-twodimensional layered materials. This renders them ideal candidates for bilayer exciton condensation, which is the topic of this publication.
The effort to achieve the condensation of excitons has a long history starting just after the discovery of BCS theory [3][4][5] . An exciton is the bound state of an electron and a hole and as such it can Bose condense. The obvious advantage of considering excitons above Cooper pairs is the strong Coulomb attraction between the electron and the hole; allowing in principle for a much higher critical temperature. To reduce the exciton lifetime problems caused by electron-hole recombination, it has been suggested to spatially separate the electrons and holes in their own subsequent layers 6,7 . This indeed has resulted in the experimental realization of exciton condensates, first in the so-called quantum Hall bilayers 8 and more recently without an externally applied magnetic field in electrically gated, optically pumped semiconductor quantum wells 9 .
The successes of exciton condensation in semiconductor 2DEG bilayer systems have led to many proposals for exciton condensation in alternative bilayer materials, such as gated topological insulators 10 or double layer graphene [11][12][13][14][15] . However, these proposals are limited to the BCS paradigm of weak coupling.
On the other hand, Mott insulators provide a com- pletely different route to exciton condensation [16][17][18] . Naively one would expect that the localization of the electrons and holes leads to a higher critical temperature, since T c is determined by the competition between the electronic kinetic energy and the electron-hole attraction. But the physics of exciton condensation in Mott insulators is in fact much richer. Instead of the picture that the electron-hole pair lives in a conduction and valence band, an exciton now consists of a double occupied and vacant site bound together on an interlayer rung, see figure 1. To estimate the binding energy, consider the in-plane charge-transfer excitons which are known to have a binding energy of the order of 1-2 eV 19 . Due to the small interlayer distances of order 1 nm we expect that a similar energy scale will set the binding of the interlayer exciton. As such, excitons in a Mott bilayer are most likely in the strongly coupled regime.
Furthermore, a single doublon-holon pair inserted into a Mott insulator leads to dynamical frustration effects 20,21 , even stronger than seen for a single hole in the t − J model 23,24 . The study of excitons in strongly correlated materials thus catches the complexity of doped Mott insulators. As we discussed elsewhere 21 the bosonic nature of the excitons actually falls short to completely eliminate all "fermion-like" signs: there are still left-over signs of the phase-string type 22 . However, it is easy to demonstrate that collinear spin order is a sufficient condition for these signs to cancel out, leaving a truly bosonic dynamics controlling the ground state and long wavelength physics. The problem thereby reduces to that of hard-core bosons (the excitons) in a sign-free spin background. This is very similar to the "spin-orbital" physics described by Kugel-Khomskii type models 25 , which can be viewed after all as describing d-d excitons interacting with spins. Also the lattice implementations 26 of the SO(5) model 27 for (cuprate) superconductivity are in this family.
Such bosonic problems can be handled with standard (semi-classical) mean field theory, and therefore the regime of finite exciton density can be addressed in an a-priori controlled manner. In most bilayer exciton setups, such as the quantum Hall bilayers or the pumped systems, there is no controllable equilibrium exciton density. In these cases one can hardly speak of the exciton density as a conserved quantity, and exciton condensation in the sense of spontaneously broken U (1) symmetry is impossible 28 . However, in Mott insulators the dopant density per layer could be fixed by, for example, chemical doping. The effective exciton chemical potential is then by definition large compared to the recombination rate. Effectively, the excitons are at finite density in equilibrium and hence spontaneous U (1) symmetry breaking is possible in the Mott insulating bilayer.
Besides the exciton superfluid phase one anticipates a plethora of competing orders, as is customary in strongly correlated materials. At zero exciton density the bilayer Heisenberg system exhibits already interesting magnetism, in the form of the antiferromagnet for small rung coupling turning via an O(3)-QNLS quantum phase transition into an "incompressible quantum spin liquid" for larger rung couplings that can be viewed as a continuation of pair singlets ("valence bonds") stacked on the rungs 29 . The natural competitor of the exciton superfluid at finite density is the exciton crystal and one anticipates that due to the strong lattice potential this will tend to lock in at commensurate densities forming exciton "Mott insulators". We will wire this in by taking also the exciton-exciton dipolar interaction into account that surely promotes such orderings. In principle there is the interesting possibility that all these orders may coexist microscopically forming an "antiferromagnetic supersolid" 30 . In this bosonic setting we can address it in a quite controlled manner, but we find that at least for the strongly coupled "small" excitons assumed here this does not happen. The reason is interesting. We already alluded to the dynamical "frustration" associated with the exciton delocalizing in the anti-ferromagnetic spin background, which is qualitatively of the same kind as for the standard "electron" t-J model. At finite densities this turns into a tendency to just phase separate on a macroscopic scale, involving antiferromagnets, exciton crystalline states and high density diamagnetic exciton superfluids, respectively.
Even though the exciton dipolar repulsion is longranged, there is no possibility of frustrated phase separation as suggested for the electronic order in cuprates [31][32][33][34][35] because the 1/r 3 interaction falls off too quickly. However, if one correctly incorporates the full exciton dipolar interaction, a variety of different exciton ordered phase may arise 36 . Here we restrict ourselves to nearest neighbor repulsion only, which allows for the formation of a checkerboard ordered exciton crystalline state.
It is disappointing that apparently in this system only conventional ground states occur. However, this is actually to a degree deceptive. The Hamiltonian describing the physics at the lattice scale describes a physics where the exciton-and spin motions are "entangled": the way in which these subsystems communicate gets beyond the notion of just being strongly coupled, since the motions of the exciton motions and the spin dynamics cannot be separated. By coarse graining this all the way to the static order parameters (the mean fields) an effective decoupling eventually results as demonstrated by the pure ground states. However, upon going "off-shell" this spin-exciton entanglement becomes directly manifest in the form of unexpected and rather counterintuitive effects on the excitation spectrum. A simple example is the zero exciton density antiferromagnet. From the rather controlled linear spin wave self-consistent Born approximation (LSW-SCBA) treatment of the one exciton problem 20 we already know that the resulting exciton spectrum can be completely different from that in a simple semiconductor. We compute here the linearized excitations around the pure antiferromagnet, recovering the LSW-SCBA result in the "adiabatic limit" where the exciton hopping is small compared to the exchange energy of the spin system, which leads to a strong enhancement of the exciton mass. In the opposite limit of fast excitons, the energy scale is recovered but the "Ising-confinement" ladder spectrum revealed by the LSW-SCBA treatment is absent. The reason is clear: in the language of this paper, the couplings between the exciton-and spin-wave modes become very big and these need to be re-summed in order to arrive at an accurate description of the exciton propagator, while our mean-field treatment corresponds with a complete neglect of these exciton-spin interactions.
The real novelty in this regard is revealed in the high density exciton superfluid phase. The spin system forms here a ground state that is a product state of pair-singlets living on the rungs. Besides the superfluid phase modes one expects in addition also the usual massive spin-triplet excitations associated with the (incompressible) singlet vacuum. The surprise is that these are characterized by a dispersion which is in part determined by the superfluid density of the exciton condensate, as we already an-nounced elsewhere 37 for which we present here the details. Counterintuitively, by measuring the spin fluctuations one can in principle determine whether the excitons are condensed in a superfluid.
Let us complete this introduction by specifying the point of departure: the Hamiltonian describing strongly bound excitons propagating through a bilayer Heisenberg spin 1/2 system. This model is derived and discussed at length in our earlier papers 20,21 and here we just summarize the outcome. Due to the strong electron-electron interactions the electronic degrees of freedom are, at electronic half-filling, reduced to spin operators s il governed by the bilayer Heisenberg model 29,38 The subscript denotes spin operators on site i in layer l = 1, 2. The Heisenberg H J is antiferromagnetic with J > 0 and J ⊥ > 0. The interlayer exciton can hop around, thereby interchanging places with the spin background. In the strong-coupling limit of exciton binding energies the exciton hopping process is described by the Hamiltonian (2) where |E is the exciton state on an interlayer rung, and |s m represent the rung spin states. Whenever an exciton hops, it effectively exchanges the spin configuration on its neighboring site. This exciton t − J model was derived earlier in Refs. 20,21 , where the optical absorption was computed in the limit of vanishing exciton density |E E| → 0. In order to study the system with a finite density of excitons, we need to enrich the current t − J model with two extra terms: a chemical potential and an exciton-exciton interaction.
The chemical potential is straightforwardly The exciton-exciton interaction requires more thought. The bare interaction between two interlayer excitons results from their electric dipole moment. Since all interlayer exciton dipole moments are pointing in the same direction the full exciton-exciton interaction is described by a repulsive 1/r 3 interaction. Hence the interaction strength decays sufficiently fast to avoid the Coulomb catastrophe responsible for frustrated phase separation 32,33 . We consider it reasonable to only include the nearest-neighbor repulsion, Here V is the energy scale associated with nearest neighbor exciton repulsion. This number can get quite high: given a typical interlayer distance 1 of 8Å and an intersite distance of 4Å the bare dipole interaction energy is 14 eV. In reality, we expect this energy to be lower due to quantum corrections and screening effects. However, the exciton-exciton interaction scale remains on the order of electronvolts and thus larger than the estimated Heisenberg J and hopping t. Let us finally consider the effects of interlayer hopping of electrons, which leads to the annihilation of excitons, This term explicitly breaks the U (1) symmetry associated with the conservation of excitons. While this term is almost certainly present in any realistic system, it is a matter of numbers whether it is relevant. In the present case of cuprates, where each layer can be doped by means of chemical substitution, we expect the chemical potential µ to be significantly larger than the interlayer tunneling t ⊥ . Consequently, the interlayer hopping is barely relevant. Throughout this publication we will discuss the effects that the inclusion of a small t ⊥ will have. The full model Hamiltonian describing a finite density of excitons in a strongly correlated bilayer is thus Let us now summarize the layout of our paper. Most of the physics of hard-core excitons on a lattice can be captured using an effective XXZ model, which is studied in section II. The ground state phase diagram of the full exciton t − J model is derived in section III, using both numerical simulations and analytical mean field theory. The excitations and the corresponding susceptibilities are discussed in section IV. We conclude this paper with a discussion on possible further lines of theoretical and experimental research in section V.

II. AN EFFECTIVE XXZ MODEL
The Hamiltonian, equation (6), has five model parameters: J, J ⊥ , t, V and µ. However, most properties of the excitons can be understood by considering the problem of hard-core bosons on a lattice. In this section we will argue that the exciton degrees of freedom can be described by an effective XXZ model. Based on some reflections on the mathematical symmetries of the full exciton t − J model, we will describe the properties of this effective XXZ model in subsection II B. We will conclude this section with an outline of the method used to obtain the excitation spectrum of the model.

A. Dynamical and symmetry algebra
Before characterizing different phases of the model we need to assess the algebraic structure of the exciton t − J model. The set of all operators that act on the local Hilbert space form the dynamical algebra, whereas the symmetries of the system are grouped together in the symmetry algebra.
To derive the dynamical algebra, it is instructive to start with the bilayer Heisenberg model which has, on each interlayer rung, a SO(4) ∼ = SU (2) × SU (2) dynamical algebra 39 . Upon inclusion of the exciton hopping term we need more operators, since now the local Hilbert space on an interlayer rung is five-dimensional (four spin states and the exciton). Consider the spin-to-exciton operator E + sm ≡ |E s m| and its conjugate E − sm = (E + sm ) † . Their commutator reads where we have introduced the operator E z sm to complete a SU (2) algebraic structure. We could set up such a construction for each of the four spin states |s m . Under these definitions the exciton hopping term, equation (2), can be rewritten in terms of an XY -model for each spin state, where the sum over sm runs over the singlet and the three triplets. Note that the exciton chemical potential, equation (3), acts as an externally applied magnetic field to this XY -model, and that the exciton-exciton repulsion, equation (4), can be rewritten as an antiferromagnetic Ising term in the E z sm operators. The dynamical algebra therefore contains four SU (2) algebras in addition to the SO(4) from the bilayer Heisenberg part. The closure of such an algebra is necessarily SU (5), which is the largest algebra possible acting on the five-dimensional Hilbert space. Hence we need a full SU (5) dynamical algebra to describe the exciton t − J model at finite density. The operators that compose this algebra are enumerated in Appendix A.
From the XY -representation of the hopping term one can already deduce that we have four distinct U (1) symmetries associated with spin-exciton exchange. The bilayer Heisenberg model contains two separate SU (2) symmetries, associated with in-phase and out-phase interlayer magnetic order. Therefore the full symmetry algebra of the model is [SU (2) Breaking of the SU (2) symmetry amounts to magnetic ordering, which is most likely antiferromagnetic (and therefore also amounts to a breaking of the lattice symmetry). Each of the U (1) algebras can be broken leading to exciton condensation. Note that next to possible broken continuous symmetries, there also might exist phases with broken translation symmetry. The checkerboard phase, already anticipated in the introduction, is an example of a phase where the lattice symmetry is broken into two sublattices.
B. What to expect: an effective XXZ model When discussing the dynamical algebra of the exciton t − J model we found that the exciton hopping terms are similar to an XY -model. The main reason is that the excitons are, in fact, hard-core bosons and thus allow for a mapping onto pseudospin degrees of freedom. Viewed as such, the exciton-exciton interaction equation (4) is similar to an antiferromagnetic Ising term and the exciton chemical potential equation (3) amounts to an external magnetic field in the z-direction. Together they form an XXZ-model in the presence of an external field, which has been investigated in quite some detail elsewhere [40][41][42][43][44][45] as well as in the context of exciton dynamics in cold atom gases 46 .
In order to understand the basic competition between the checkerboard phase and the superfluid phase of the excitons, it is worthwhile to neglect the magnetic degrees of freedom and study first this effective XXZ-model for the excitons only. The transition between the checkerboard and superfluid phases is known as the 'spin flop'transition 40 . Keeping the identification of the exciton degrees of freedom as XXZ pseudospin degrees of freedom in mind, let us review the basics of the XXZ Hamiltonian where E + = |1 0| = E x + iE y creates a hard-core bosonic particle |1 out of the vacuum |0 . This model has a built-in competition between t > 0, which favors a superfluid state, and V > 0, which favors a crystalline state where all particles are on one sublattice and the other sublattice is empty. The external field or chemical potential µ tunes the total particle density. The ground state can now be found using mean field theory. It is known that for pseudospin S = 1 2 models in (2 + 1)D the quantum fluctuations are not strong enough to defeat classical order and therefore we can rely on mean field theory, as supported by exact diagonalization studies 44 .
To find the ground state we introduce a variational wavefunction describing a condensate of excitons, The mean-field approximation amounts to choosing ψ i constant and θ i only differing between the two sublattices. We find the following mean-field energy Let's rewrite this in terms of θ = θ A + θ B and ∆θ = When |µ| ≥ 1 2 (V z + zt) the ground state is fully polarized in the z-direction. This means either zero particle density for negative µ, or a ρ = 1 for the positive µ case. Starting from the empty side, increasing µ introduces a smooth distribution of particles. This phase amounts to the superfluid phase of the excitons. The particle density on the two sublattices is equal and the total density is given by At the critical value of the chemical potential a first order transition occurs towards the checkerboard phase: the spin flop transition. In the resulting phase, which goes under various names such as the antiferromagnetic 63 , solid, checkerboard or Wigner crystalline phase, the sublattice symmetry is broken. The resulting ground state phase diagram is shown in figure  2a, where we also show the dependence of the particle density on µ. At finite temperatures in (2 + 1)d there can be algebraic long-range order. At some critical temperature a Kosterlitz-Thouless phase transition 47 will destroy this long-range order. The topology of the phase diagram however can be obtained using the finite temperature mean field theory for which we need to minimize the mean field thermodynamic potential 48 Expectation values are and the parameter m needs to be determined selfconsistently. The resulting phase diagram is shown in figure 2b, which is of the form discussed by Fisher and Nelson 41 . The first order quantum phase transition at µ c turns out to be non-trivial, a point which is usually overlooked a.

b.
FIG. 2: a. The ground state phase diagram of the XXZ model, equation (10). The graph shows the mean field particle density E z as a function of µ, with model parameters t = 1 and V = 2t. One clearly distinguishes the fully polarized phases for large µ, the superfluid phase with a linear E z vs µ dependence and the crystalline checkerboard phase with E z = 0. In between the checkerboard and the superfluid phase a non-trivial first order transition exists, with a variety of coexistence ground states with the same ground state energy. The insets show how the (E x , E z )-vectors look like in the different phases. b. Finite temperature phase diagram of the XXZ model with the same parameters. The background coloring corresponds to a semiclassical Monte Carlo computation of E z , the solid lines are analytical mean field results for the phase boundaries. We indeed see the checkerboard phase and the superfluid phase, as well as a high-temperature non-ordered 'normal' phase.
in the literature. A trivial first order transition occurs when there are two distinct phases with exactly the same energy. In the case presented here, there is a infinite set of mean field order parameters all yielding different phases yet still having the same energy. A simple analytic calculation shows that the energy of the ground state at the critical point is E c = −V z/8. Now rewrite the mean field parameters ρ A and ρ B into a sum and difference For each value of ∆ ρ with |∆ ρ | ≤ (1/2) we can find a value of ρ such that the mean field energy is exactly −V z/8. This has interesting consequences. If one can control the density instead of the chemical potential around a first order transition, in general phase separation would occur between the two competing phases. From the mean field considerations above it is unclear what would happen in a system described by the XXZ Hamiltonian, equation (10). All phases would be equally stable, at least on the mean field level, and every phase may occur in regions of any size. Such a highly degenerate state may be very sensible to small perturbations. We consider it an interesting open problem to study the dynamics of such a highly degenerate system, and whether this degeneracy may survive the inclusion of quantum corrections.
In the introduction we mentioned the existence of interlayer hopping, equation (5). Qualitatively the t ⊥ is irrelevant, which can be seen in the XXZ pseudospin language where it takes the form of a tilt of the magnetic field in the x-direction, As a result the phase diagram is shifted but not qualitatively changed. The effect of the t ⊥ on the excitation spectrum is briefly discussed in section IV B.

C. Excitations of the XXZ model
Of direct experimental relevance are the elementary excitations of a phase. The dispersion of these excitations can be computed using the 'equations of motion'-method based on the work of Zubarev 49 . We present the formalities of this method in Appendix B. In this subsection we briefly show the essence of this technique, applied to the XXZ model. Later, in section IV, we will compute the excitations for the full exciton t − J model.
The key ingredients of this Zubarev-approach are the Heisenberg equations of motion, where δ runs over all nearest neighbors. These equations cannot be solved exactly, and one relies on the approximation controlled by the mean field vacua. That is, we neglect fluctuations of the order parameters, so that products of operators on different sites are replaced by 49,50 where . . . denotes the mean field expectation value. By such a decoupling the Heisenberg equations of motion become a coupled set of linear equations which can be solved easily. In the homogeneous phase we thus obtain, after Fourier transforming, We find an analytical expression for the excitations in the superfluid phase, where γ k = 1 2 (cos k x + cos k y ). For small momenta this excitation has a linear dispersion, conform to the Goldstone theorem requiring a massless excitation as a result of the spontaneously broken U (1) symmetry. Exactly at µ = µ c the dispersion reduces to ω k = zt 1 − γ 2 k , hence the gap at k = (π, π) closes thus signaling a transition towards the checkerboard phase.
At the critical point and in the checkerboard phase, we need to take into account the fact that expectation values of operators differ on the two sublattices. The Heisenberg equations of motion now reduce to six (instead of three) linear equations, which can be straightforwardly solved. For now we postpone the discussion on the dispersion of elementary excitations to section IV, where the full exciton t−J model will be considered using the technique discussed here.

III. GROUND STATE PHASE DIAGRAM
In the previous section we have seen that the effective XXZ model predicts the existence of both an exciton superfluid phase and a checkerboard phase, separated by a first order transition. Now we derive the ground state phase diagram for the full exciton t − J model given by equation (6).
We will proceed along the same lines as in the previous section, starting with a variational wavefunction. Numerical simulation of this wavefunction creates an unbiased view on the possible inhomogeneous and homogeneous ground state phases. This serves as a basis to further analyze the phase diagram with analytical methods. The analytical mean field theory also allows us to characterize the three homogeneous phases: the antiferromagnet, the superfluid and the checkerboard crystal. Finally, combining the numerical and analytical mean field results we obtain the ground state phase diagram, see figure 7.
A. Variational wavefunction for the exciton t − J model Recall that the local Hilbert space consists of four spin states |s m and the exciton state |E . We therefore propose a variational wavefunction consisting of a product state of a superposition of all five states on each rung. For the spin states we take the SO(4) coherent state 39 which needs to be superposed with the exciton state, to obtain the total variational (product state) wavefunction This full wavefunction acts as ansatz for the numerical simulations. Note that the homogeneous phases can be described by this wavefunction with the parameters χ, θ, φ, ψ and ρ only depending on the sublattice. Given this wavefunction, the expectation value of a product of operators on different sites decouples, The only nonzero expectation values of spin operators are for S i = s i1 − s i2 and it equals wheren i is the unit vector described by the angles θ and φ. This variational wavefunction therefore assumes interlayer Néel order of magnitude sin 2χ i , which enables us to correctly interpolate between the perfect Néel order at χ = π/4 and the singlet phase χ = 0 present in the bilayer Heisenberg model. The exciton density at a rung i is trivially given by ρ i .

B. Simulated annealing
Given the variational wavefunction, we can use simulated annealing to develop an unbiased view on the possible mean field ground state phases. Therefore we start out with a lattice with on each lattice site the variables θ i , χ i , φ i , ψ i and ρ i and with periodic boundary conditions. The energy of a configuration is We performed standard Metropolis Monte Carlo updates of the lattice with fixed total exciton density. The fixed total exciton density is imposed as follows: if during an update the exciton density ρ i is changed, the exciton density on one of the neighboring sites is corrected such that the total exciton density remains constant.
The main results of the simulation are shown in figure 3, for various values of the hopping parameter t and exciton density ρ. We performed the computations on a 10 × 10 lattice. Notice that even though true long-range order does not exist in two dimensions, the correlation length of possible ordered phases is larger than the size of our simulated lattice. The other parameters are fixed at J = 125 meV, α = 0.04 and V = 2 eV. The Heisenberg couplings J = 125 meV and α = 0.04 are obtained from measurements of undoped YBCO-samples 1,51 , which we consider to be qualitatively indicative of all strongly correlated electron bilayers. The dipolar coupling is estimated at 2 eV, following our discussion in the introduction.
For each value of ρ and t we started at a high temperature T = 0.1 eV, to slowly reduce the temperature to 10 −5 eV while performing a full update of the whole lattice 10 million times. We expect that by such a slow annealing process we obtain the true ground state of equation (33), devoid of topological defects. Once we arrive at the low temperature state, we performed measurements employing 200.000 full updates of the system.
We measured six different order parameter averages: • The Néel order parameter defined by where we first sum over all spin vectors and then take the norm.
• The checkerboard order, defined as the difference in exciton density between the sublattices divided by the maximal difference possible. The maximal difference possible equals Min(ρ, 1 − ρ), so  (34), the checkerboard order parameter, equation (35), the superfluid density, equation (36), the phase coherence, equation (37), and the ratio signaling phase separation according to equation (39), 0 means complete phase separation, 1 means no phase separation. Notice that the prominent line at ρ = 0.5 signals the checkerboard phase.
• The superfluid density is given by the expectation value of the exciton operator. Here we do not make a distinction between singlet exciton condensation or triplet exciton condensation. Therefore • Now the superfluid density is not the only measure of the condensate, we can also probe the rigidity of the phase ψ. Therefore we sum up all the phase factors on all sites, If the phase is disordered, this sum tends to zero. On the other hand, complete phase coherence in the condensate phase implies that this quantity equals unity.
• Finally, we considered a measure of phase separation between the checkerboard phase and the superfluid phase. If the exciton condensate and the checkerboard phase are truly coexisting, then the maximal superfluid density attainable would be where ∆ ρ = 1 N i (−1) i ρ i . If there is phase separation however, the actual superfluid density is less than this maximal density. Therefore we also measured the ratio Ratio = Superfluid density Max SF density (39) to quantify the extent of phase separation. When this ratio is less than 1 this indicates phase separation.
The results for a full scan for the range 0 < ρ < 1 and 0 < t < 2.5 eV are shown in figure 3. In figures 4 and 5 we have displayed typical exciton density configurations for various points in the phase diagram. In combination these results suggest that there are three homogeneous phases present in the system: the antiferromagnet at low exciton densities, the exciton superfluid at high exciton hopping energies and the checkerboard crystal at halffilling of excitons. However, for most parts of the phase diagram the competition between the three phases appears to result in phase separation.
Let us investigate the phase separation in somewhat more detail. In our earlier work we found that the motion of an exciton in an antiferromagnetic background leads to dynamical frustration 20,21 . In other words: excitons do not want to coexist with antiferromagnetism. The introduction of a finite density of excitons will therefore induce b: More complicated phase separation, such as the halter form depicted here, are higher in energy. c: Starting at high temperatures with the configuration a, we slowly lowered the temperature. The resulting configuration shown here is a local minimum. d: Using the same slow annealing as for c starting from configuration b. The local energy minimum obtained this way is lower in energy than the configuration c. We conclude that even though macroscopic phase separation has the lowest energy, there are many local energy minima without macroscopic phase separation. phase separation. For large t, we find macroscopic phase separation between the antiferromagnet and the exciton superfluid, see figure 4a. At low exciton kinetic energy the excitons will crystallize in a checkerboard pattern as can be seen in figure 4b.
Close to half-filling the role of the dipole repulsion V becomes increasingly relevant. The first order 'spin flop' transition we discussed in section II B implies that there will be phase separation between the superfluid and the checkerboard order. Figures 4c and d show this phase separation. Finally there is a regime where the condensate, the checkerboard order and the Néel order are all present. However, given the dynamical frustration on the one hand and the spin-flop transition on the other hand, we again predict phase separation. A typical exciton configuration in this parameter regime is shown in figure 4e.
These simulated annealing results suggest that phase separation dominates the physics of this exciton system. To check whether the numerics are reliable we inspected directly the energies of the various homogeneous mean field solutions, using the Maxwell construction for phase separated states. The constructed phase separated configurations and their energies are shown in figure 5. The lowest energy configuration (5a) has macroscopic phase separation between the checkerboard and the antiferromagnetic phase. Intermediate states with one blob of excitons (5c) are slightly higher in energy than states with two blobs of excitons (5d). However, even though macroscopic phase separation has the lowest energy, configurations with more blobs have more entropy. Consequently for any nonzero temperatures complete macroscopic phase separation is not the most favorable solution. This is indeed seen in the numerical simulations: annealing leads to high-entropy states such as figure 5d rather than to the lowest energy configuration.
We thus conclude that the dominant phases are the antiferromagnet, the superfluid and the checkerboard. The competition between these three phases leads to phase separation in most parts of the phase diagram. The unbiased Monte Carlo simulations show the direction in which further analytical research should be directed: we will use mean field theory to characterize the three homogeneous phases.

C. Mean field theory and characterization of the phases
Given the fact that we are dealing with a hard-core boson problem, we know that mean field theory is qualitatively correct. A remaining issue is whether one can tune the exciton chemical potential rather than the exciton density in realistic experiments. Since we are prescient about the many first-order phase transitions in this system, we will perform the analysis with a fixed exciton density (the canonical ensemble). Using the Maxwell construction and the explicit µ vs. ρ relations, we can transform back to the grand-canonical ensemble.
The numerical simulations suggest that the only solutions breaking translational symmetry invoke two sublattices, and so forth for χ, θ, ψ and φ. This broken translational symmetry allows for the antiferromagnetic and exciton checkerboard order. Evaluation of the energy E = Ψ|H|Ψ of the variational wavefunction, equation (31), directly suggests that we can set θ = ψ = φ = 0 on all sites. 64 We are left with four parameters ρ A , ρ B , χ A and χ B , and as it turns out it will be more instructive to rewrite these in terms of sum and difference variables, The mean field energy per site is now given by which has to be minimized for a fixed average exciton density ρ with the constraint |∆ ρ | ≤ min(ρ, 1 − ρ). The resulting mean field phase diagram for typical values of J, J ⊥ and V , and for various t, ρ, is shown in figure 6.

Antiferromagnetic phase
As long as the exciton density is set to zero, the mean field ground state is given by the ground state of the bilayer Heisenberg model, The Néel order is given by and the energy of the antiferromagnetic state is The introduction of excitons in an antiferromagnetic background leads to dynamical frustration effects which disfavors the coexistence of excitons and antiferromagnetic order 20,21 . In fact, the numerical simulations already ruled out coexistence of superfluidity and antiferromagnetism.

Exciton condensate
For large exciton hopping energy t it becomes more favorable to mix delocalized excitons into the ground state. Due to the bosonic nature of the problem this automatically leads to exciton condensation. The delocalized excitons completely destroy the antiferromagnetic order and the exciton condensate is described by a superposition of excitons and a singlet background, Here we wish to emphasize the ubiquitous coupling to light of the superfluid. The dipole matrix element allows only spin zero transitions, and since the exciton itself is S = 0 the dipole matrix element is directly related to the superfluid density, The dipole matrix element thus acts as the order parameter associated with the superfluid phase. In most bilayer exciton condensates, such as the one in the quantum Hall regime 8 , this order parameter is also nonzero in the normal phase because of interlayer tunneling of electrons. One can therefore not speak strictly about spontaneous breaking of U (1) symmetry in such systems; there is already explicit symmetry breaking due to the interlayer tunneling. In strongly correlated electron systems the finite t ⊥ is small compared to the chemical potential µ. As discussed in the introduction, the Mott insulating bilayers now effectively allow for spontaneous U (1) symmetry breaking, and the above dipole matrix element acts as a true order parameter. Note that the irrelevance of interlayer hopping t ⊥ implies that this order parameter is, unfortunately, not reflected in photon emission or interlayer tunneling measurements. The exciton condensate is a standard two-dimensional Bose condensate. The U (1) symmetry present in the XY -type exciton hopping terms is spontaneously broken and we expect a linearly dispersing Goldstone mode in the excitation spectrum, reflecting the rigidity of the condensate. We will get back to the full excitation spectrum in section IV.
The energy of the singlet exciton condensate is and the exciton density is given by

Checkerboard phase
Whenever the exciton hopping is small, the introduction of excitons into the system leads to the 'spin flop' transition towards the checkerboard crystalline phase. As shown in the context of the XXZ model, this phase implies that one sublattice is completely filled with excitons and the other sublattice is completely empty. On the empty sublattice, any nonzero J ⊥ will guarantee that the singlet spin state has the lowest energy. Hence the average exciton density is here ρ = ∆ ρ = 1/2 and the energy of the checkerboard phase is given by It is interesting to note that the checkerboard phase is in fact similar to a Bose Mott insulator: with the new doubled unit cell we have one exciton per unit cell. The nearest neighbor dipole repulsion now acts as the 'on-site' energy preventing extra excitons per unit cell.

Coexistence of antiferromagnetism and exciton condensate
Within the analytical mean field theory set by equation (45) there exists a small region where antiferromagnetism and the exciton condensate coexist. There the energy of the homogenous coexistence phase is lower than the energy of macroscopic phase separation of the antiferromagnet and the condensate, as obtained using the Maxwell construction. However, within numerical simulations we found no evidence of coexistence. Instead, we found microscopic phase separation, which hints at a possible complex inhomogeneous phase. We therefore conclude that the homogeneous mean field theory discussed here is insufficient to find the true ground state.

Exciton Mott insulator
Finally, when the exciton density is unity we have a system composed of excitons only. In the parlance of hard-core bosons this amounts to a exciton Mott insulator. This rather featureless phase is adiabatically connected to a standard electronic band insulator: the system is now composed of two layers where each layer has an even number of electrons per unit cell. The energy of the exciton Mott insulator is, trivially

D. Phase separation
In this mean field theory most of the phase transitions are first order, with the exciton density varying discontinuously along the transition. The critical values of µ or t/J for the first order transitions are The transitions towards the coexistence region from the antiferromagnet or the condensate are second order. Additionally, the transition from the condensate to the exciton Mott insulator is second order. The critical values of t/J or µ at these second order transitions are The subscripts indicate the phases: antiferromagnetic phase (AF), coexistence phase (CO), exciton condensate (EC), exciton Mott insulator (EI), checkerboard phase (CB). For any nonzero α the first order transitions from the antiferromagnetic or coexistence phase towards the checkerboard phase are 'standard' in the sense that at the critical value of µ there are only two mean field states with equal energy. This is also true for the transitions from the antiferromagnet to the exciton condensate except at a single point. At the tricritical point separating the coexistence phase, the antiferromagnetic phase and the exciton condensate, we can set the parameters χ = 0, ∆ ρ = 0 and ∆ χ given by the value in the coexistence phase. Now the energy becomes independent of the exciton density ρ. Similarly, at the critical value of describing the transition between the checkerboard phase to the singlet exciton condensate, we can choose the mean field parameters χ = 0, ∆ χ = 0 and With these parameters, the energy becomes independent of ρ. This implies that the mean field theory predicts highly degenerate states at the critical values of µ, similar to the one we found in the XXZ model. The phase separation that thus occurs can be between an infinite set of possible ground states that have all a different exciton density. Coincidentally, the numerical simulations indicate that around the two 'degenerate' critical points indeed all the three phases are present. While the macroscopic phase separated state might have the lowest energy, figure 5 suggests that more complicated patterns of phase separation are likely to occur. The degeneracy of the critical points on the level of mean fields theory might be responsible for richer physics in these special regions of the phase diagram.

E. Conclusion
Combining the simulated annealing results of figure 3 with the analytical mean field results of figure 6 we arrive at the definitive mean field phase diagram of the exciton t − J model in figure 7. There are three main phases: the antiferromagnet at zero exciton density, the checkerboard crystal at exciton density ρ = 1/2 and the superfluid at high hopping energy t. For most parts of the phase diagram, phase separation between these three phases occurs in any possible combination. The competition between these three phases leads generally to macroscopic phase separation.
Finally, within the limitations of the semi-classical Monte Carlo approach we deduce an estimate of the transition temperature towards the superfluid state. Given a typical point in the phase diagram where the exciton condensate exists, at t = 2.5 eV and ρ = 0.18, we find a Kosterlitz-Thouless transition temperature of approximately 700 Kelvin, see figure 3c. This number should be taken not too seriously, as the exciton t − J model might not be applicable at such high temperatures given possible exciton dissociation. Additionally, at high temperatures the electron-phonon coupling becomes increasingly important, which we neglect in our exciton t − J model. Nonetheless, our estimate suggests that exciton superfluidity may extend to quite high finite temperatures.

IV. COLLECTIVE MODES AND SUSCEPTIBILITIES
Each phase of the excitons in the strongly correlated bilayer has distinct collective modes, that are in principle measurable by experiment. In order to obtain the dispersions of the collective modes we employ the technique of the Heisenberg equations of motion, introduced in the context of the XXZ model in section II C and further formalized in Appendix B. In the case of the exciton t − J model the set of equations is larger and analytical solutions can in general not be obtained. Whenever this is the case we compute the dispersions numerically.
Quantities of direct experimental relevance are the dynamical susceptibilities. We are for instance interested in the absorptive part of the dynamical magnetic susceptibility, defined by Here |ψ 0 is the ground state of the system and |n are the excited states with energy E n . It appears unlikely that bilayer exciton systems can be manufactured in bulk form which is required for neutron scattering, while there is a real potential to grow these using thin layer techniques. Therefore the detection of the dynamical spin susceptibility forms a realistic challenge for resonant inelastic Xray scattering (RIXS) 52 measurements with its claimed sensitivity for interface physics 53 .
Furthermore we are interested in the charge dynamical susceptibility (68) which is directly related to the polarization propagator. We use the operator E 00 (q) because this amounts to the interlayer dipole matrix element. Therefore, this charge dynamical susceptibility expresses the excitonic excitations. It can be observed by optical absorption experiments 54 at q = 0. Finite wavelength measurements may be obtained using the aforementioned RIXS 52 technique, or using electron energy loss spectroscopy (EELS) 55,56 . The method we use to compute the susceptibilities, based on the Heisenberg equations of motion method, is also described in Appendix B.
The three dominant phases we encountered in our mean field analysis will have distinct magnetic and optical responses. Let us briefly summarize our main findings with respect to the collective excitations. The results for the antiferromagnetic phase are shown in figures 9 to 11. This limit of vanishing exciton density has been studied with in far greater rigor than our current Zubarev method is capable of 20,21 . We can therefore compare the results of the Zubarev method with a full resummation of spin-exciton interactions using the self-consistent Born approximation (LSW-SCBA). It turns out that for small exciton kinetic hopping t the non-interacting equationsof-motion method yields reliable results. For large t one needs the full SCBA code to correctly reproduce the dynamical frustration effects of excitons in the antiferromagnetic background.
The collective modes of the exciton condensate are shown in figures 12 and 13. Due to the absence of dynamical frustration and the presence of a spin-gap we expect that these results survive in a fully interacting computation. In fact, here the modes of the simple hard-core boson system discussed in section II can be used as a template. Just as for the phase diagram, the qualitative features of XXZ model are still of relevance for the more complicated t − J model. Nonetheless, in this condensate phase the interplay between excitonic and magnetic degrees of freedom gives rise to a rather counterintuitive effect. We find that the exciton superfluid density can be detected directly in a measurement of the magnetic excitations, as we already announced elsewhere 37 .
In contrast, in the checkerboard crystalline phase the spin and exciton degrees of freedom are once again decoupled. In the remainder of this section we will elaborate further on these results for each phase separately. Throughout the following discussion, the model parameters are J = 125 meV, α = 0.04, V = 2 eV and a varying t and ρ. In order to visualize the susceptibilities we have convoluted χ with a Lorentzian of width 0.04 eV. The color scale of the susceptibility plots is in arbitrary units. A. Antiferromagnetic phase: a single exciton In the limit of zero exciton density we recover the wellknown bilayer Heisenberg physics 29 . As discussed in section III C, the spins tend to order antiferromagnetically. The excitations spectrum thus contains a Goldstone spin wave with linear dispersion around Γ and a similar mode centered around (π, π). In addition, the bilayer nature is reflected in the presence of two longitudinal spin waves with a gap of order Jz and a narrow bandwidth of order J ⊥ . The excitation spectrum and the corresponding magnetic dynamical susceptibility is shown in figure 9. Since the spin modes of the bilayer antiferromagnet are independent of any exciton degrees of freedom, we will not discuss these any further.
The dynamics of an isolated exciton in an antiferromagnetic background has been studied extensively by means of a linear spin-wave self-consistent Born approximation technique (LSW-SCBA) 20,21 . The noninteracting equations of motion method used in this paper, amounts to the complete neglect of exciton-spin interactions, while these are on the foreground of the (resummed) LSW-SCBA computation. However, the mere existence of LSW-SCBA results allows us to compare it with our current non-interacting calculations. Let us therefore first go through the LSW-SCBA results. There we need to distinguish between two limits: the adiabatic limit with t J shown in figure 10, and the antiadiabatic limit where t J shown in figure 11. Consider a single exciton in an antiferromagnetic background. Now if this exciton hops to a neighboring site, it will leave behind two spins that are ferromagnetically aligned with their neighbors. This process is called dy- . This is to be expected: in the adiabatic regime spins react much faster than the exciton motion and the exciton still moves freely dressed by a spin polaron, reducing its bandwidth to order t 2 /J. namical frustration and limits severely the motion of an exciton. In the adiabatic limit (t J) this causes the exciton bandwidth to be drastically reduced to an order t 2 /J. In addition, the magnetic background acts as a confining potential leading to small but detectable ladder states at higher energies.
At the other hand, in the anti-adiabatic regime t J exciton hopping will destroy the antiferromagnetic order as it will be surrounded by a cloud of frustrated spins. The quasiparticle picture completely breaks down and the spectral weight of the exciton is redistributed to a wide incoherent spectral bump. The ladder spectrum arising from the effective confinement will still be visible, though smeared out.
The equations-of-motion method however ignores the effects of spin-exciton interactions such as dynamical frustration. It treats the excitons as well-defined quasiparticles. As such we can already guess beforehand that the non-interacting results will be reliable in the adiabatic regime. Indeed, in the equations-of-motion method we find four exciton modes corresponding to either the singlet E + 00 or m = 0 triplet exciton E + 10 operator, just as in the LSW-SCBA. When α → 0 we can write out an  figure 10 we find four exciton bands (a.), renormalized with respect to the free hard-core boson results (b.). However, upon inclusion of the interaction the free susceptibility (c.) gets extremely renormalized (d.). The large exciton kinetic energy together with the relatively spin dynamics create an effective potential for the exciton: the exciton becomes localized and the confinement generates a ladder spectrum. Note that thus in the antiadiabatic regime the free results (a., c.) cannot be trusted.
analytical expression for the non-interacting dispersions, where each branch is twofold degenerate. This degeneracy is lifted when α = 0, leading to a splitting of order α which is largest around Γ and M . In the limit of t J the dispersions, equation (69), indeed result in an effective exciton bandwidth of order t 2 /J, conform the fully interacting theory as can be seen in figure 10. The natural question then arises: how is it possible that in the present non-interacting theory the exciton bandwidth depends on the spin parameter J? For sure, the effective exciton model introduced in section II has no such renormalization as is shown in figure 10. There the exciton bandwidth fully depends on zt.
However, it is important to realize that the exciton operators E + s0,i do not commute with the antiferromagnetic order parameter operator S z i . As a result the mean field energy of exciting an exciton is shifted either up or down (depending on the sublattice) yielding a gap between the two exciton branches of O(Jz). Now for small t, propa-gation of the exciton requires that one has to 'pay' the energy shift Jz to move through both sublattices. As a result the effective hopping is reduced by a factor t/J. Therefore the exciton bandwidth renormalization, seen in the full LSW-SCBA, is already present at the mean field level.
For large t/J however we will pay a price for the convenience of the non-interacting equations of motion method. At the mean field level one still expects the dispersions to be described by equation (69). However, upon inclusion of the interaction corrections this picture breaks down completely. The bandwidth of the non-interacting exciton is of order zt, whereas in the interacting theory an incoherent ladder spectrum of the same width arises. Thus for large t/J the non-interacting results cannot be trusted. However, this only applies to the antiferromagnetic phase due to the presence of dynamical frustration. In general, it appears that the non-interacting results are qualitatively correct in the absence of gapless modes that need to be excited in order for an exciton to move. This condition is naturally met for the other two phases. We therefore expect that exciton-spin interactions only lead to qualitative changes in the antiferromagnetic phase.
By simple selection rules one can already conclude that the singlet exciton mode couples to light. As a consequence this is the mode that is visible in the charge dynamical susceptibility, which is related to the polarization propagator. The exciton excitations are shown in figures 10d (for t < J) and 11d (for t > J).
Finally, note that at the transition from the antiferromagnetic phase to the checkerboard phase the gap in the exciton spectrum vanishes at (π, π).

B. Superfluid phase
The mode spectrum of superfluid phase, as shown in figures 12 and 13, is characterized by a linearly dispersing Goldstone mode associated with the broken U (1) symmetry. This superfluid phase mode has vanishing energy at the Γ point, where we find the inescapable linear dispersion relation The speed of the superfluid phase mode is the same as for the XXZ model in equation (28) up to a rescaling of the t and V parameters. Indeed, this speed is proportional to the superfluid density √ ρ SF = ρ(1 − ρ).
This mode can be seen in the charge susceptibility, figures 12e and f. The Goldstone mode has a gap at (π, π) which decreases monotonically with increasing exciton density. Precisely at the first order transition towards the checkerboard phase this gap closes. This mode softening at (π, π) is reminiscent of the roton in superfluid Helium: the wavelength of the roton is the same as the lattice constant of solid Helium. Next to the Goldstone mode there are two triplet excitations, shown in figure 13, each one three-fold degenerate. The degeneracy obviously arises from the standard triplet degeneracy m = −1, 0, +1. The two branches however distinguish between exciton-dominated modes and spin-dominated modes, let us discuss them separately.
The spin-dominated modes have a gap of order ∆ S = Jz α(1 + α − ρ), which is similar to the triplet gap in As the exciton condensate is spin singlet, we assume that the excitation spectrum is governed by propagating triplet modes. These modes have a gap of order J ⊥ and a bandwidth of order Jz. c,d. In contrast to the simple Heisenberg results, the actual triplet modes have enhanced kinetics 37 . The modes are split in a spin-dominated branch with small gap and large bandwidth proportional to the superfluid density (e,f.); and an excitondominated branch with a large gap and a small bandwidth (g,h.). the bilayer Heisenberg model for large α. However, the bandwidth of these excitations scales with t rather than with J, as would be customary in a system without exciton condensation (see figures 13a and b). We discussed this in great detail in recent work 37 , so let us briefly review these results. In the absence of a excitons the motion of triplets is governed by the Heisenberg superexchange yielding a bandwidth of order J. Now introduce Fock operators e † = |E 0| and t † = |1m 0|, so that the exciton-triplet exchange equation (2) This is an interaction term, thus seemingly irrelevant to the bandwidth of the triplet. However, when the exciton condensation sets in the operator e † obtains an expectation value, in fact e † = √ ρ SF where ρ SF is the condensate density. Therefore the higher order exchange term yields a quadratic triplet hopping term and the bandwidth of the triplet excitations becomes of order ztρ SF . Now remember that the exciton hopping energy t resulted, in second order perturbation theory, from the ratio t 2 e /V where t e is the electron hopping energy and V is the nearest neighbor Coulomb repulsion 20,21 . The Heisenberg coupling however was given by J = 2t 2 /U where U is the onsite Coulomb repulsion. Since for obvious reasons U > V , we find that the triplet bandwidth is enhanced whenever exciton condensation sets is. This enhancement is clearly visible in the spin susceptibility χ T , which allows for an experimental probe of the exciton superfluid density.
The other branch of triplet excitations is dominated by triplet excitons, and is therefore barely visible in the spin susceptibility and not visible in the exciton susceptibility (which only shows singlet excitons). That it is indeed dominated by triplet excitons can be inferred from computing the matrix elements of the operator E 1m , which are shown in figures 13g and h. Furthermore, the gap ∆ E = (V z+tz)ρ−µ is a function of exciton model parameters only. The bandwidth of this mode is of order O(zt), relatively independent of the exciton density. As a result, for large superfluid densities the exciton-dominated modes cross the spin-dominated triplet modes. One can directly see this in the excitation spectrum for ρ = 0.27 as shown in figure 13d.
We can compare the triplet spectrum to the mode spectrum of the singlet phase of the bilayer Heisenberg model. When J ⊥ J the ground state consists of only rung singlets. The excitation towards a triplet state, shown in figures 13a and b, has a gap Jz α(α − 1) and a bandwidth of order Jz, which is considerably smaller than the O(zt) bandwidth in the condensate. However, because the topology of the triplet mode is the same we expect that the effect of the spin-exciton interactions is the same in the bilayer Heisenberg model as for the superfluid. Since earlier LSW-SCBA showed no changes in the spectrum due to interactions, we infer that the non-interacting results for the superfluid are reliable.
To conclude our discussion of the excitations of the superfluid phase let us consider the influence of the interlayer tunneling. In the context of the XXZ model we noticed that interlayer tunneling has no qualitative influence on the phase diagram itself. However, the presence of a weak interlayer tunneling may act as potential pinningthe phase 57 opening a gap in the superfluid phase mode spectrum of order O( t ⊥ (V + t)). Persistent currents can still exist, but one needs to overcome this gap in order to get the exciton supercurrent flowing.

C. Checkerboard phase
The third homogeneous phase of the exciton t − J model is the checkerboard phase. In this phase the unit cell is effectively doubled with one exciton per unit cell. This state is analogous to a Bose Mott insulator. The trivial excitations are then the doublon and the holon: create two bosons per unit cell which costs an energy V z − µ or to remove the boson. The latter will generate a propagating exciton mode, with dispersion (73) There are two such propagating modes: one associated with the singlet exciton and one with the triplet exciton. Precisely at the transition towards the superfluid phase, one of these exciton waves becomes gapless. Note that the arguments that lead to the bandwidth renormalization in the antiferromagnetic phase also apply here, leading to an exciton bandwidth of order t 2 /V . The dispersions and the corresponding charge dynamical susceptibility can be seen in figure 14.
In the spin sector one can excite a localized spin triplet on the empty sublattice. The triplet gap is set by the interlayer energy J ⊥ , and the dispersion is flat because this triplet cannot propagate, as can be seen in figure 14b.

V. CONCLUSIONS AND DISCUSSION
We have studied the possibility of exciton condensation in strongly correlated electron bilayers. Starting from the description of the Mott state, with localized electrons, an exciton is defined as an interlayer bound state of a double occupied and vacant site. In the strong coupling limit, as of relevance to laboratory systems based on Mott insulators, the physics of such a system is described by the exciton t − J model, equation (6).
We constructed the ground state phase diagram (figure 7), based on both numerical simulations and analytical mean field theory. Three distinct phases are dominant: the antiferromagnetic phase, the checkerboard phase and the exciton condensate. For most parts of the phase diagram however, macroscopic phase separation will occur between these three phases.
Measurements of the spin and charge susceptibilities may discern in which one of the three main phases a specific system is in. The antiferromagnetic phase is characterized by a spin wave centered at (π, π), whilst in the exciton condensate the triplet bandwidth acts as a probe for the superfluid density (see figure 13 and Ref. 37 ). In the checkerboard phase the spin degrees of freedom are reflected only in a localized triplet excitation at low energy.
The charge dynamic susceptibility shows distinct qualitative behavior depending on the phase. In the antiferromagnet the spin-exciton interactions play an important role 20,21 . The superfluid phase is characterized by the visibility of the condensate Goldstone mode, whereas the checkerboard phase has propagating exciton waves with bandwidth O(t 2 /V ). Note however that since we expect phase separation to occur for most model parameters, realistic samples will likely display features from all phases in its susceptibilities.
Our theoretical work presented here is largely based on the assumption of strong coupling. In this limit, the excitons behave as local hard-core bosons. If the exciton binding energy is less dominant, the exciton will extend over more lattice sites and thus probably enable coexistence phases. On the other hand we expect that spinexciton interactions destabilize the coexistence phases, since these interactions generally lead to frustration effects. One could also wonder what happens if one includes longer-ranged interactions for the excitons, with the possibility of exciton stripes and incommensurate charge ordered phases 36 . Next, we are dealing with first order phase transitions where small changes may have severe consequences. Combining all these effects may lead to significant changes in the phase diagram, most notably in the regime where we predict phase separation.
Within the context of the strongly coupled exciton t−J model, a weaker exciton binding energy can be incorporated via interaction and hopping terms for the next nearest neighbors, next next nearest neighbors, etcetera. This might lead to complex ordered phases such as stripes 58 . Such phases are found in many strongly correlated electron systems 31,34 , and studying these in the context of the simple bosonic exciton t − J model might shed new light on the more troublesome fermionic t − J model.
In addition to stripy behavior other non-trivial exciton density profiles may occur when one considers a density imbalance between the electrons and holes. Semiconductor imbalanced systems are predicted to exhibit Fulde-Ferrell-Larkin-Ovchinnikov density modulated phases 59,60 . It is worthwhile to investigate whether such phase can exist in strongly correlated electron systems.
Next to an improvement of the phase diagram, we can also improve the susceptibilities by including the effect of exciton-spin interactions. Similar to our earlier work 20,21 on the interaction between excitons and spins in the limit of a single exciton, one could perform a diagrammatic expansion of these interactions. We expect that, apart from our earlier results in the antiferromagnetic phase, inclusion of spin-exciton interactions will not qualitatively alter the excitation spectra.
Experimentally, the close coupling of p-and n-doped Mott insulators is still relatively ill explored. However, important advances in complex oxide thin film growth, by techniques such as Molecular Beam Epitaxy (MBE) and Pulsed Laser Deposition (PLD) equipped with in situ monitoring tools such as Reflective High Energy Electron Diffraction (RHEED) are making it possible now to grow multilayers of perovskite oxides -of which many are Mott insulators -with unit cell precision. A complicating factor in fabricating multilayers of p-and n-doped perovskites, like the cuprate family from which also the high-Tc superconductors are derived, are the oftentimes conflicting (de)-oxygenation requirements. Optimized deposition and post-anneal procedures have made it possible however to make thin film contacts between n-and p-doped superconducting cuprates 61 , which is now further being explored in our labs to create and study the parallel n-p combinations resembling the theoretical model. An interesting additional system that can be included in this endeavor is the 2-dimensional electron gas that is formed at the interfaces between selected oxide bandinsulators such as SrTiO 3 and LaAlO 3 . In this respect it is noteworthy that in specific configurations, in particular a 1 unit cell SrTiO 3 layer on top of a 2 unit cell LaAlO 3 layer grown on TiO 2 -terminated SrTiO−3, a system of a closely coupled 2-dimensional electron gas and a 2-dimensional hole gas has been realized 17 .
Finally we note that some cuprate high-Tc materials appear to have an intrinsic stacking of electron-doped and hole-doped CuO 2 layers, such as Ba 2 Ca 3 Cu 4 O 8 F 2 62 , where one could look for excitonic effects. lattice sites. Such products can be decoupled within the mean field approximation as 49,50 where i and j are different lattice sites. Upon Fourier transforming lattice position into momentum and time into energy, we thus obtain a set of linear equations for the operators, ω q A (q, ω) = M (q)A (q, ω).
The spectrum of excitations is simply found by solving this eigenvalue equation for the matrix M (q). In order to find the matrix elements n|A (q)|0 that enter in susceptibilities we need will introduce the following scheme. Assume that the Hamiltonian is of the form where the sum over q runs over momenta, and n indicates the different excited states. Now α † qn is a creation operator, and irrespective of whether we are dealing with fermions or bosons we have the following equations of motion That is: every eigenvector of M (q) corresponding to a negative eigenvalue can be identified as a creation operator for one of the elementary excitations. However, the eigenvalue equation itself is not enough because it does not yield the proper normalization of α † . Since we have the eigenvector solution α † qn = U n A (q) (B6) we can write out the (anti)commutation relation for α † qn in terms of the (anti)commutation relations for the A (q). Upon requiring that on the mean field level the operators α † qn obey canonical commutation relations, that is for bosons we obtain a proper normalization for the new creation operators. We can invert the normalized matrix U n to express A (q) in terms of the creation operators α † qn . Finally, using n |α † qn |0 = δ nn we can compute the wanted matrix element for A (q).
As an example of this technique we can compute the matrix element | n|S + (q)|0 | 2 for the antiferromagnetic Heisenberg model on a square lattice. The mean field ground state is the Néel state, which leads to the following equations of motion, where the subscript A and B denote the two different sublattices, and γ q = 1 2 (cos q x + cos q y ). We quite easily infer that the eigenvalues are and thus we have one eigenvector corresponding to a creation operator, and one to an annihilation operator. If we define then the commutation relations tell us that the eigenvector matrix U must satisfy 1 = [α, α † ] = −2u 2 11 S z A − 2u 2 12 S z B = −u 2 11 + u 2 12 . (B11) The initial S + q operator, which enters in the spin susceptibility, can be expressed in terms of the eigenvector matrix as Some straightforward algebra now yields n|S + (q)|0 2 = 1 2 which is the same susceptibility one can obtain by using the Holstein-Primakoff linear spin wave approximation. The approximation scheme we introduced here can therefore be viewed as a generalization of the linear spin wave approximation.