Optical signatures of antiferromagnetic ordering of fermionic atoms in an optical lattice

We show how off-resonant light scattering can provide quantitative information on antiferromagnetic ordering of a two-species fermionic atomic gas in a tightly-confined two-dimensional optical lattice. We analyze the emerging magnetic ordering of atoms in the mean-field and in random phase approximations and show how the many-body static and dynamic correlations, evaluated in the standard Feynman-Dyson perturbation series, can be detected in the scattered light signal. The staggered magnetization reveals itself in the magnetic Bragg peaks of the individual spin components. These magnetic peaks, however, can be considerably suppressed in the absence of a true long-range antiferromagnetic order. The light scattered outside the diffraction orders can be collected by a lens a with highly improved signal-to-shot-noise ratio when the diffraction maxima are blocked. The collective and single-particle excitations are identified in the spectrum of the scattered light. We find that the spin-conserving and spin-exchanging atomic transitions convey information on density, longitudinal spin, and transverse spin correlations. The different correlations and scattering processes exhibit characteristic angular distribution profiles for the scattered light and, e.g., the diagnostic signal of transverse spin correlations could be separated from the signal by the scattering direction, frequency, or polarization. We also analyze the detection accuracy by estimating the number of required measurements, constrained by the heating rate that is determined by inelastic light scattering events. The imaging technique could be extended to the two-species fermionic states in other regions of the phase diagram where the ground state properties are still not fully understood.


I. INTRODUCTION
The two-species fermionic Hubbard model is one of the most studied models with strong correlation effects in condensed matter physics, particularly since Anderson proposed in 1987 [1] that the Hubbard model (or its strong coupling limit, the t-J model) is the minimum model to describe the physics of the high temperature superconductors [2]. It is well understood that at halffilling, i.e., at an average of one spin-1/2 fermion per site, the Hubbard model has the ground state which is a Mott insulator with antiferromagnetic (AFM) ordering. Away from half-filling, quantitatively accurate predictions are scarce, but it is generally believed that a d-wave superconductor can be present in some regions of the phase diagram, hence the relevance of the Hubbard model for the high temperature superconductors.
With the rapid advance in cooling and trapping technologies of neutral atoms in optical lattices, there is the exciting prospect that the phase diagram of the fermionic Hubbard model can be explored via quantum simulation (rather than numerical simulation) in an optical lattice set up with ultracold two-species fermionic atomic gas [3][4][5]. Such a system has been cooled down in a laboratory to the Mott insulator regime [6,7] and recent experiments have demonstrated evidence of AFM correlations [8,9], * francisco.cordobesaguilar.2009@live.rhul.ac.uk owing to exchange coupling between the atomic spin states. Magnetic ordering has also been observed in doublon formation in tilted lattice systems [10,11]. Accurate diagnostic tools for ultracold atoms in lattices form key elements in emulating strongly-correlated physics of condensed-matter systems. Currently, advanced imaging provides a microscopic scanning technology of atoms with a single-site resolution [12][13][14][15][16][17][18][19] in which case each atom may resonantly scatter thousands of photons while being detected. The analogy between x-ray (or neutron) diffraction of crystalline structure in solids and offresonant light scattering from an ultracold atomic gas in an optical lattice provides an alternative route. The periodic crystalline order leads to constructive interference and the emergence of sharp diffraction maxima, or Bragg peaks, revealing the underlying lattice structure. Importantly, diffuse background of scattered light in between the diffraction maxima can convey information on fluctuations of atomic positions. The far-field off-resonant imaging therefore provides a powerful probe of stronglycorrelated ultracold atom systems in optical lattices, being sensitive, e.g., to correlations, excitations, and temperature.
Spontaneously scattered off-resonant light was first proposed as a diagnostic tool for correlations of ultracold atoms in an optical lattice in Ref. [20]. This study considered a two-species gas where the hopping element of the atoms between adjacent sites acquired an artificially constructed phase factor from the Peierls substitution. It was then shown that the resulting topological properties of edge states, which originate from nontrivial spin correlations, can be mapped onto fluctuations of scattered light and directly detected [20,21]. Several other studies have addressed signatures of atom statistics in an optical lattice from the scattered light with [22,23] and without [24][25][26][27][28][29][30][31][32][33] an additional optical cavity.
A two-species fermionic Hubbard model exhibits a Mott insulator state where the fluctuations of the total on-site atom number are suppressed [6,7]. At sufficiently low temperatures the effective exchange coupling can generate magnetic ordering of the atomic spins within the Mott state. At half-filling with equal spin populations this leads to alternating checkerboard pattern of atom densities of individual spin components. The resulting doubling of periodicity of atom densities leads to the emergence of additional magnetic Bragg peaks of scattered light that have been proposed as a detection mechanism of AFM ordering [34]. The extra Bragg peaks in the optical signal of a two-species system were experimentally observed using an artificially constructed stripe pattern of atomic densities [35].
Here we analyze the far-field diffraction pattern for offresonantly scattered light from a two-species fermionic atomic gas in an optical lattice. We consider a tightlyconfined two-dimensional (2D) Hubbard model at halffilling with equal spin populations. We show how optical diagnostics can provide quantitative information on properties of strongly correlated states in a lattice. The time-ordered correlation functions that can be calculated using the Feynman-Dyson perturbation series are mapped onto the fluctuations of the scattered light and detected in the optical signal. We evaluate the relevant correlation functions in the mean-field theory (MFT) and in the random phase approximation (RPA) for the AFM ordered state. MFT provides information on singleparticle excitations but fails to capture the effect of quantum fluctuations and collective excitations that are approximately incorporated in RPA. In the RPA calculations we follow the approach of Ref. [36].
The scattered light intensity may be separated into elastic and inelastic components. In the elastic scattering process the atom scatters back to its original state. The elastic part produces a diffraction pattern from a nonfluctuating atom density, analogous to that of N s × N s diffracting apertures when the lattice has N s × N s sites. The overall envelope of the diffraction pattern is determined by the lattice site wave function of the atoms. The magnetic ordering specifies the relative atom densities of the two spin components and appears as additional Bragg peaks. The measurement of the magnetic Bragg peaks is possible, provided that the lattice spacing a > λ/ √ 2, where λ denotes the probe wavelength. Our analysis confirms that detecting the magnetic peaks by scattering light from a single spin component alone constitutes the most accurate observable for the staggered magnetization of the AFM order. However, the absence of true long-range AFM order can significantly suppress the magnetic Bragg peaks. We demonstrate this by considering a phenomenological Ising model for the staggered magnetization when the system exhibits a finite correlation length.
The inelastic scattering processes, on the other hand, are those in which an atom scatters between two different quasimomentum states. The inelastically scattered light conveys information on correlations between the atoms and results in diffuse scattering of light outside the diffraction orders, generating fluctuations of the diffraction pattern. We collect the light in the near-forward direction with a lens and block the diffraction maxima, so that the amount of elastically scattered light entering the detector is suppressed [24]. This method is then extended to detection of light in the direction approximately perpendicular to the propagation direction of the incident field. In Ref. [24] detecting the near-forward light was shown to provide an experimentally feasible technique for measurements of temperature of fermionic atoms in a lattice. We find that in an AFM ordered two-species state the scattered light in the near-forward direction is not only sensitive to the temperature of the atoms but provides also a suitable probe of density and longitudinal spin correlations (from spin conserving atomic transitions). The scattering processes in which the atomic spin state changes, on the other hand, are prominent in the scattered light perpendicular to the light propagation direction and can be employed in detection of transverse spin correlations. Furthermore, we estimate the detection accuracy of the magnetic ordering in different measurement configurations. This is done by calculating the number of required experimental realizations of the lattice system to detect the order parameter above the shot noise of light at a desired accuracy. In each experimental realization the scattered light heats up the atomic sample and the total number of inelastic scattering events is constrained to be a small fraction of the total atom number. We find that a strong lattice and trap confinement is beneficial for the measurements in suppressing scattering of atoms to higher energy bands as compared with the inelastic lower energy band scattering.
We also calculate the spectrum of the scattered light and show how it reveals the excitation of the system that could be measured using optical heterodyne techniques [37]. The differences between the MFT and RPA treatments are especially prominent in the spectrum of transverse spin correlations: The low-energy collective mode excitations manifest themselves as a well-separated peak (at sufficiently strong interaction energy) from the gapped single-particle excitations.
The optical setup that blocks the Bragg diffraction maxima in the measured signal can considerably reduce the scattered light that is insensitive to correlations, therefore improving the signal-to-noise efficiency. Similar separation of less important part of the signal, for instance, in atom shot-noise correlation measurements [38,39] would be challenging. We estimate the detection accuracy of the AFM ordering in the scattered light intensity by calculating the number of experiments needed to achieve a given measurement accuracy of the magnetic order parameter, when the heating rate of the atoms limits the total number of possible scattering events.
The proposed technique differs from the simple angleresolved imaging of the diffraction peaks, since we collect the inelastically scattered light with a lens that has a relatively large numerical aperture (NA). By analyzing the spectrum of the scattered light using the realistic optics setup, we show that the essential features of the excitation spectrum can be captured even when the light is collected over a range of scattering angles. A large lens provides a stronger signal even when the number of inelastically scattered photons remains a small fraction of the total atom number. In less demanding detection scenarios (e.g., in large lattices and at high temperatures) even a single experimental realization of the lattice system can then be sufficient to determine, e.g., the approximate temperature of the atoms in the lattice. This contrasts with the single-site microscopy [12,40] in which case typically several thousands of photons are scattered from each atom and the light is also simultaneously used to cool down the atoms. The microscopy also requires a scanning of the lattice sites that in very large lattice systems may become less practical. The off-resonant light scattering, on the other hand, becomes more efficient a method when the size of the system increases. Other advantages of the off-resonant imaging include the possibility for spectral measurements of excitations and the access to correlation functions that include combinations of different spin states via spin-exchanging scattering processes. Finally, the diffractive far-field imaging does not need to be limited to probing the atoms only by light, but also matter-wave probes are possible [41].
The remainder of the paper is organized as follows: Section II presents a summary of the key results. We then start in Sec. III by introducing the basic formalism of the lattice system. We continue in Sec. IV, where the MFT and RPA results for the lattice system are presented. The scattered light as a diagnostic tool is introduced in Sec. V. In Sec. VI we present results for the scattered light intensity for 40 K atoms. The specific experimental setups for the optical detection and the estimates for the measurement accuracy are considered in Sec. VII. The scattered spectrum is studied in Sec. VIII and some concluding remarks are made in Sec. IX. Finally, a diagrammatic description of the RPA susceptibilities is presented in Appendix A and the finite temperature MFT susceptibilities are given in Appendix B.

II. SUMMARY OF KEY RESULTS
In this Section we briefly highlight the main findings of the paper. We study a two-species fermionic atomic gas trapped in a tightly-confined 2D optical lattice [Eq. (9)]. We assume equal spin populations and that on average, there is one atom per site in the lattice. The atoms are probed by incoming laser light propagating perpendicular to the lattice. The scattered light is collected by a lens positioned at various angles relative to the propagation direction, see Fig. 6.
We study the Mott insulator state of the Hubbard model where fluctuations of total on-site atom number are suppressed. We focus on the AFM ground state, characterized by a checkerboard pattern of ordering (period doubling). The associated AFM order parameter is the staggered magnetization m, defined as [Eq. (19)] is the real space spin operator. The magnetic ordering wavevector is Q = (π/a, π/a) (a denotes the lattice spacing) and r j is the coordinate of the center of the lattice site j. We employ MFT in Sec. IV B to diagonalize the Hamiltonian around the AFM order parameter m using the Bogoliubov transformation. This leads to the gap equation and description of single-particle excitations of the system. However, MFT does not include quantum fluctuations (specifically, spin waves). Crucially, such quantum fluctuations modify qualitatively the correlation functions of spin and density operators, which then leads to large changes to the optical response (see Secs. VI and VIII). Hence, in Sec. IV D and Appendix A, we describe the calculation of the correlation functions by a partial summation of the Feynman-Dyson perturbation series in RPA that incorporates the correct spin wave physics at strong coupling.
We show in Sec. V A how quantum statistical correlations of the atoms can be mapped onto fluctuations of the scattered light. Measurements on the scattered light therefore convey information about the correlated phases of the ultracold atoms in the lattice. The formalism describing the relationship between optical signal (the intensity and spectrum of the scattered light) and the atomic correlations is described in Sec. V: see Eqs. (66)-(68) for the scattered intensity, and Eqs. (78)-(80) for the scattered spectrum.
The scattered light intensity contains contributions from the elastic and inelastic scattering events. The elastic part generates the diffraction pattern of the atomic lattice structure. If the detected signal cannot distinguish the two spin components, the light provides almost no information about the AFM order. But with spinspecific imaging, the emerging AFM order and the period doubling can be identified as additional Bragg peaks.
Specifically, we find that the elastic component of the scattered light intensity is given by In this expression the dependence on level structure, polarization, and scattering direction are coded in M g3g4 g2g1 [Eq. (64)], which contains the dipole matrix element [see Eqs. (55) and (56)]. ∆k is the change in the momentum of the scattered photons, Eq. (54), while∆k is the same momentum projected into the lattice plane. Here the Debye-Waller factor α ∆k depends on the lattice site wave function and determines the overall envelope of the diffraction pattern, see Eqs. (63) and (65). The density operator in momentum space for species g iŝ whereĉ k,g is the annihilation operator at momentum k for hyperfine state g.
For the AFM state this can be evaluated to be where f g = 1/2 is the atomic filling factor of species g at half-filling and η(↑) = 1, η(↓) = −1 [Eq. (22)]. Here |u∆ k | 2 [Eq. (69)] generates the diffraction pattern from the periodic atom density [Eq. (76)]. The second term of Eq. (3), when substituted into Eq. (1), is responsible for these new peaks due to the period doubling in the AFM state. This effect was first analyzed in Ref. [34].
In the absence of true long-range AFM order the staggered magnetization, however, may vary in space. We study the effects of short-range AFM order on the optical response by introducing a phenomenological Ising model that demonstrates a significant suppression of the magnetic Bragg peaks when the correlation length is not much larger than the size of the lattice.
While elastically scattered light gives information about the AFM ordering of the system, the inelastic component directly probes the quantum atomic correlations of the spin operators (we follow the notation convention from Ref. [36]) or the total density operator We explicitly show how the time-ordered many-body correlation functions obtained in RPA [Eqs. (45)- (49)] are mapped onto the experimental observables of the optical signal [see, e.g., Eq. (35)]. We then find that the inelastically scattered light intensity for the two-component 40 K gas is given by Eq. (97), where the equal time 2 × 2 matrix response function [cf.
Eqs. (94) and (95)] O i q can be any of the operators in Eqs. (4), (5) and (6). The subscript c in the expectation value denotes a connected correlation function for which q = 0. We have also introduced the two-component vector [Eq. (96)] u u u † k = (u * k , u * k+Q ). The matrix notation reflects the fact that in the AFM ordered state with period doubling, the original Brillouin Zone is split up into the RBZ, and the zone outside the RBZ which is connected by Q to the RBZ [ Fig. 1].
In Eq. (7), the scattering contributions in which spin is conserved are proportional to density and longitudinal spin correlations. The spin-exchanging transitions (see Fig. 7) generate the term depending on the transverse spin correlations. The two processes exhibit very different angular distribution of the scattered light, as shown in Eq. (93), and can also be separated in frequency. Thus, another key result of this paper is that specific quantum correlations can be separated in the detected signal.
For the detection of inelastically scattered light we consider two different experimental configurations: One has the scattered light collected by a lens in the near-forward direction with the central diffraction peak blocked, see Fig. 16(b), and the other one has light collected in the perpendicular direction, Fig. 16(c). The full intensity and spectrum are then evaluated numerically by integrating over the scattering directions (the momenta) ∆k collected by the lens.
We estimate the number of experimental measurements needed to obtain a given relative accuracy in the determination of the magnetization or temperature [24]. This minimum number of measurements is determined by two factors. On one hand, photon shot noise dictates a minimum difference in the number of photons detected for two close AFM magnetization values. On the other hand, the total number of photons collected is set by the experimental duration. This duration is constrained by the heating rate of the system due to inelastic scattering of atoms by light. In order to satisfy these two constraints, a minimum number of experimental realizations has to be achieved, see Sec. VII B and Figs. 17-20. Armed with these measurement accuracy data, we can propose near-optimal parameters for the experimental configurations mentioned. The best magnetization measurement accuracy can be achieved with spin-specific detection where the photons are collected near the perpendicular direction or around the direction of the magnetic Bragg peak, while the temperature of the atoms can be measured in the near-forward direction [ Fig. 16(b)]. In the perpendicular direction [ Fig. 16(c) and 20], it is preferable to detect only the light scattered from spinexchanging transitions. The light scattered by spinconserving or spin-exchanging transitions can be separated owing to the different frequency and polarization (Sec. VII C) of the scattered photons.
We calculate the spectrum of the scattered light in Sec. VIII. We obtain an expression similar to the scattered light intensity, with the static correlations replaced by dynamic ones S S S ij (q, ω) [Eqs. (116) and (35) and Sec. VIII]. The spectrum is calculated for a realistic optical setup where the light is collected using a finiteaperture lens, Fig. 23. We find that at moderately large to large U/J, the collective excitations generate a sharp peak in the low energy part of the spectrum, which is separated by an energy gap of order U from a more broad feature generated by single particle excitations. The location and width of these features can give a quantitative measure of the AFM order.

III. OPTICAL LATTICES AND THE HUBBARD MODEL
Optical lattices, generated by a standing wave (periodic) laser potential, provide ideal tunable systems where almost every system parameter can be changed independently. In typical experimental situations, the trapping potential for the atoms is a superposition of the lattice and an external, approximately harmonic, trap. For the case of sufficiently weak external potential, the harmonic trap may be ignored. Experimentally, it is also possible to produce entirely homogeneous lattices; a first step towards this has recently been demonstrated for a Bose-Einstein condensate in a uniform trap [42]. Here we will ignore, for simplicity, any modulations of the uniform lattice potential. The effects of the additional harmonic potential, e.g., in the context of the Mott insulator states has already been addressed by several studies [6,7]. For a typical 2D square lattice, the periodic potential then reads (r = (x, y, z)) V (r) = s x E R sin 2 πx a + s y E R sin 2 πy a + 1 2 mω 2 z z 2 . (9) In our case we choose the lattice depth s = s x = s y , similar to the 2D lattice experiments with a disk-like lattice [12]. Here the frequency of the harmonic confinement in the z direction is denoted by ω z and the lattice light recoil energy, E R , is defined by We define an effective wavenumber for the optical lattice in terms of the lattice spacing a by When the effective wavenumber coincides with the laser wavenumber, we have k l = 2π/λ l and a = λ l /2, where λ l denotes the wavelength of the laser beam generating the lattice. In the case of accordion lattices the lattice spacing is manipulated by optical components and can be considerably increased [43][44][45].
The Hamiltonian that describes the atoms in the optical lattice can be written in terms of second quantization operators. These operators are related to the field operatorsΨ g (r) for the hyperfine state g, via the expansion in the complete set of Wannier functions w n,j (r) = w n (r − r j ) representing a localized state at position r j and band n [46]Ψ g (r) = n,j w n,j (r)ĉ njg , whereĉ njg is the fermionic annihilation operator for hyperfine state g at site j = (j x , j y ). The spatial variation of the density along the lattice is represented by the modulation of the expectation values ĉ † njgĉ njg . In this study we assumed that the potential in the xy plane can be described by the uniform lattice potential alone and any deviations, e.g., due to harmonic trap may be neglected. The spatial variation of Eq. (12) between different lattice sites is therefore entirely encapsulated in the phase factors ofĉ njg .
We consider only the regime where U, T are much smaller than the bandgap, so that the higher bands in the initial ground state are not populated. Hence, only the n = 0 state in Eq. (12) is included in the analysis of the ground-state properties, and we drop the band index for now. (Later, we will consider the effect of light scattering of atoms into higher bands in Sec. V C.) In order to describe the lattice system we introduce the one-band Hubbard Hamiltonian [47] Here µ is the chemical potential and the hopping amplitude J only includes tunneling of atoms between the nearest-neighbor sites as indicated by j 1 j 2 in the summation. The on-site interaction strength is denoted by U . It can be modified by changing the spatial confinement or via a Feshbach resonance of the s-wave scattering length a s [48]. In this paper we consider the ground state of the system consisting of atoms trapped in two different hyperfine states, for which we use the "spin" labels | ↓ and | ↑ , in analogy to electrons in solids. For deep lattices, the lattice potential can be approximated locally close to the site minimum as a harmonic potential with frequencies Then, the ground state Wannier function becomes where l i is the oscillator length. With this approximation, In the limit s ≫ 1, the hopping matrix element can be obtained from the 1D Mathieu equation as Thus J and U are readily changed by tuning the laser intensity and magnetic fields in the experiment.

IV. APPROXIMATE SOLUTION OF THE HUBBARD MODEL AT HALF-FILLING
A. Low energy physics of the Hubbard model In the experimentally relevant large U/J limit, and at half-filling (i.e., on average one atom per site), the ground state of the Hubbard model is a Mott insulator. This is a state where the energy cost of a doubly occupied site (∼ U ) is too large compared to both the hopping energy J and the temperature T when k B T U . In the Mott insulator state the on-site total atom number fluctuations are suppressed. Immediately below the onset of the Mott transition there is, however, very little energy cost in mixing the relative populations of the two spin states and the on-site relative atom number may still fluctuate. At even lower temperatures, below a new characteristic scale determined by the Néel temperature, the spins in a square lattice form into an AFM pattern. Classically, with equal numbers of ↑ and ↓ atoms, this represents a checkerboard pattern where spin ↑ and spin ↓ atoms occupy alternating sites. In the alternating density pattern, virtual hopping of atoms between nearest neighbor sites becomes energetically favorable. Second-order perturbation theory at large U/J then leads to an effective spin exchange interaction ∼ 4J 2 /U between the atoms in the neighboring sites, defining the Néel temperature scale. At such low energies, the system is described by the (spin only) Heisenberg model [49].
At zero temperature, it is known that the ground state of the 2D Hubbard model (or the effective Heisenberg model) has true long-range AFM order; for a review, see [50]. The order parameter for the AFM state is the staggered magnetization, defined as where the spin operator in real space isŜ z j = (ĉ † j↑ĉ j↑ − c † j↓ĉ j↓ ) [cf. the momentum space version in Eq. (4)]. For the checkerboard pattern in a 2D square lattice, the ordering momentum is This in turn corresponds to the expectation value of the number operator for each spin component at half-filling where The AFM state breaks lattice translation symmetry such that the new unit cell in real space is doubled in size (to contain both spin components). Hence, in momentum space, the Brillouin zone is halved to become the RBZ. In Fig. 1, we show one choice of the RBZ, which is bounded by the four lines k x ± k y = ±π/a. The filled circles denote the momentum states belonging to the RBZ, and the ordering vector Q links a state within the RBZ to one outside (and vice versa). Note that in order to avoid double counting, half the states on the bounding lines belong to the RBZ and the other half are outside of the RBZ.
The AFM order can approximately be characterized by a MFT, and results so obtained can be used to calculate the effect of the ordering on off-resonant light scattering. However, MFT contains only single particle excitations (of order U at large U/J). Importantly, MFT does not include quantum fluctuations that give rise to collective excitations at low energies. Hence MFT fails to capture the spin waves of the AFM state. Schrieffer et al. [36] and others [51][52][53] have demonstrated how the RPA can partially incorporate quantum fluctuations and describe the spin wave excitations of the system. Thus, we will use the RPA at zero temperature to study how collective modes modify the scattered light.
On the other hand, at any non-zero temperatures, there is no true long-range AFM order in the thermodynamic limit in the Hubbard model or the Heisenberg model [54]. This loss of long-range order is due to the enhanced quantum and thermal fluctuations in 2D. Instead, there is at most quasi-long-range order; for a review, see [55]. For example, in the Heisenberg model, the spin correlation function Ŝz iŜ z j decays with the distance r = |r i − r j | as exp(−r/ξ AFM ), for r >> ξ AFM . ξ AFM is the AFM correlation length in 2D given by where c 0 ∼ 0.26, b 0 ∼ 0.2, J H is the Heisenberg exchange coupling and a is the lattice spacing [49,56]. In contrast, true long-range order is signaled by an additional constant term in the spin correlation function, which is just m 2 . In the absence of true long-range order, the existence of the length scale ξ AFM leads to the following physical picture: If the system size L is such that L ≫ ξ AFM , the system is made up of small domains of size ∼ ξ AFM within which there exists AFM ordering. The order parameter m in different domains, however, are uncorrelated, so that there is no net AFM order overall. Such short-range order can be masked, if the system size is small compared to ξ AFM and only consists of a single domain (see [55] for the actual form for the spin correlation function at r ξ AFM ). We can estimate roughly the temperature at which such a finite size effect can become significant in ultracold atom lattice systems. The crucial ingredient is the strong exponential dependence of ξ AFM on T in Eq. (23). The Heisenberg exchange coupling J H can be related to the Hubbard model parameters by J H ≈ 4J 2 /U in the large U/J limit. At U/J = 6, we find that ξ AFM ∼ 1000a for k B T /J ∼ 0.1, or T /2J ∼ 0.05 where the Fermi temperature is ∼ 2J/k B at half-filling [2J is half the bandwith of the dispersion relation]. On the other hand, at T /2J ∼ 0.1 the correlation length ξ AFM ∼ 20a is already smaller than current typical optical lattice size of ∼ 30 sites in each dimension.

B. Mean-field Hamiltonian
In order to use the staggered magnetization as the MFT order parameter for the Hubbard Hamiltonian, we writen jg = n jg + (n jg − n jg ). It is then assumed that in the MFT Hamiltonian, terms second order in the fluctuation (n jg − n jg ) are small. Hence, the interaction in Eq. (13) can be rewritten as: U jn j↑nj↓ ≈ U j (n j↑ n j↓ +n j↓ n j↑ − n j↑ n j↓ ).
As usual, the real space Hamiltonian can be simplified by transferring to momentum space viâ The summation over momenta is defined for the whole Brillouin zone: (k x , k y ) = 2π Nsa (j x , j y ) with −N s /2 ≤ j x,y ≤ N s /2 − 1, where we assume N s sites in each dimension. If, for simplicity, we assume a translationally invariant system with periodic boundary conditions, we only need to consider coupling between the momentum states k and k + Q in the momentum space representation of the Hamiltonian. Then it is useful to explicitly split up the fermion operatorĉ kg defined over the whole Brillouin zone, to two operatorsĉ kg andĉ k+Qg , where k is now defined only for the RBZ. They are then collected into a two-component Nambu spinor (in analogy to superconductivity) Using the definition of the staggered order parameter [Eq. (21)] and substituting Eq. (24) and (25), the MFT Hamiltonian H can be written as a 2×2 matrix where we have defined the coefficient and the order parameter (also called gap parameter) Also, at half-filling we have used µ = U/2 to simplify the equation. Note that the Hamiltonian factors into separate spin sectors due to the choice of ordering in the z direction. The MFT Hamiltonian (26) can be diagonalized by a canonical (Bogoliubov) transformation with appropriately chosen v k,g and u k,g . The solution is The last equation holds because ∆ g = η(g)∆ [see Eq. (22)]. The subscript of the new operators in Eq. (29) refers to the energy bands in RBZ. In the new basis, the Hamiltonian reads (32) In summary, the MFT ansatz [Eq. (19)] leads to the MFT Hamiltonian Eq. (26) that couples a particle at k to a hole at k + Q. The Bogoliubov rotation then transforms the original one-band Hubbard model, defined over the full Brillouin zone, to a two-band model with energies ±E kg . To accommodate the same number of states, the momenta in the two-band model are defined over the RBZ only. At temperature T = 0, the ground state has all the negative energy modes filled at halffilling, i.e., band 1 is filled. At finite temperatures, the usual Fermi-Dirac distribution describes the occupation of the two bands. The MFT is found by minimizing the total energy with respect to m at a given temperature using the Hamiltonian (32). The resulting order parameter equation is At half-filling at T = 0 with U ≫ J the solution saturates towards m = 1/2. The MFT critical temperature T C,MFT can be obtained from solving Eq. (33) with m = 0. The value of staggered magnetization m is shown in the (U, T ) space in Fig. 2.
As mentioned at the beginning of this Section, in 2D, there is strictly no true long-range order except at T = 0 [54]. However, one can define instead a cross-over temperature T X below which there is at least short-range order, see Borejsza and Dupuis [57,58]. In the moderate to large U regime, T X is defined as the temperature at which the AFM correlation length equals the lattice spacing a: ξ AFM (T X ) = a. It turns out that [57] at large U/J, k B T X ∼ 4J 2 /U , the Néel temperature scale. On the other hand, at small U/J, T X ∼ T C,MFT where the critical temperature T C,MFT is exponentially small in J/U . A detailed analysis of the cross-over phase diagram can be found in Ref. [57]. Figure 3 (based on Fig. 2 of [58]) shows a schematic phase diagram for the Hubbard model in 2D at half-filling.

C. Mean-field susceptibilities at zero temperature
One key result of this paper is the explicit connection between the physical observables of inelastic scattered light intensity and spectrum (coded in static and dynamical structure factors), and susceptibilities calculated for the Hubbard model. We first indicate the chain of relations from structure factors to susceptibilities in Sec. IV C 1, and then proceed to give the MFT susceptibilities in Sec. IV C 2 and the RPA ones in Sec. IV D.

Time-ordered and retarded correlation functions
We shall show in Sec. V that the inelastic scattered intensity [Eq. On the other hand, anticipating the calculations of Sec. IV D, the RPA method involves a Feynman-Dyson perturbation series that requires the use of time-ordered correlation functions (often called susceptibilities). The Fourier transform in time of time-ordered correlation function is defined as where T represents the time ordered product, andÔ i q (t) can beρ q (t) orŜ i q (t). Fortunately, linear response theory [59] allows to connect time-ordered correlation functions to the response functions needed for scattered intensity and spectrum, as follows. First, at temperature T = 0, the dynamic response function S ij (q, q ′ ; ω) [Eq. (82)] is related to the retarded susceptibility χ ij R via where the indices i, j = ρ, z, +, −. Note the factor of 2 in Eq. (35) which is there to compensate for the unconventional factor of 2 in Eq. (34). The superscript R in the susceptibility χ ij R denotes retarded correlation functions that corresponds to the normally ordered physical observables of the scattered light intensity or spectrum. Next, these retarded susceptibilities can be analytically continued from time-ordered susceptibilities, resulting in [59] Re For the MFT (Sec. IV C 2), and for the RPA, (Sec. IV D), all susceptibilities are time-ordered (and Fourier transformed in time).
For future reference, we also point out that the static response function [Eq. (72)] can be obtained from the dynamic response function by integrating over ω,

Mean-field susceptibilities
Because of the RBZ structure of the AFM state, the MFT susceptibilities defined in Eq. (34) can also be written in a 2×2 matrix form, .
The subscript (0) in all the susceptibilities signifies that these correlation functions are calculated within MFT. A similar 2×2 structure can also be written for the RPA susceptibilities. We shall use a boldface χ χ χ to denote the susceptibility matrix. Here, and for the rest of the Section, unless it is explicitly indicated that this is not the case, we use the notation that the momentum q belongs to the RBZ only. Schrieffer et al. [36] and others [51] have calculated the zero temperature susceptibilities using MFT. Here we will present the results; see Appendix A for an outline of the derivation.
For the 2D Hubbard model in a square lattice at halffilling in the AFM ground state, the MFT density susceptibility is Here in Eq. (40) the momentum q belongs to the full BZ for both sides of the equation. We have explicitly shown the convergence factor +iδ (δ > 0) appropriate for timeordered susceptibilities. For the calculations in this paper we set δ a small but finite value which determines the frequency resolution in calculated spectra. Since the total density does not distinguish between the spin components, there can be no component in the susceptibility matrix which transfers momentum Q between atoms, hence the diagonal nature of the matrix in Eq. (39). The form of the susceptibility in Eq. (40) is similar in structure to susceptibilities for BCS superfluidity. It turns out that at the MFT level, the longitudinal spin susceptibility is equal to the density susceptibility This however, is no longer true when RPA is used to compute the susceptibility, cf. Eqs. (45) and (46).
The transverse spin susceptibility, on the other hand, is sensitive to the coupling of atoms that differ in momentum by Q. This then leads to two distinct components in the transverse spin susceptibility matrix with nonzero off-diagonal components Similarly to Eq. (40) the momentum label q in Eqs. (43) and (44) is valid for the full BZ on both sides of the equation and χ +− (0) (q, q+ Q; ω) = χ +− (0) (q+ Q, q; ω). Note that the transverse spin susceptibilities are different to the longitudinal one. This is due to isotropy in spin space being broken in our MFT: we have assumed in Eq. (19) that the AFM ordering occurs in a specific spin direction (along z axis).
We also generalize these MFT susceptibilities to finite temperatures in Appendix B. At finite temperatures, there are more available scattering processes as the lower effective band (band 1) is no longer fully filled, leading to several additional terms in the expressions for the susceptibilities.

D. RPA susceptibilities
The MFT results of the previous subsection capture only the ground state order and single particle excitations. The latter have energies of order U at U/J ≫ 1, and hence only describe high energy excitations. How-ever, the low energy physics is that for spin, coming from the effective Heisenberg exchange interaction at large U/J. In particular, there are low energy quantum fluctuations that lead to gapless collective excitations, corresponding to the spin waves of the Heisenberg antiferromagnet. The simplest approximate theory that can capture these collective modes is the RPA [36], which we will therefore employ here.
The basic idea behind the RPA is that certain terms in the perturbation expansion in U in the Dyson equation can be summed to infinite order as a geometric series (see Appendix A 2). Formally, at least for the nonmagnetically ordered state, RPA can be justified as a series expansion in the small parameter k F a s where k F is the Fermi wavevector, and a s is the scattering length related to U [60]. However, it is interesting that the RPA for the AFM ordered state also captures the collective modes at large U/J. These collective modes show up as bosonic gapless modes in the transverse spin susceptibility χ +− RPA . In Appendix A 2, we outline the derivation of the RPA susceptibilities, here we only state the relevant results.
For the density correlation, and for the longitudinal spin correlation, Note the only difference between these two is the sign of the term proportional to U in the denominator. Both have the form recognizable from summing a geometric series; indeed they have the same form as the more familiar RPA result for the non-ordered interacting Fermi gas.
The transverse spin susceptibility is more complicated: the AFM state doubles the unit cell and halves the Brillouin zone, coupling a spin g particle at k to a spin g hole at k + Q [see Eq. (26)]. This gives rise to the nondiagonal matrix form of Eq. (42) for the MFT susceptibility, which then feeds into the RPA one. The RBZ structure can be accommodated in a 2×2 matrix version of Dyson's equation to give (see Appendix A 2) Or in the explicit form given by [51], Here in Eqs. (48) and (49), as in Eqs. (40), (43) and (44), the momentum q belongs to the full BZ on both sides of the equation and χ +− Fig. 4. Figure 4 Figure 4(b) shows the real and imaginary parts of the denominator of Eq. (48). When the real part of the denominator crosses zero, a new pole emerges. We can estimate the energy of this pole as follows. The Heisenberg model has the spin wave dispersion given by [49] with γ q = (cos q x a + cos q y a)/2. The effective coupling constant J H is related to the Hubbard model parameters via J H = 4J 2 /U in the large U/J limit. Then, in Fig. 4 with q = (π/a, 0) and U = 5J, the pole occurs at energy ω q ∼ 8J 2 /U ∼ 1.6J. It can be seen that this is an overestimate, because we are not strictly in the U ≫ J limit. The imaginary part of both MFT and RPA transverse spin susceptibilities are compared in Fig. 5(a) to illustrate how the RPA renormalizes the single particle excitations and induces collective modes. The RPA transverse susceptibility notably differs from the MFT result owing to the collective excitations. In Fig. 5(b), we show the corresponding RPA renormalization of the MFT density and longitudinal spin susceptibilities. The different signs in the denominators in . Each term contributes with a Lorentzian shaped peak at ω = E k − E k+q . (b) Imaginary (solid) and real (dashed) parts of the denominator in the RPA case χ +− RPA (q, q; ω) [Eq. (48)]. The inset in (b) zooms in close to the origin and shows the zero-crossing in the real part.
Eqs. (45) vs. (46) indicate that the RPA corrections have opposite effects on the MFT density and longitudinal spin susceptibilities, as shown in Fig. 5(b).

RPA correction to the AFM order parameter
One can also incorporate quantum fluctuations and the effect of collective excitations in the calculation of the AFM order parameter m by including the RPA corrections to the MFT result of solving the MFT Eq. (33). Indeed, quantum fluctuations are known to strongly sup-  (40) and (41). However, the RPA renormalizes these two susceptibilities differently.
press m (by about 40%) [49] in the Heisenberg model. For the half-filled Hubbard model, Schrieffer et al. [36] have computed the RPA corrections to m numerically (see their Fig. 8). We will use their RPA-corrected data for computing the elastic scattering intensity in Secs. VI B 1 and VII B.

A. Scattered intensity
In the previous Section we showed how the timeordered and normally-ordered correlation functions can be calculated for the AFM ordered fermionic atoms in an optical lattice. Here we show how these quantities are related to measurements on the light scattered from the atoms. In optical lattice systems the atoms can be confined by highly anisotropic trapping potentials in which the atom dynamics is restricted to 1D or 2D. This makes the atomic samples particularly suitable for imaging. For the incident light tuned off from the atomic transition resonance frequency, the sample is then optically thin and the quantum statistical correlations of the atoms can be mapped onto fluctuations of the scattered light [20][21][22][23][24][25][26][27][28][29][30][31][32][33]. Measurements on the scattered light therefore convey information about the many-body state of the atoms and can be employed in the diagnostics of the correlated phases of the ultracold atoms in the lattice. Moreover, it was shown in Ref. [24] that by collecting the scattered light into the forward direction by a lens with the diffraction maxima blocked, an experimentally feasible thermometer for fermionic atoms can be realized. The sensitivity of the thermometer was analyzed by comparing the shot noise of the scattered light with fluctuations in the far-field diffraction pattern that arise from thermal correlations of the atoms. In this paper we study quantum z y x FIG. 6. Schematic illustration of the light scattering setup. The atoms are confined in the optical lattice close to the xy plane. The incident light field with the wavevector k1 propagates perpendicular to lattice in the positive z direction. The wavevector of the scattered light is denoted by k2 with the scattering direction determined by the coordinates θ and φ correlations in AFM-ordered strongly-correlated phase of fermionic atoms in an optical lattice. In this Section we introduce the formalism describing the relationship between optical signal (the intensity and spectrum of the scattered light) and the atomic correlations.
We assume that the atoms in the lattice are illuminated by light that can be approximated by a monochromatic plane wave with the frequency Ω in , propagating perpendicular to the lattice (in the positive z direction). The setup is illustrated in Fig. 6. We write the positive frequency electric field component E + in as whereê in and k 1 = kê z (k = Ω in /c) denote the polarization and wavevector, respectively, of the incoming light. In an optically thin sample, the dynamics of the electronically excited atomic state may be adiabatically eliminated and the scattered field amplitude E + sc (r, t) is proportional to the transition amplitude of atoms between the initial and final hyperfine electronic ground states g and g ′ [61], (52) Here the scattered field at r is evaluated in the far radiation zone, so that |r − r ′ | ≃ r −n · r ′ witĥ and the origin is located inside the atomic sample. The integration is over all the radiating atomic dipole sources at the positions r ′ . The field is scattered in the direction k 2 and the change of the wavevector of light upon scattering is, see Fig. 6, In Eq. (52) the effects of the level structure are incorporated in Λ g ′ g defined by Here the atomic transition dipole matrix elements between the ground state g and the excited state e are where e|σg denotes the corresponding Clebsch-Gordan coefficient. The summation runs over the circular polarization vectorsê σ (σ = −1, 0, 1), and D is the reduced dipole matrix element. The latter is related to the Weisskopf-Wigner radiative resonance linewidth γ by where ǫ 0 is the vacuum permittivity. In Eq. (52) we defined the prefactor C by Here δ denotes the detuning of the incident light frequency Ω in from the atomic resonance frequency ω 0 . Typically the wavenumber for the probe light is not the same as the effective wavenumber of the optical lattice potential. In order to suppress the spontaneous emission due to the lattice light lasers, the lattice lasers are considerably more off resonant. The probe and the lattice lasers may also be tuned to different transitions and the lattice spacing can be modulated by optical components, resulting, e.g., in the accordion lattices where the lattice spacing may vary for a given lattice light laser frequency [43][44][45]. We investigate the effect of different ratios κ between the probe light wavenumber k and the effective lattice light wavenumber k l [Eq. (11)], so that The intensity of the scattered light is given by where c denotes the speed of light in vacuum. By substituting the field amplitudes from Eq. (52) in the expression for the light intensity, we obtain the intensity in terms of the atomic correlation functions. The atoms are initially assumed to occupy the ground state in the lowest energy band. In the tight-binding regime the atomic field operators in the correlation function may be expanded in the series of the Wannier functions [Eq. (12)] whereĉ jg denotes the annihilation operator for the atoms in the electronic ground state g and the lattice site j = (j x , j y ). We obtain for the scattered light intensity Here we have defined where I in denotes the intensity of the incoming light, and In deriving Eq. (62) we have assumed that the spatial overlap between different sites is negligible. If the lattice potential is approximately independent of the different ground state levels g, the Debye-Waller factor α ∆k exhibits a simple format of the Fourier transform of the lattice site density that can be evaluated by means of the Wannier functions [Eq. (15)], where the oscillator length l i is defined by Eq. (16). In the expression (62) for the scattered light intensity, the spatial variation of atomic correlations is encapsulated in the operatorsĉ ig . The result is general and also includes the cases where the translational invariance of the lattice is broken owing to finite-size effects. In this work, we neglect any additional potential superposed with the lattice that would lead to a nonuniform density distribution, so that ĉ † jgĉ jg is here constant. The spatial profile in Eq. (62) is therefore solely determined by the phase factors ofĉ ig 's. The simple relationship (62) between the scattered light intensity and the atomic correlation functions is a consequence of the weak off-resonant coupling of light. For near-resonant light the coupling is strong even in a 2D lattice and results in excitations of collective polarization modes [62].
We will consider the scattering processes of the atoms to higher energy bands in Sec. V C. Although such processes result in the photon frequencies that are shifted by the energy difference between the bands and can be filtered from the signal, they can contribute to the heating rate of the atoms in the lattice.
In the specific analysis of the optical signatures of the atomic correlations it is beneficial to separate in the scattered intensity the contributions from the elastic and inelastic scattering events. In the following study we define the elastic scattering processes as those in which the atom scatters back to its original momentum state. We evaluate Eq. (62) in terms of the elastically and inelastically scattered light intensities I e (∆k) and I i (∆k), respectively, Here∆k denotes the change of wave vector of light on the xy plane and the density operatorρ kg for the spin state g is given by Eq. (6). We have defined and the static structure factor All the q = q ′ = 0 terms in the previous equation are included in the elastic part. This corresponds to incorporating only the connected Feynman diagrams in the correlation function of the static structure factor (indicated by the subscript c) and the disconnected ones in the relevant expansion are precisely those that go into the elastic part, see Ch. 13.4 of [63]. The finite-size effects of the lattice contribute in Eq. (70) via the u k factors. In the limit of a large lattice, we may approximate S g3g4 g2g1 by translationally invariant atomic mode functions, so that the summation over the sites approaches a delta function u k → N 2 s δ k,0 and In a typical 40×40 lattice we are investigating, taking the continuum limit changes the integrated inelastically scattered light intensity by less than 2%. In order to calculate the intensity of the scattered light [Eq. (68)] and the corresponding static structure factor [Eq. (70)] it is useful to define a static response function as In Sec. IV C 2 and IV D we showed how both the RPA correlations based on the Feynman-Dyson perturbation series as well as the MFT correlations for the AFM state in the large lattice limit can then be efficiently calculated in a compact form by evaluating the diagonal components S g3g4 g2g1 (q, q). In this paper we approximate the inelastically scattered light intensity by these diagonal expressions, while still including finite-size contribution from the diffraction pattern via |u∆ k−q | 2 . The corresponding intensity expression reads In Sec. VI B we will show how S g3g4 g2g1 (q, q) can be related to the density as well as longitudinal and transverse spin correlation functions [Eq. (95)].
The general Eqs. (66)- (68) give the scattered light intensity for an arbitrary lattice system. The scattered light carries information about the atomic correlation functions. The lattice structure generates the diffraction pattern and the overall envelope of the pattern is produced by the Debye-Waller factors that depend on the profile of the atomic wave functions on individual sites. The dependence of the scattered light on the polarization, atomic level structure, and the scattering direction is incorporated in M g3g4 g2g1 . The elastic scattering produces a diffraction pattern from a non-fluctuating atom density in the lattice where the Wannier site wave functions play the role of the diffraction slit profile [24]. The inelastic scattering processes are those in which an atom scatters from one quasimomentum state to another different state. The inelastic scattering is sensitive to the fluctuations of the atoms and reflects the underlying statistical correlations between the atoms. It produces scattered light into angles outside the diffraction orders, generating fluctuating shifts in the diffraction pattern that result from the atomlattice system absorbing recoil kicks from the scattered photons [24].
The role of the elastically scattered light intensity is easiest to analyze in the case of a uniformly filled lattice (when the translation symmetry of the lattice is not broken by any of the atomic species). In that case the density operator expectation value reads where f g is the atomic filling factor of species g (the total number of atoms of species g divided by the total number of sites, N 2 s ).∆k is defined after Eq. (68). We then find that the elastically scattered light intensity, is determined by the Bragg diffraction pattern of the lattice |u∆ k | 2 , weighted by the contributions from the atomic level structure via Λ gg , and modulated by the Debye-Waller factor α ∆k . In the specific case of a 2D square lattice we obtain the familiar diffraction pattern of a 2D square array of N s × N s diffracting apertures For uniformly filled lattice the elastic part contains no information about the atomic correlations in the system. Since the atom statistics is mapped onto the inelastically scattered light intensity according to Eq. (68), it is beneficial to block the elastically scattered light before the measurement [24]. We will explain this procedure in detail in Sec. VI B 1. For the AFM state studied in Sec. IV B, lattice translation symmetry is broken, resulting in the Brillouin Zone being halved in size. We shall show in Sec. VI B 1 that in this case, there are new diffraction peaks that correspond to this new lattice periodicity, and is therefore a key signature of the AFM order. This effect was first analyzed in Ref. [34].

B. Scattered spectrum
In the previous Section we established the relation between the scattered light intensity with the equal time atomic correlations [Eqs. (66)- (68)]. We now analyze the spectrum of the scattered light and show how it conveys information about the excitation spectrum of the atoms in the lattice. The scattered light spectrum may be obtained as a Fourier transform of the two-time correlation function of the scattered electric field [61] S(∆k, ω) = A dt e iωt E − (r, 0)E + (r, t) , where A denotes the normalization factor. We can then write the scattered spectrum in terms of the two-time correlation functions of the atoms in the optical lattice. The spectrum can be separated into an elastic and inelastic component and we obtain where A ′ ≡ AB/(2ǫ 0 c). In the last equation we have introduced the dynamic structure factor, which is analogous to the static case of Eq. (70) The elastic component corresponds to a peak at ω = 0. The subscript c indicates the connected diagrams for which ω = 0, revealing the excitations of the system [see Sec. VIII]. Analogously to Eq. (72) we define the dynamical response functions as We approximate Eq. (80) in a similar fashion as Eq. (73) C. Atom losses to higher bands In Sec. V A we calculated the optical intensity signal for probing the atoms in an optical lattice. This consists of the elastic scattering processes in which the final state of the atoms is the same as the initial state as well as the inelastic scattering processes within the lowest energy band where the quasimomentum state of the atoms changes. The atoms that initially occupy the lowest energy band may also undergo scattering to higher bands. Owing to the energy splitting between the adjacent bands, which is of the order of 2s 1/2 E R [E R was defined in Eq. (10)], the photons that scatter to higher bands are frequency-shifted from the optical signal and could be filtered out. This is because the maximum recoil kick absorbed by the atom within the lowest band on the xy plane is k [the recoil component on the lattice plane satisfies |∆q| = k sin(θ), see Eq. (54)], corresponding to the energy shift of κ 2 E R [Eqs. (10) and (60)], which is less than the energy difference between the bands [64]. The scattering to higher bands, however, provides a loss mechanism which we will estimate when calculating the measurement accuracy of AFM correlations of the atoms.
Here we extend the analysis of the loss rates of Ref. [28] to our multi-level formalism. In the evaluation of the scattered intensity the following correlation functions involving states in higher energy bands yield nonvanishing contributions whereĉ mjg denotes the annihilation operator for the atoms in the band m, site j, and ground state g, Eq. (12). The nonvanishing contribution from the empty excited energy band m = 0 results from ĉ mjg3ĉ † mjg2 by the creation of an atom at (m, j) followed by the annihilation of an atom at (m, j). We concentrate on the case that the interactions do not mix the spin states, so that only the term g 3 = g 4 is nonvanishing. The coefficients G are defined in terms of the Wannier function integrals This can be used to simplify Eqs. (84) and (85). We find that the contribution to the scattered intensity of this process reads where I hb refers to the scattered light intensity resulting from the scattering events where the atoms end up in the higher energy bands. Here N = N 2 s is the total number of atoms in the lattice. This result remains valid for both fermions and bosons as long as the assumption of unpopulated higher bands is valid.

VI. OPTICAL SIGNATURES OF MAGNETIC ORDERING IN SCATTERED INTENSITY
A. Two-species atomic gas of 40 K In the previous Section we presented the general expressions for the dependence of the scattered light on the atomic correlation functions in an optical lattice system. Next we analyze the 2D square lattice system of two-species fermionic atoms introduced in Sec. III, with equal population of N 2 s /2 atoms of both species. As a specific example we consider 40 K, which has been used in an experimental realization of fermionic Mott insulator states in lattices [6][7][8] and recently AFM ordering [8]. We consider the two electronic ground states The incident field is assumed to be σ − polarized, so that the two ground states are coupled to the electronically excited states |1 = |4P 3/2 , F e = 11/2, m F = −11/2 , |2 = |4P 3/2 , F e = 11/2, m F = −9/2 .
The level scheme and the transitions are illustrated in Fig. 7. The atoms in | ↓ undergo a cycling transition in which case they are only excited to |1 , decaying back to the original state | ↓ . The atoms in | ↑ are excited to |2 from where they can decay either back to | ↑ or to | ↓ . The latter represents a spin-exchanging transition. We will find that the transitions in which the spin is changed convey information about transverse spin correlation functions, while those associated with the scattering processes in which the spin is conserved are proportional to density and longitudinal spin correlation functions. Specifically, the different transitions can be identified with the different RPA susceptibilities in Eqs. (45), (46) and (47), respectively. In our system we vary the lattice size between N s = 16 and 40 sites along each direction. We consider two lattice heights, 7.8E R and 25E R . The trap frequency perpendicular to the lattice is chosen as ω z = 10E R / . For 40 K we take experimentally realistic [6][7][8] values for the incoming light λ = 766.5nm and I in = 5W/m 2 , and assume it to be detuned from the atomic resonance δ = 20γ [Eq. (59)]. This yields in Eq. (63) Br 2 ≈ 1615 photons/s. We vary the ratio between the lattice spacing and the wavelength of the incident light by changing κ [Eq. (60)] and take κ = 0.66, 1.05 and 1.5. All these correspond to subwavelength lattice spacing, but the additional magnetic peak due to period doubling may only be observed for κ = 1.5 > √ 2.

B. Scattered intensity
In Sec. V A the total scattered intensity was separated into the elastic component I e (∆k) and an intraband inelastic component I i (∆k), see Eqs. (67) and (68). Furthermore, in Sec. V C we have taken into account the inelastic scattering of atoms to higher bands in Eq. (87), I hb (∆k), that is not assumed to contribute to the detected signal but affects the heating rate of the atoms. The total scattered intensity is thus We now analyze each of these components for the specific case of a two-species 40 K. Applying the level structure of Fig. 7 we find the explicit expressions for all the scattered elastic and inter-band intensity components Here the Debye-Waller factor α ∆k is given in Eq.
We now consider the scattered inelastic intraband intensity component I i . Due to the broken translation symmetry in the AFM state, the RPA susceptibilities of Eqs. (45), (46) and (47) .
(94) The static response functions that appear in each element of the previous matrix are defined in of Eq. (72). To relate to the RPA density and spin susceptibilities [Eqs. (45), (46), and (47)] we first note that using Eqs. (6), (5) and (4), we can define the matrix of static response functions for the density, longitudinal and transverse spin operators where η(g) is defined in Eq. (22). These static response functions can in turn be related to the time-ordered correlation functions (susceptibilities) of Eqs. (45), (46) and (47) of Sec. IV, using the relationship in Eqs. (36) and (35) together with frequency integration from Eq. (37). The diffraction factors u q defined in Eq. (69) can also be accommodated into the RBZ structure by defining Hence, the inelastic component of the intensity of light scattering off atoms in the AFM state of Eq. (73) can be written as In the following we will analyze the different scattering contributions. The spectrum and the excitations will be studied in Sec. VIII.

Elastic scattering
This Section is devoted to studying the elastic component of the scattered light. We show that the emergence of AFM ordering in the system is directly observable in the elastically scattered light intensity as this results in magnetic Bragg peaks in the scattered light signal. For the two-species system studied here (Fig. 7) where f g = 1/2 is the atomic filling factor of species g at half-filling and η(g) is defined in Eq. (22). Comparing to Eq. (74), the new term is proportional to the AFM order parameter m, Eq. (19), and is centered at the ordering vector Q = (π/a, π/a). The order parameter m can be obtained by solving the implicit MFT [Eq. (33)]. The MFT phase diagram of Fig. 2 shows its dependence on T and U . However, as discussed in Sec. IV D, quantum fluctuations around the MFT solution have significant effects not only on susceptibilities, but also on the AFM order parameter m: the ordering can decrease by up to ∼ 40% at large U . Thus, we shall use the RPA-corrected m values as computed by Schrieffer et al. [36]. Substituting Eq. (100) into Eq. (98), we see that in addition to the usual diffraction term u * ∆ k u∆ k centered at ∆k = 0, there is a new magnetic Bragg peak centered at∆k = Q, proportional to m 2 (see Fig. 8). This new peak can be detected by collecting scattered light around ∆k ∼ Q, which, according to the definition of ∆k in Eq. (54), corresponds to θ B = arcsin √ 2 κ and φ B = π/4 (see Sec. VII A). The position of the magnetic peak depends on the ratio κ [Eq. (60)] between the probe light wavevector and the effective wavevector of the lattice light, see Eqs. (54). Hence it will only be observable if κ ≥ √ 2. Magnetic Bragg peaks were first studied in optical lattices in Ref. [34] and experimentally observed for an artificially prepared density pattern in Ref. [35].
The dependence of the magnetic Bragg peak on the staggered magnetization m is illustrated in Fig. 8(a) that shows the elastically scattered intensity from a lattice  (102) The resulting magnetic Bragg peak is now strongly enhanced, see Fig. 8  tween the scattered light from the two species. In fact, the magnetic Bragg peak is observable only because the dipole transition matrix elements between the two species are not equal and, according to Eqs. (98) and (100), is weaker than the single species response by the factor of (2/11) 2 . If the transition matrix elements were the same, the spin-independent imaging would only probe the total density without revealing the AFM order.

Elastic scattering in the presence of short-range correlations
So far we have considered long-range AFM order where the staggered magnetization m is constant throughout the lattice. In Sec. IV A, we discussed how there is no genuine long range AFM order in the Hubbard model (or its strong coupling limit, the Heisenberg model), at any finite temperature in the thermodynamic limit. It is beyond the scope of this paper to calculate the finite temperature short-range effects, but we can simulate shortrange ordering effects in a phenomenological manner, as follows: At finite temperatures, the spins are correlated up to the AFM correlation length ξ AFM which leads to domains of size ∼ ξ AFM . We introduce the spatial variation of the AFM order parameter by letting m i to depend on the site index i, with the amplitude fixed at the T = 0 magnetization value m RPA . For simplicity, we assume that m i 's in different domains are not oriented in random directions, but that there exists a preferred axis, generated, e.g., by a small imbalance in the Fermi levels of the two species. We therefore introduce a sign s i = ±1 (Ising variable) for m i 's that fluctuates from site to site The configuration of s i is modelled using the nearestneighbor Ising Hamiltonian J Ising is the coupling strength, and K = J Ising /k B T . For J Ising < 0, the ground state of the system at T = 0 is an AFM Neél state. In contrast to the Heisenberg model, the Ising model has a nonzero critical temperature. The transition temperature for a finite system with periodic boundary conditions can be estimated as [65] K c (N s = 40) ≈ −0.437 .
We approximate the the correlation length ξ Ising (K) for a finite-size system (N s = 40) and for T > T c (40) [K c (40) < K] by [66] ξ Using the Wolff algorithm [67], we can numerically generate a specific configuration for the Ising variable s i , with a correlation length determined by the temperature in the Ising model. To ensure that the number of up and down spins is equal, we impose a constrain |m F | = | i s i | < 0.01 on the ferromagnetic Ising order parameter. Figure 9 shows two of the generated configurations at two different temperatures. A specific configuration of m i is then used in Eq. (99) to calculate the single-species (↓ atom only) elastic part of the scattered intensity of Eq. (101) as Note that this expression does not contain the diffraction peak in the forward direction, as we are only interested in the magnetic Bragg peaks whose contribution is significant around the perpendicular direction. The numerically calculated scattered light intensity at different temperatures is shown in Fig. 10(a)-(c). The average value of the order parameter m is smaller than 0.02 in all the cases. The average scattered light intensity grows at temperatures close to the transition [ Fig. 10(d)], and the magnetic Bragg peaks emerge as the correlation length increases. The large fluctuations between different realizations in the simplified model (Fig. 10) are likely to significantly overestimate the fluctuations of the true AFM state.

Inelastic intraband scattering
As was demonstrated in the previous Section, the elastic part of the scattered light intensity conveys information about the atom density in the lattice and generates the diffraction pattern of the atomic lattice structure. If the detected signal cannot distinguish the two spin components, the light provides almost no information about the AFM order. On the other hand, when the contributions of the two spin components can be separated in the scattered light, the emerging AFM order and the period doubling can be identified as additional Bragg peaks. In addition to the elastic signal, one may also study the inelastically scattered light intensity [Eq. (68)]. The inelastic scattering processes are proportional to static structure factors S g3g4 g2g1 (∆k) [Eq. (70)] that represent scattering events in which an atom is excited from a quasimomentum state q and scatters to a different quasimomentum state q ′ . Atoms absorb recoil kicks from the scattered photons. The recoil events depend on the statistical correlations between the atoms, generating fluctuating shifts in the diffraction pattern and significant scattering outside the diffraction orders. In the process, the atomic correlations are mapped onto the properties of the emitted light. Inelastically scattered light in a singlecomponent fermionic gas in a lattice can reveal thermal correlations [24] and has in a two-component case previously been proposed as a detection method for topological order of the atoms [20].
The inelastically scattered light intensity for the twocomponent 40 K gas is given by Eq. (97). The scattering contributions in which the spin is conserved are proportional to density and longitudinal spin susceptibilities. The spin-exchanging transitions (see Fig. 7) generate the term depending on the transverse spin susceptibility. The two processes exhibit very different angular distribution of the scattered light as shown in Eq. (93). In the spin- conserving processes the emitted photons are generated by the σ + transition in which case the intensity in the forward direction is twice the intensity in the perpendicular direction. The spin-exchanging process, on the other hand, produces scattered photons via the π transition which is oriented parallel to the propagation direction of the incident field. Therefore the scattering reaches its maximum in the perpendicular direction (θ = π/2) and entirely vanishes in the forward direction. We have calculated the angular distribution of the inelastically scattered light for different values of the on-site interaction strength U , and hence the staggered magnetization m of the AFM ordering. In Fig. 11 we show the results based on MFT at T = 0. In this case the intensity is obtained using the MFT static structure factor [Eqs. (B1) and (B2)] in Eq. (68). The corresponding MFT susceptibilities required in the calculation are provided by Eqs. (40)- (44), (35), (36), and (37). The different angular distribution of the different scattering contributions is clearly visible in Fig. 11. We find that in MFT the density and longitudinal spin susceptibilities are more sensitive to the variation of U than the transverse susceptibility.
The intensity calculations based on MFT fail to cap- The difference between the two approaches in the scattered intensity distribution is illustrated in Fig. 12. The low-energy collective modes are notable in the transverse spin correlations, corresponding to the spin-flip transitions, but in the case of the spin-conserving scattering processes the two approaches yield almost identical intensity distributions. The scattered light intensity distributions based on the calculation of the atomic correlations within RPA at T = 0 for different U are shown in Fig. 13. The variation of the signal as a function of U and the AFM order are most notable in the perpendicular direction, compared with the MFT case of Fig. 11. Finally, to produce an example how the scattered signal also depends on the temperature of the system, we show the calculated light intensity distributions for different T based on the MFT Eq. Although the collective excitations in this example are ignored, the comparison between the MFT and RPA T = 0 results suggests that the scattered intensity is less sensitive to low-energy collective excitations in the near-forward direction where the the spin-conserving transitions are dominant. A more accurate description of the temperature dependence would require a finitetemperature version of RPA which is beyond the scope of the present study.
We find in Fig. 14(a) a significant dependence of the light intensity on the temperature T of the atoms in the near-forward direction that is analogous to the temperature sensitivity of a single-component noninteracting fermionic gas [24]. (Note that in this figure, lower T is represented by increased staggered magnetization m, for a given fixed value of U = 5.3J.) The suppression of small-angle scattering at low T can be understood in terms of the Fermi blocking: the scattering events in which an atom would recoil to an already occupied state are forbidden and in MFT the small-momentum recoil events can take atoms out of the Fermi sea only near the Fermi surface. Owing to the sensitivity of the signal to temperature, optical diagnostics could be used as a thermometer of the atoms also for an interacting twocomponent case. We study the detection accuracy of this method in Sec. VII B.

Inelastic losses to higher bands
The detected light intensity consists of the elastic and inelastic intraband components that were calculated for 40 K in Secs. VI B 1 and VI B 3. In addition, the atoms can scatter to higher bands as demonstrated in Sec. V C. The interband scattering can be separated from the detected signal (owing to the different frequency of the photons), but still contributes to the heating rate of the atoms. For the level structure of 40 K we can write the intensity of We show the angular distribution in Fig. 15. The inelastic interband scattering is proportional to 1 − α ∆k , where α ∆k denotes the Debye-Waller factor [Eq. (65)]. On the other hand, the intraband inelastic scattering is proportional to α ∆k . Consequently, in deep lattices atoms are more strongly confined and the loss rate to higher bands can be suppressed, while the inelastic intraband scattering, that provides information on correlations, can be enhanced.

A. Optical components for light detection
As discussed in Sec. VI, the elastically and inelastically scattered light from the atoms can provide different information on the AFM order parameter m, the temperature of the system, and collective and single particle excitations. Here in this subsection, we discuss in sequence how one might optimize the experimental configurations for i) detection of the extra intensity peak due to AFM ordering and the measurement of the AFM order parameter m, ii) temperature measurement, and iii) measuring the effects of atomic correlations on inelastically scattered light.
The elastically scattered light intensity in the presence of the AFM ordering contains extra peaks that result from the period doubling of the atom density when atoms in only one of the two spin states are measured [Eqs. (98) and (100)]. We consider detection of this emerging Q peak by a small optical lens when the lens is placed at the appropriate angle so as to maximize the contribution of the peak [ Fig. 16(a)]. This then can be used to to measure directly the AFM order parameter m.
In order to detect atomic correlations or to measure temperature of the system from inelastically scattered light, we consider two experimental configurations: i) a lens is placed in the near-forward direction, as depicted in Fig. 16(b), and ii) a lens is centered close to the perpendicular direction, see Fig. 16(c). The total light intensity collected by the lens (L) with a given NA can be obtained by integrating the scattered intensity, Eqs. (91)-(97), over the solid angle (dΩ = sin θdθdφ) determined by the corresponding scattering angles, where ∆k denotes the change of momentum upon scattering [Eq. (54)] and the index α = e, i, hb refers to the elastic, inelastic intraband or inelastic interband (higher band) components, respectively Eqs. In (a) the elastically scattered light corresponding to the emerging additional Bragg peak generated by the AFM ordering is collected by a small lens. In (b) the setup is closely related to that of Ref. [24]. The two lenses have focal lengths f1 and f2. The light scattered in the near-forward direction is first collected by lens 1. In the focal plane the scattered light is selectively stopped by a block in order to suppress the intensity of the elastically scattered light at the detector. The shape and the size of the block can be optimized for different measurements of collective excitations or of temperature. In (c) the scattered light is collected near the perpendicular scattering direction of θ = π/2, to measure transverse spin correlations.
in temperature leads to enhanced inelastic scattering in the near-forward direction, see Fig. 14. However, a lens placed near the forward direction will also capture the much stronger elastic scattering signal, see, e.g., Fig. 15. Now, the elastically scattered light generates a diffraction pattern [Eq. (76)] and, for the subwavelength lattice spacing we consider, only the zeroth order Bragg peak is collected by the lens. Thus, in order to maximize the proportion of inelastically scattered light in the mea- surement, we can block the high-intensity regions of the elastically scattered light by placing an appropriately designed block on the focal plane of the lens, see Fig. 16(b). This setup is similar to the one proposed in Ref. [24] to measure the temperature of an ideal single-species fermionic gas in a lattice. In that case the inelastically scattered light varied as a function of the temperature in the near-forward direction resulting from the enhancement of the Pauli blocking effect at low temperatures.
In the two-species interacting case within the MFT we observe the same behavior as shown in Fig. 14. We correspondingly achieve the optimized detection efficiency by selecting a narrow cross-shaped block that covers the Bragg peak and high-intensity regions along the principal axis of the lattice, analogously to Ref. [24], as shown in Figs. 19(a) and (b).
We now consider how to optimize the detection of AFM correlations in the inelastically scattered light. At zero temperature, the effect of the AFM ordering on the angular distribution of the scattered light calculated with RPA is displayed in Fig. 13. The staggered magnetization most noticeably changes the signal away from the forward direction. Optimal shape and size of the block can be estimated from the angular distribution of scattered light that shows both the elastic and inelastic contributions (see Fig. 15). Since near-forward scattering provides only little information on magnetization in this case, it is beneficial to consider a block that combines a narrow cross block with a circular block located at the center. The light can then be collected with a lens of large NA. For a lattice of 40×40 sites the block sizes used are listed in Table I. In Fig. 16(b) we show an example circular block with the angular size Θ Circ-Block ≈ 0.20 rad.
As has been discussed in Sec. VI, spin-exchanging scattering processes dominate the inelastic signal near the perpendicular direction [ Fig. 16(c)]. For κ < √ 2 the elastic scattering contribution to the collected signal is negligible (Fig. 15) and no block is needed. On the other hand, for κ ≥ √ 2 the elastic component near the perpendicular direction strongly depends on the AFM order parameter m [see Fig. 8, and Sec. VI B 1, in particular Eqs. (98) and (102)]. As we will show in next subsection, this enhances the sensitivity of the signal to changes in m.

B. Measurement accuracy
We will analyze the accuracy of the optical measurements of the AFM correlations in the lattice when the light is collected by a lens. We will follow the procedure introduced in Ref. [24] where the optical detection accuracy of temperature in a single-species fermionic atomic gas in an optical lattice was calculated. In an inelastic scattering event an atom scatters to a different quasimomentum state owing to the photon recoil kick. Inelastic scattering leads to heating of the atomic gas and perturbs the many-body state of the atoms. In order to limit the effect of heating, the number of inelastic scattering events in each experimental realization of the lattice system should be limited to a small fraction of the total number of atoms in the lattice. We set the maximum number of allowed inelastic scattering events to be W , so that W/N 2 s is sufficiently small. In the example analysis we take W/N 2 s = 0.1. We specify the fraction of inelastically scattered photons that are collected by the lens and constitute the measured signal for a given magnetization m by η(m), Here I tot i (m) and I tot hb denote the total rate of inelastic intraband and interband scattering events, respectively. We assume that the scattered light corresponding to interband transitions is filtered out of the signal so the corresponding rate is excluded from the numerator. If, for simplicity, we assume a 100% photon detector efficiency, the number of detected inelastically scattered photons in each experimental realization of the lattice system is given by N i c (m) = η(m)W . If the lattice system is prepared and the experiment is repeated τ times, we find that the total number of detected photons is N e c (m) = where N e c (m) denotes the total number of detected elastically scattered photons in a single experimental realization of the lattice system and I L e (m) is the scattering rate of elastically scattered photons that are collected by the lens. In optical diagnostics of AFM ordering one would need to distinguish in the scattered light signal two fermionic states in the lattice that exhibit different magnetic orderings m 1 and m 2 . After τ experimental realizations the difference in the total number of detected photons between the two ordered states is τ [N c (m 2 ) − N c (m 1 )]. The minimum requirement for the two states to be distinguishable is that this difference is at least equal to the photon shot-noise τ N c (m 2 ), so that τ [N c (m 2 ) − N c (m 1 )] τ N c (m 2 ). Thus the minimum number of experimental realizations τ min required to distinguish between the optical signal from two magnetization states m 1 and m 2 approximately satisfies In the rest of this Section, we present results for the rate of detected photons as a function of the staggered magnetization m and the number of experimental realizations τ required to determine changes in the AFM order parameter m. The changes are calculated with respect to a reference value, m ref , and for a given relative accuracy ∆m/m ref [Eq. (113)]. At T = 0, the inelastic scattering is calculated using the RPA susceptibilities. In  (19). Except the scaling analysis of τ with lattice size, all the results in this Section are for a lattice of size 40×40.
We first study the AFM order parameter measurement accuracy when a lens (NA= 0.2) is used to collect the light along the direction of the emerging magnetic Bragg peak [ Fig. 16(a)]. We also study the configuration in which a lens (NA= 0.4) is placed perpendicular to the propagation direction of the incident probe laser [ Fig. 16(c)]. In both cases, the ratio between the wavenumber of probe light and the effective lattice light κ = 1.5 > √ 2 [Eq. (60)]. This allows the NA= 0.4 lens to collect the signal also from two magnetic Bragg peaks [Eq. (100)]. In Sec. VI B 1, we have shown that the elastic signal is strongly enhanced when only the ↓ atoms are detected. Figure  In both lens configurations the number of experimental realizations τ needed to achieve a given accuracy drops dramatically when only ↓ atoms are detected. Some specific values of τ for single species detection are presented in Table II. These values are (much) lower than all other experimental configurations to be presented later [see Table III]. We conclude that when available, the single-species detection scheme provides the most accurate determination of the AFM order. Experimentally, the single-species imaging can be realized by transferring the other spin component to a different hyperfine state [35].
In Fig. 18 we show the effect of lattice depth s [Eq. (9)] and κ on the detection accuracy of m. The scattered light is collected by a lens of NA=0.8 in the forward direction, corresponding to the experimental arrangement of Fig. 16(b). We find that a deeper lattice generally enhances the scattering rate, but the effect of s on τ in the studied cases is negligible [ Fig. 18(a),(b)]. Some example values are shown in Table III. The number of required experimental realizations for a 40×40 lattice drops rapidly  when the desired accuracy is reduced and the staggered magnetization is increased.
Similarly, increasing κ enhances the rate of detected photons [Figs. 18(c),(d)]. Generally however, the detection accuracy is lower for larger values of κ. This is so because for larger κ, there are more inelastic scattering events particularly near the ordering wavevector. Such inelastic signal contributes to inelastic losses and heating, but is not captured by the forward direction lens   Fig. 18(b). The parameters of the block are given in Table I B 3), we only provide a qualitative analysis of finite temperature effects using MFT, without taking into account the collective excitations included in RPA. Using a narrow cross-shaped block of width 0.048 rad, the light can be collected near the forward direction where the temperature strongly affects the scattering rate (Fig. 14). Lower T corresponds to stronger magnetization values and fewer collected photons, as shown in Fig. 19(a), and the temperature changes in m can be accurately detected [ Fig. 19 Finally in Fig. 20, we show the RPA results for the accuracy in m RPA when the scattered light is collected perpendicular to the propagation direction of the incident laser [ Fig. 16(c)]. The light scattered from the spinconserving and spin-exchanging transitions may be separated owing to the different frequency of the scattered photons (or the polarization, see the next Section). If the transitions are not separated, the measurement accuracy in the perpendicular direction is significantly lower than, e.g., for forward direction measurements. There is a notable improvement in the detection accuracy in the perpendicular direction when only the spin-exchanging scattering processes (representing the transverse spin correlations) are selected (for the angular distribution of the scattered light for the different components, see Fig. 13). For instance, for κ = 1.5, 10% accuracy for m ref ≃ 0.19 can now be achieved after 50 realizations for NA= 0.4. By increasing the size of the lens to NA= 0.5 this can be further improved to 40 realizations. In general, for the perpendicular direction the large κ = 1.5 case gives the highest number of scattered photons because of the strong enhancement of spin-exchanging scattering processes [cf. Fig. 21(a) vs Fig. 13(c)]. Also the τ values are substantially lower for the κ = 1.5 case [ Fig. 20(b),(d)].
Our example calculations are for a 40×40 lattice. Smaller values of τ can be obtained for larger lattices. The number of required experimental realizations of the lattice system τ is approximately inversely proportional to the number of sites τ ∝ N −2 s . We have simulated lattices sizes between N s = 16 and N s = 40, and our results confirm this scaling to be qualitatively accurate for reasonably large lattice systems N s 25 for both forward and perpendicular directions measurements. For smaller systems the choice of the block size and shape can result in larger variations owing to the dependence of the width of the diffraction peak on the lattice size.

C. Distinguishability of transitions by light polarization
In general the scattered light corresponding to different transitions are separated in frequency space and could be identified by filtering the relevant frequencies. If the sufficient frequency resolution is not achievable the transitions could still be partially distinguished by the polarization of the scattered light. The polarization of the scattered light from the spin-conserving and spin-exchanging transitions depend on the scattering direction and generally they are not orthogonal. For any given scattering direction we may project, e.g., the light scattered from the spin-exchanging transition to the direction of the polarization of the scattered light from the spin-conserving transitions. This provides the optimum value how much for the given scattering direction the light from the spinconserving transition can be distinguished from the total FIG. 21. Angular distribution of the inelastically scattered light intensity along the direction φ = π/4 when the scattered light is projected along specific polarization directions. We show different values of the interaction strength U at T = 0. The calculations are based on RPA, κ = 1.5 and the other parameters are as in Fig. 11. (a) The total scattered light intensity; (b) the total scattered intensity from the spin-conserving transition; (c) the total scattered intensity from the spin-exchanging transition; In (d) and (e) the light has been projected along the direction of the scattered light from the spin-conserving transition. (d) projected component of the scattered intensity from the spin-conserving transition which in this case equals to (b); (e) projected component of the scattered intensity from the spin-exchanging transition. In (f) and (g) the light has been projected along the direction of the scattered light from the spin-exchanging transition. (f) projected component of the scattered intensity from the spinconserving transition; (g) projected component of the scattered intensity from the spin-exchanging transition which in this case equals to (c). scattered intensity. In order to analyze this we modify the polarization vector of Eq. (55) Λ g ′ g by where the polarization of the scattered is given bŷ In a studied example case of Fig. 21 we show the angular dependence of the scattered light intensity and compare this with the projected intensities along the polarization of the scattered light from either of the two transitions. The relative part of the signal from the density and longitudinal spin correlations is enhanced in Figs. 21(d) and (e), while the contribution of the transverse spin correlations is particularly strong in Figs. 21(f) and (g) in the direction around θ ≃ π/2.

VIII. DIAGNOSTICS OF EXCITATIONS FROM SCATTERED SPECTRUM
In this Section we calculate the spectrum of scattered light and show how intraband inelastic scattering can reveal the single-particle and collective excitations in an AFM ordered lattice system. In Sec. V B, the scattered spectrum is split up into an elastic [Eq. (79)] and an inelastic part [Eq. (80)]. The elastically scattered light has no nontrivial spectral structure, consisting of only the zero frequency part, by definition. Inelastic losses due to scattering to higher bands occur at high frequency (on the order of the band gap) and can therefore be filtered out and ignored, as this part does not contain information about the state probed.
The inelastic spectrum of Eq. (80), or the more specific form of Eq. (83), contains useful information about the excitation spectrum of the system being probed. In Sec. IV D, we show how the MFT contains only single particle excitations in the density, longitudinal and transverse spin susceptibilities. The RPA partially takes into account quantum fluctuations around the AFM ordered state and can capture the collective excitations (spin waves) emerging in the transverse spin susceptibility. RPA also renormalizes the single particle excitations. (See Figs. 4 and 5 for comparisons of the RPA and MFT.) Via linear response theory (see Sec. IV C 1), the various RPA susceptibilities can be related to the dynamic structure factor Eq. (81), or the dynamic response function Eq. (82). The spectrum of inelastically scattered light reads This spectrum has been derived using the same procedures as the analogous expression for the intensity in Sec. VI B, generalizing the definitions of Eq. (95) to the frequency-dependent case here. We can translate the dynamical response function of Eq. (82) to the density, longitudinal and transverse spin dynamical response functions S S S ij (q, ω) (i, j = ρ, z, +, −) using a frequency dependent version of Eq. (95). In turn, these dynamical response functions are related via Eq. (35) to the various RPA spin and density susceptibilities Eqs. (45), (46) and (47) of Sec. IV, which are then used to compute the scattered spectrum. We first compare in Sec. VIII 1 the physical quantities of interest, the angle-resolved spectrum [Eq. (116)] and susceptibilities [Eqs. (45), (46) and (47)], with the angle-integrated spectra [Eq. (117)] that corresponds to a measurement of the spectrum by a lens over a range of scattering angles. We then analyze in Sec. VIII 2 the features corresponding to single particle excitations and in Sec. VIII 3, the collective mode peak.

Comparison between angle-resolved spectrum and angle-integrated spectrum
The scattering rate of off-resonantly illuminated atoms is generally low. Experimentally, the number of measured photons can be increased by collecting the inelastically scattered photons over a range of scattering angles using a lens (see Sec. VII A). The spectrum, however, changes with the scattering angle, since each ∆k represents a different argument of the dynamical structure factor. In this Section we analyze how much of the spectral structure can be extracted when the measured light is collected by a lens with a large NA and the spectrum is integrated over the corresponding range of scattering angles. We therefore define where L = F, P indicates the experimental configuration: F for the lens with a given NA pointing in the forward direction [ Fig. 16(b)] and P for the perpendicular direction [ Fig. 16(c)]. In Fig. 22, we compare the angular integrated spectrum [Eq. (117)] for two different lens sizes with the  Fig. 16(b)]. The angle-resolved quantities, the spectrum S(∆k, ω) and total susceptibility χ(∆k, ∆k; ω) are computed at a wavevector close to the axis of the lens (θ ∼ 0, φ ∼ 0). (b) shows the case for a lens in the perpendicular direction [ Fig. 16(c)]. S(∆k, ω) and χ(∆k, ∆k; ω) are computed at a wavevector ∆k close to the axis of the lens at (θ ∼ π/2, φ ∼ π/4). (b) shows the collective mode at low energies and the inset shows the single-particle excitations at energies above the gap.
It can be seen that the small lens measurement reproduces quite closely the spectrum S(∆k, ω), in both forward and perpendicular directions. However, compared with the large lens, the small lens leads to the absolute magnitude of the signal being down by nearly a factor of 200 in the forward direction and about a factor of 10 in the perpendicular direction. But even the large lens captures well at least the position of the collective mode peak (at low energies) and single particle peak (at high energies). The collective mode peak is somewhat broadened by the large lens. (These peaks are analyzed shows the renormalized single-particle excitations starting at ω ≈ U/J. All the curves have been normalized to the maximum of the mRPA ≃ 0.30 case. in detail in the next subsections.) Note that unlike the total susceptibility χ(∆k, ∆k; ω), the spectrum S(∆k, ω) contains the dipole matrix elements M g3g4 g2g1 , Eq. (93), that skews the signal towards the forward direction for spinpreserving transitions, and towards the perpendicular direction for spin-exchanging transitions (see Sec. VI B 3). Hence one key finding here is that different lens positions can be used to select for different types of transitions, to separate out the collective modes versus the single particle excitations.

Single particle excitations
In Fig. 23, we plot the RPA spectrum collected by a NA = 0.5 lens in the forward direction [ Fig. 23(a)] and in the perpendicular direction [ Fig. 23(b)]. A broad peak at high energies can be seen in both lens directions. This peak is due to gapped single particle excitations already included by MFT. The form of the MFT susceptibilities of Eqs. (40), (41), (43) and (44) indicates that the single particle excitations spectrum is non-zero only for 2∆ ≤ ω ≤ 2 √ ∆ 2 + 16J 2 . For U ≫ J, the gap can be approximated by 2∆ = 2mU ≈ U as m saturates to 1/2. Indeed, it can readily be seen in Fig. 23 that the gap opens up progressively and more sharply with larger U , and there is signal only between the aforementioned bounds. To compute numerical values, we give a finite value to the infinitesimal imaginary part using δ = 0.07J, in the denominators of the susceptibilities of Eqs. (40), (43) and (44). Thus, individual delta functions coming from the susceptibilities are then spread out into Lorentzian functions, representing the finite frequency resolution of spectral measurements.

Collective modes
The other main feature in Fig. 23(a) and (b) is the sharp peak at low energies. This originates from light scattering off the spin-exchanging transitions that dominate the transverse spin susceptibility χ +− RPA (q, q; ω). In turn, such transitions can be induced by the excitation of gapless collective modes, the spin wave excitations above the AFM ground state (see Sec. IV D and Fig. 5). Notice that the sharp collective mode at low energy separates cleanly in frequency from the single particle excitations only from moderately large U/J onwards. A rough criterion for this separation is that the collective mode bandwidth [∼ 2J H , see Eq. (50) and the text after] should be smaller than the single particle gap 2∆. At small U/J, the single particle gap is small and the collective mode peak merges with the single particle broad peak, as can be seen for the case of U ≈ 2.6J in Fig. 23.
As discussed in Sec. IV, for large U/J, the Hubbard model reduces to the Heisenberg model at low energies. In Fig. 24 we show the angular and frequency dependence of the scattered light along the φ = π/4 line for U = 8.3J (m RPA ≃ 0.30) and compare it with the spin wave dispersion relation of the Heisenberg model, Eq. (50). As a visual aid, the line connecting the maximum of the spectrum for each θ point is also drawn. It is indeed similar to the dispersion relation curve. Such measurement could in principle be realized with a small NA lens scanning the θ direction. Note that this figure only shows the low energy spectrum, cutting off the higher energy single particle excitations shown in Figs. 22 and 23. Despite the similar shape, the RPA derived dispersion occurs for smaller frequencies compared with the Heisenberg model dispersion. This mismatch disappears for values of U 25J.

IX. CONCLUDING REMARKS
We have studied off-resonant imaging of AFM correlations in a two-species fermionic atomic gas in a tightlyconfined 2D optical lattice. The AFM ordering represents a checkerboardlike alternating density pattern of the two species that effectively doubles the lattice periodicity of each spin component. This can be revealed in emerging magnetic Bragg peaks of elastically scattered light when only one spin component is detected.
The density correlations as well as the longitudinal and transverse spin correlations of the atoms are mapped onto the fluctuations of the scattered light where they can be detected in the inelastically scattered light. We have shown how the standard Feynman-Dyson perturbation theory of interacting many-body systems can then be related to the experimental observables of the intensity and the spectrum of the scattered light. Our specific example concerned RPA of the AFM ordering that corresponds to a partial summation of the diagrammatic perturbation series. The general principle could be adapted to other strongly-correlated states, indicating how offresonant imaging can provide a powerful diagnostic tool in interacting ultracold atom systems.  26. Diagrammatic representation of the MFT transverse spin susceptibility matrix, see Eq. (38). The sums over internal momentum k, frequency ν have not been written explicitly.

Mean-field susceptibilities
The MFT susceptibilities correspond to computing Eq. (34) in the matrix form Eq. (38) without any interaction vertices V 0 and V Q . As an example, we calculate the first diagram of the diagonal element of the susceptibility 27. Diagrammatic representation of the MFT densitydensity susceptibility matrix, see Eq. (38). The sums over internal momentum k, frequency ν, and spin have not been written explicitly.
FIG. 28. Diagrammatic representation of the RPA series for the density susceptibility. The first order term is the MFT bubble. The sums over internal momenta k, k ′ , k ′′ , and frequencies ν, ν ′ , ν ′′ have not been written explicitly.ḡ stands for the opposite spin of g.
matrix shown in Fig. 26, The two momenta in the left-hand side of Eq. (A4) are a consequence of the translation invariance with periodic boundary conditions that we have imposed on the system (Sec. V A) and Tr Tr Tr denotes the trace of the matrix. The overall minus sign comes from anticommuting the Nambu spinors, and is an example of the diagram rule for a fermion loop leading to a (−1) factor. In frequency space, this becomes Tr Tr TrG G G 0 ↓ (k, ν)G G G 0 ↑ (k+q, ω+ν) .
This leads to a simple diagonal structure for the RPA susceptibilities for these quantities, as we shall see next. However, the matrix structure of χ χ χ +− (0) does lead to a full matrix equation for its RPA susceptibility.

RPA susceptibilities
For the RPA susceptibilities, the interaction vertices need to be inserted n times for the nth order diagram. . In (c), we show the RPA series for the transverse spin susceptibility. The sums over internal momenta k, k ′ , k ′′ , and frequencies ν, ν ′ , ν ′′ have not been written explicitly.