{'fraction_non_alphanumeric': 0.040642266473838395, 'fraction_numerical': 0.01071680201728038, 'mean_word_length': 4.533032416906032, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spin-Mediated Mott Excitons

Motivated by recent experiments on Mott insulators, in both iridates and ultracold atoms, we theoretically study the effects of magnetic order on the Mott-Hubbard excitons. In particular, we focus on spin-mediated doublon-holon pairing in Hubbard materials. We use several complementary theoretical techniques: mean-field theory to describe the spin degrees of freedom, the self-consistent Born approximation to characterize individual charge excitations across the Hubbard gap, and the Bethe-Salpeter equation to identify bound states of doublons and holons. The binding energy of the Hubbard exciton is found to increase with increasing the N{\'e}el order parameter, while the exciton mass decreases. We observe that these trends rely significantly on the retardation of the effective interaction, and require consideration of multiple effects from changing the magnetic order. Our results are consistent with the key qualitative trends observed in recent experiments on iridates. Moreover, the findings could have direct implications on ultracold atom Mott insulators, where the Hubbard model is the exact description of the system and the microscopic degrees of freedom can be directly accessed.

For example, it is known that Mott insulators are often antiferromagnetic at low temperature, but very little work has been done to understand how and to what extent the presence of such order affects exciton properties. The qualitative role of magnetization is sketched in Fig. 1 -charges remain bound so as to minimize the number of spins disrupted by their motion -but a quantitative description has been lacking.
Recent experiments have begun to investigate this question. In Refs. [22,23], pump-probe experiments were performed on the Mott insulator Na 2 IrO 3 both with and without magnetic order (controlled by varying temperature or applying an intermediate pulse). The authors concluded that the binding energy and exciton mass are both enhanced by the presence of magnetization. Ref. [24] similarly observed that the binding energy increases with the spin-spin interaction strength in cuprates. See also Ref. [25], which found that the relaxation time in Mott insulators decreases with increasing The same question can apply to Mott insulators in synthetic quantum systems, such as ultracold gases. By loading fermionic atoms into an optical lattice and tuning their interactions, the Fermi-Hubbard model can be synthesized experimentally [26][27][28][29][30]. Unlike condensed matter systems, such as the iridates, neutral fermionic atoms in an optical lattice are genuinely described by the Hubbard Hamiltonian, without any additional effects arising from longer-range Coulomb interactions, phonons, etc. Researchers have quite recently begun investigating the interplay of spin and charge degrees of freedom in this setting [31][32][33][34][35][36][37] (note that here the "charge" excitations are not actually charged).
In this paper, we perform a theoretical study of the role of magnetic order in Mott excitons, the first such to our knowledge. As depicted in Fig. 1 by Coulomb interaction in two aspects. First, the spinmediated interaction is far from instantaneous, and second, the individual charges are themselves renormalized by spin fluctuations. We shall demonstrate that both effects are necessary ingredients in the trends reported here.
Given the complexity of the problem, our analysis requires multiple stages. We first use slave particles to isolate spin and charge degrees of freedom, then describe the spin dynamics by mean-field theory, calculate the dispersion of charges self-consistently, and finally characterize excitonic states via the Bethe-Salpeter equation. Many of the steps in this program are analogous to those in Ref. [38], which studied charge dynamics in the Hubbard model. The good agreement between the results of Ref. [38] and alternate numerical methods lends support to the present approach.
Our key finding is that larger magnetization leads to an increased binding energy of the Hubbard exciton but a decreased mass. This observation is in some tension with interpretations of recent experiments [23]. It also stands in contrast to conventional Coulomb-mediated excitons, where the binding energy and mass are proportional to each other.
Note that the formation of Mott excitons is closely related to the Cooper pairing of holes in high-T c superconductors. Similar treatments of hole-hole binding can be found in the corresponding literature [39][40][41][42][43]. Nonetheless, there are differences between the two problems, as we discuss below. Moreover, to the best of our knowledge, a study of how the bound state properties change as a function of magnetization has not yet been carried out.
In the following Sec. II, we describe the steps of our analysis in detail. Results are presented in Sec. III, and Sec. IV concludes.

II. FORMALISM & METHODS
Our starting point is the 2D Fermi-Hubbard model, which by now needs no introduction: where σ ∈ {↑, ↓} and ij denotes nearest-neighbor sites on a square lattice. c iσ is the usual electron annihilation operator and n iσ ≡ c † iσ c iσ . We shall consider the system at half filling in the U t limit. It is well-known that in this limit, the Hubbard model features two types of excitations, associated with the transport of charge and spin respectively [44][45][46][47]. Furthermore, the charge excitations can be either positive or negative, corresponding to sites with zero or two electrons, and their creation comes with a large energy cost of order U . By analogy with conventional semiconductors, we thus expect this system to support well-defined excitons in the dilute-charge limit. However, long-wavelength spin excitations do not come with an energy cost, and their presence plays a significant role in determining the exciton properties.
There are many formalisms with which to study the Hubbard model [48][49][50][51]. Since our focus is on the motion of only a few charges within a background of spin excitations, the slave-particle formalism is particularly well-suited [38,52,53]. The steps of our calculation are as follows: i) Express the Hamiltonian in terms of slave particles -doublons, holons, & spinons -and reduce to the t-J model following the standard procedure [54].
ii) Treat the Heisenberg interaction within semiclassical and mean-field approximations, while neglecting the back action of doublons and holons on the magnetic order.
iii) Calculate the dispersion of individual doublons and holons in the magnetic background via the selfconsistent Born approximation.
iv) Calculate exciton properties using the Bethe-Salpeter equation.
The major limitation of this program is our approximate description of the magnetic order. Thus we do not claim to have quantitatively accurate results, especially at small magnetization. That said, we do expect that the qualitative trends seen here are accurate, including near the equilibrium value of magnetization, for which mean field theory is known to work reasonably well (see Ref. [38] and references therein).

A. Slave particles
In the slave-particle formalism, we express the electron operator as (with σ = ±1) where d i and e i are fermionic operators and s iσ is bosonic. One can confirm that Eq. (2) is consistent with the commutation relations. A site with a d particle is to be interpreted as a site with two electrons (a "doublon"), a site with an e particle is to be interpreted as an empty site (a "holon"), and a site with an s σ particle is one with a single electron having spin σ (a "spinon"). See Fig. 2. The physical content of Eq. (2) is then clear: removing an electron of given spin is equivalent to replacing the doublon with the opposite spinon if the site is doubly-occupied and replacing the spinon with a holon if the site is singly-occupied (otherwise the state is annihilated). Note that since every site is in one of the four states -empty, spin-up, spin-down, doubly-occupiedthere must be exactly one of the fictitious particles on each site: Sketch of the slave-particle formalism. Each of the four possible occupancies on a site corresponds to a different type of fictitious particle.
The original Hamiltonian clearly preserves this relationship. Substituting Eq. (2) into Eq. (1), we have that Note that the first line preserves the number of doublons and holons, whereas the second line does not. At large U , the second line of Eq. (4) can be treated by perturbation theory in t/U . The method as applied here is standard, and can be found in, e.g., Ref. [54]. We obtain the t-J model: where J ≡ 4t 2 /U . Strictly speaking, Eq. (5) should include additional next-nearest-neighbor terms, as well as a direct interaction between nearest-neighbor doublons and holons, but these are commonly neglected.

B. Magnetic ordering
The second line of Eq. (5) is precisely the antiferromagnetic Heisenberg Hamiltonian, expressed in terms of Schwinger bosons (here the spinons s iσ ) [54]. We are specifically considering Néel-ordered ground states, with positive magnetization on the A sublattice and negative on B. It is thus convenient to express s i↑ (s i↓ ) in terms of s i↓ (s i↑ ) for i ∈ A (B), using Eq. (3) -this is equivalent to using Holstein-Primakoff rather than Schwinger bosons [54]: We neglect d † i d i and e † i e i in using Eq. (3) because we are interested in the dilute-charge limit. From here on, we shall simply write s i in place of s i↓ (s i↑ ) for i ∈ A (B). On both sublattices, the s i boson represents a fluctuation relative to perfect Néel order.
Note that Eq. (6) cannot be an exact equality because it does not respect the fact that [s iσ , s † iσ ] = 1. It is more of a semiclassical approximation, valid in the limit of large spin S. Of course, the case S = 1/2 under consideration here is far from large, but all other analytical techniques of which we are aware for treating long-range order via slave particles (such as Bose-Einstein condensation of the Schwinger bosons [55,56]) have the same regime of validity. We refer for Ref. [54] for more details.
Inserting Eq. (6) into H tJ does not yield a solvable Hamiltonian on its own. Thus to progress further, we make a mean-field approximation by expanding H tJ to first order in where m is the Néel magnetization, we find that The second line is diagonalized by passing to momentum space and performing a Bogoliubov transformation, leading to the final Hamiltonian (neglecting constant terms) where the sum is over the 2D Brillouin zone (N is the number of lattice sites) and β q ≡ u q s q + v q s † −q is the transformed spinon operator. With where the frequencies ω q and vertices M kq entering into Eq. (8) are given by Normally one would determine m self-consistently from the ground state spinon occupation: This is known to give m ≈ 0.3 for a 2D square lattice [38]. We shall instead treat m as an independent parameter -this is an approximate (albeit crude) means of estimating exciton properties as a function of magnetization.

C. Self-consistent Born approximation
The self-consistent Born approximation (SCBA) gives the doublon-doublon and holon-holon propagators by the integral equation in Fig. 3. The same equation holds for each propagator separately. This approximation is expected to be accurate in the dilute-charge limit, where the charge dynamics is strongly affected by spinons but not vice-versa.
In terms of the doublon/holon self-energy Σ k ( ), Fig. 3 translates to (after a frequency integration) The quasiparticle spectrum k is given by the solution to Σ k ( k ) = k . Eq. (14) can be solved quite efficiently. Note that all ω q are positive [57], thus Eq. (14) in fact expresses Σ k ( ) in terms of the self-energy at lower frequencies.
We start at sufficiently negative , below which we approximate Σ k ( ) ≈ (1 + 2m)t 2 /2N q M 2 kq /( − ω q ), and then compute the self-energy at incrementally higher frequencies in terms of the previous values. To help avoid numerical errors, we add a small imaginary part (namely 0.1i(1 + 2m)J) to .
While one could proceed using the full Σ k ( ), it has been found that the quasiparticle dispersion can be wellapproximated by the form [58] This expression has a clear physical interpretation: t 1 is the amplitude for performing a two-step hop along the diagonals of the lattice, and t 2 is the amplitude for a two-step hop along the principal axes (see Fig. 4). Thus in what follows, we shall use for the single-particle propagators the simpler expression

D. Bethe-Salpeter equation
We next consider the two-particle Green's function (T denotes time ordering) and its Fourier transform G k d ke;k d k e ( d , e ; d , e ). Due to translational invariance, G depends only on differences in position and time, which we choose to parametrize by the relative coordinates with relative times defined analogously. The corresponding momenta are We will use absolute and relative momenta interchangeably, depending on notational convenience, with = + FIG. 5. The integral equation which determines the twoparticle Green's function, within the ladder approximation. The hatched square is the Green's function, and all other symbols are as in Fig. 3.
Eq. (19) always giving the relationship between the two. G k d ke;k d k e ( d , e ; d , e ) will often be written as G kk ;K ( , ; E).
Within the ladder approximation, G is determined by the integral equation of Fig. 5. Written out, where G k ( ) is given by Eq. (16) and the vertices M kq are as in Eq. (8).
Since our goal is to identify bound states, we reduce Eq. (20) to the Bethe-Salpeter equation. The details of this approach can be found in Ref. [59]. We assume that G has an isolated pole in the total energy E, near which it has the form where the "wavefunction" ψ k ( ), its time-reversed partner ψ k ( ), and the bound state energy E b remain to be determined. Inserting this ansatz into both sides of Eq. (20) and equating the residues at E b on each side, we obtain a non-linear eigenvalue problem (the Bethe-Salpeter equation): (22) E b and ψ k ( ) are given by the solution to Eq. (22). Note that they will depend on the center-of-mass momentum K. ψ k ( ) is the bound state wavefunction in a quite literal sense: it is the Fourier transform of where |b denotes the bound state and |0 denotes the ground state. Note that t = 0 is of particular interest, since it gives the amplitude for simultaneously observing the holon at site 0 and the doublon at site j. Thus to simplify the problem, we integrate Eq. (22) over , and furthermore, make the ansatz with Ψ k independent of . The explicit factor of E b − k d − ke is included so that Ψ k is the equal-time wavefunction, i.e., Ψ k = ψ k (t = 0). This ansatz allows us to perform the integral straightforwardly, giving a closed equation for E b and Ψ k : Eq. (25) is the two-particle Schrodinger equation albeit with an energy-dependent potential. We find the values of E b at which it has a non-zero solution, and record the corresponding eigenvector.
Strictly speaking, Eq. (24) is not a valid ansatz for ψ k ( ), i.e., it does not solve the frequency-dependent Eq. (22). However, it has a clear physical interpretation. The Fourier transform ψ k (t) gives the wavefunction for inserting the doublon and holon separated by time t (see Eq. (23)). The poles coming from the single-particle π π k x π π k y k /t propagators in Eq. (24) correspond to the phase factor acquired by the remaining particle during that interval, and our ansatz amounts to neglecting any other time dependence. This approximation has been applied previously to study holon-holon binding [42], and we expect it to be qualitatively accurate for our purposes.
Eq. (25) and those preceding it differ from the equations for holon-holon binding in two respects. First, the holon-holon equations must include exchange terms not found here. Second, due to the relative phase between the doublon-spinon and holon-spinon vertices, the effective potential in Eq. (25) would have the opposite sign for the holon-holon problem.

A. Single-particle properties
We first review the behavior of individual quasiparticles, determined within the SCBA as described above. Although these calculations have been reported previously, e.g., in Refs. [52,58], it will be useful to reproduce them here. Fig. 6 shows the quasiparticle dispersion throughout the Brillouin zone, with the magnetization set to the equilibrium value for concreteness. As noted above, it can be well-approximated by a next-nearest-neighbor hopping model with amplitude t 1 for moving along the diagonals of the lattice and amplitude t 2 for moving along the principal axes (see Fig. 4). The form of the dispersion is not sensitive to the value of magnetization. However, the effective hopping amplitudes, which we determine empirically by fitting the computed spectrum to Eq. (15), do depend on m as shown in Fig. 7.
Some features of the dispersion can be explained by a simple Hartree-Fock approximation to the original Hamiltonian, in which the Hubbard interaction is replaced by n i↑ n i↓ + n i↑ n i↓ . Assuming Néel order for n iσ , the Hamiltonian becomes a tight-binding model on a bipartite lattice with dispersion using that t U . Up to a constant shift, the second line is equivalent to Eq. (15) for the special case t 1 = 2t 2 . Note in particular that t 1 , t 2 < 0. Thus Hartree-Fock correctly predicts that the band minimum is within the lines k k + k y = ±π. It also correctly suggests that the bandwidth should be significantly reduced to O(J) instead of O(t). However, it incorrectly claims that the dispersion is degenerate along the entire magnetic Brillouin zone boundary. The more sophisticated SCBA resolves this degeneracy, identifying four minima at (k x , k y ) = (±π/2, ±π/2).
Recent work on magnetic polarons in the t-J model [60,61] has made clear that this behavior can be understood by the charge excitations forming spinon-charge bound states (the "polarons"), held together by strings of displaced spins (much as we have sketched in Fig. 1 but with individual charges). The charges are forced to move on the timescale set by their slower spinon partners, namely O(1/J), and are subject to the bipartite lattice felt by the spinons. The dispersion results shown above confirm this string picture quite nicely if one interprets them as being for the polaron as a whole. With that in mind, it is rather striking that these three approaches -SCBA, Hartree-Fock, and the string picture -all lead to consistent conclusions.
Returning to Fig. 7, we see that the bandwidth W increases noticeably as the magnetization increases. Equivalently, the single-particle mass decreases. Within the framework of our calculation, the explanation is clear: a doublon/holon can move only if a spinon takes its place (see Eq. (5)), and since one spinon factor is always in the direction of the Néel magnetization, the doublon/holon hopping term is proportional to (1 + 2m)/2.

B. Exciton properties
We now turn to the exciton properties as functions of magnetization, using the Bethe-Salpeter equation. All of the quantities presented here are straightforward to compute from the energy E b and wavefunction ψ k given by Eq. (22). Fig. 8 shows the energy of the lowest internal state −π −π/2 π/2 π k e, x −π −π/2 π/2 π k e, y Ψ k FIG. 9. Exciton wavefunction in (relative) momentum space at total momentum (Kx, Ky) = (π, π), for the same system as in Fig. 8.
as a function of the center-of-mass momentum K. As was the case for the single-particle dispersion, the shape of the exciton dispersion is not particularly sensitive to the magnetization. Note that the bottom of the band is not at the origin but rather at (K x , K y ) = (π, π). The wavefunction of the (π, π) state is shown in Fig. 9. We see that it has p-wave symmetry, unlike what one would expect for a two-holon bound state.
The binding energy, mass, and radius of the exciton are plotted versus magnetization in Fig. 10. We see that as one increases the magnetization m, the mass decreases while the binding energy and size increase. It is interesting to compare these trends with what one would expect for a conventional exciton formed via Coulomb attraction. In that situation, a decrease in mass is associated with an increase in radius and a decrease in binding energy. Here, we find a similar relationship between radius and mass, but the binding energy instead scales inversely with mass.
Eq. (25) can be simplified further in the large-t/J limit. We will see that E b scales as t, whereas k and ω q are asymptotically smaller [58]. Thus we can neglect the single-particle and spinon dispersions, leaving the equation the importance of this retardation, we have compared the results in Fig. 10 to what would be obtained through the static approximation (setting ω = 0 in the kernel of Eq. (22)). The static approximation would predict significantly different results, as seen in Fig. 11: the binding energy would instead decrease slightly with magnetization and the mass would increase slightly. Thus the retardation of the effective interaction is an essential ingredient to the behavior seen here.
Similarly, one can ask whether the trends observed in Fig. 10 are due primarily to changes in the spinon behavior or rather due to the single-particle mass, which itself decreases with magnetization. We have repeated the above calculations under "fixed-mass" conditions, in which the single-particle parameters t 1 and t 2 are kept fixed (to their values at m = 0.5) as we vary the magnetization. Fig. 11 shows that each of the three observables responds differently. The binding energy becomes more sensitive to magnetization, indicating that the quasiparticle and spinon properties play antagonistic roles. On the other hand, the exciton mass becomes less sensitive -the change to the effective interaction suppresses the mass by itself. Finally, the exciton radius shows the reverse behavior to before, instead decreasing with magnetization (although the size remains quite small in absolute terms). The recent pump-probe experiments in Refs. [22,23] have investigated how excitons are influenced by magnetic order in the Mott insulator Na 2 IrO 3 . Our results support their interpretation in some aspects but not in others. In Ref. [22], the authors observe an increase in the fraction of bound excitations when below the Néel temperature, which they attribute to an increase in the exciton binding energy. Fig. 10 shows that magnetic order does indeed increase the binding energy. On the other hand, Ref. [23] demonstrates that the relaxational dynamics following a pump are slower in the presence of magnetic order. This is attributed to the mass increasing with magnetization, yet we have observed the opposite (consistent with past works calculating the dependence on J/t [53,58,62]). Given the highly non-equilibrium nature of the experiments, as well as the approximations inherent in an analytical approach, further investigation is clearly needed.
Finally, let us compare the present calculation of doublon-holon binding to that of holon-holon binding, which is obviously of significant interest in its own right [51,63]. Clearly the two have much in common, yet there are two important differences. First, the integral equation which determines the two-particle Green's function (Fig. 5) has an additional exchange term due to the indistinguishability of the holons. Second, even the direct term comes with an extra minus sign, i.e., the effective interaction is of opposite sign. The sign can be removed by redefining the hole operator on one sublattice, but the additional phase may modify further results depending on the application. It is important to keep these distinctions in mind when relating the present results to the high-T c literature.

IV. CONCLUSION
We have studied the role that magnetic order plays in the formation of excitons within Mott insulators, using the Hubbard model as a concrete Hamiltonian. The binding energy increases in the presence of (antiferromagnetic) magnetization, whereas the exciton mass decreases. The size of the exciton increases slightly, yet the radius is never more than a lattice spacing. Using the standard classification, these are Frenkel excitons regardless of magnetic order.
In addition, we have established that the trends observed here require a detailed understanding of the manybody dynamics in these systems. Retardation effects in the effective spinon-mediated interaction are essential. Furthermore, the constituent charge and spin excitations are each affected separately by the background magnetic order, in ways cooperative for some exciton properties but antagonistic for others.
It must be noted that despite the complexity, there are significant limitations to our approach. In particular, we have made approximations in the spirit of linear spin-wave theory, which is only justified at large spin S and full magnetization (neither of which we assume here). Thus we do not expect these results to be quantitatively accurate -we instead view this analysis as expressing our physical intuition regarding Mott excitons in the language of slave particles, from which we can make sharp predictions to be verified or falsified by more systematic investigations.
As an outlook, the predictions made here will be important when analyzing recent and future experiments on the optical properties of strongly correlated electronic materials. The existing experiments are quite complex, and require interpretations of their own. Our results agree with those interpretations in some respects but disagree in others. A complete understanding of the systems will require numerous approaches, both experimental and theoretical, including but not limited to the one described here.
Particularly promising are the recent experiments on fermionic atoms in optical lattices [34,35,37]. Since ultracold gases do not have many of the complicating features found in condensed matter systems, we expect that this will be a valuable direction to explore further. It is also likely that our conclusions, being based on the single-band nearest-neighbor Hubbard model, are more applicable to those systems than to materials such as the iridates. Importantly, current quantum gas microscopes allow one to directly create localized doublons and holons via optical tweezers, and reliably measure the spin correlation functions [36]. Such an unprecedented direct access to the system microscopics will provide a powerful way of investigating many-body excitons.

V. ACKNOWLEDGEMENTS
We would like to thank Eugene Demler, Immanuel Bloch, Andrew Allocca, Zachary Raines, Yang-Zhi Chou, Fabian Grudst, and Annabelle Bohrdt for stimulating discussions and suggestions. This work was supported by the NSF Physics Frontier Center at the Joint Quantum Institute, the NSF DMR-1613029 and US-ARO contracts W911NF1310172 (T.