Quantum Noise Analysis of Spin Systems Realized with Cold Atoms

We consider the use of quantum noise to characterize many-body states of spin systems realized with ultracold atomic systems. These systems offer a wealth of experimental techniques for realizing strongly interacting many-body states in a regime with a large but not macroscropic number of atoms where fluctuations of an observable such as the magnetization are discernable compared to the mean value. The full distribution function is experimentally relevant and encodes high order correlation functions that may distinguish various many-body states. We apply quantum noise analysis to the Ising model in a transverse field and find a distinctive even versus odd splitting in the distribution function for the transverse magnetization that distinguishes between the ordered, critical, and disordered phases. We also discuss experimental issues relevant for applying quantum noise analysis for general spin systems and the specific results obtained for the Ising model.


I. INTRODUCTION
A. Strongly correlated systems of cold atoms A natural regime for ultracold akali gases is when interactions are weak. For typical experimental densities of n ≈ 10 14 atoms per cubic centimeter and scattering lengths of a few nanometers, interparticle distances are two orders of magnitude larger than the scattering length. This regime was the focus of most experimental work done in the first decade since the observation of Bose-Einstein condensation [1,2,3]. Excitement in this field was fueled by the possibility of turning gedankenexperiments into reality. This includes interference of independent condensates [4], creation of atom lasers [4,5], and condensation of multicomponent systems [6,7,8,9,10,11,12]. The Gross-Pitaevski equation captures all essential features of systems in this regime [13,14]. This is a mean-field approximation that only requires the knowledge of a single wavefunction describing a macroscopically occupied state describing the condensate. Conceptually, the many-body aspect is important only in turning quantum coherence into a macroscopic phenomenon. However, solutions to the non-linear equations may be non-trivial with the vortex lattice [15,16,17] being a prime example.
The current focus for the field of ultracold atoms is shifting towards creation and study of systems with strong interactions and correlations. The driving force of this evolution is remarkable experimental progress which provide powerful tools for realizing strongly interacting systems of cold atoms. Magnetically tuned Feshbach resonances (see [18,19] and references therein) can be used to make scattering lengths comparable to interparticle distances. This allowed creation of molecules and observation of pairing in fermionic systems. Confining particles in optical lattices is another powerful tool for amplifying the role of interactions. Atoms in a periodic potential created by a standing wave laser beam can exhibit both enhanced interactions and suppressed kinetic energy. The relative strength of the two energies can be tuned by changing the intensity of the laser beam. This technique was used to observe the superfluid to Mott transition in which the behavior of atoms changed from wavelike in the superfluid phase to particle like in the Mott phase [20]. Reduced dimensionality is another useful tool for making strongly correlated systems. The strength of an effective one-dimensional interaction can be changed by tuning either the particle density, transverse confinement, or scattering length. This approach was used in several experiments that explored the crossover from weakly interacting quasicondensate to the Tonks gas of hard core bosons [21,22]. There are many other promising techinques for obtaining strongly interacting systems of quantum degenerate particles including rotating condensates [23], ultracold polar molecules (see [24] and references therein), and Rydberg atoms [25,26].
One of the primary motivations for studying strongly correlated systems of cold atoms is the expectation that they will realize important and not yet fully understood quantum many body states. The target list of modern experiments includes the Fulde-Ferrel-Larkin-Ovchinnikov phase of unbalanced fermion mixtures with attraction [27,28], supersolid phase [29,30,31], quantum magnetic phases of atoms in optical lattices (ordered or spin liquids) [32,33,34,35], fermions with d-wave pairing in lattice systems with repulsive interactions [36], and quantum Hall phases [37,38,39]. An important open question is detection and characterization of these states.
Analysis of correlation functions is a common way of characterizing many-body states in condensed matter physics. Neutron and light scattering experiments can be used to measure spin and charge response functions. Spin and charge density wave orders give rise to additional δ-function peaks in the static correlation functions, hence they immediately show up as new Bragg peaks in elastic neutron scattering experiments. Conductivity measurements at finite frequency (microwave, infrared, optical, etc.) provide access to current-current correlation functions [40]. They are a useful tool for detecting states with a gap in the quasiparticle spectrum, such as superconducting and various density wave states. Photoemission experiments measure electron spectral functions and can be used to study excitation spectrum or observe electron fractionalization in one dimensional systems [41,42,43]. So what we want for strongly correlated systems of cold atoms is to have a tool box for measuring correlation functions.
Certain experimental techniques used in cold atoms provide natural analogues of condensed matter experiments. For example, Bragg spectroscopy is similar to light and neutron scattering in solid state physics and can be used to analyze the dynamic response functions [44]. RF spectroscopy in cold atom systems is analogous to finite frequency conductivity measurements [45,46,47,48]. There are also several proposals for designing experiments that mimic solid state methods [49,50]. Additional techniques such as time of flight and interference experiments [1,2] or the possibility of single atom detection [51,52] are also available which are not feasible with traditional condensed matter systems. The richness of physical phenomena that we expect from strongly correlated systems of cold atoms provide strong motivation to think about new methods and approaches that would provide "smoking gun" signatures of various many-body phases.

B. From noise analysis to correlation functions in interacting cold atoms systems
The usual philosophy of performing measurements in physics is to minimize the noise by averaging over many experiments. However, there are notable exceptions. In mesoscopic solid state physics and quantum optics, one is often interested in the noise rather than the average signal. For example, shot noise was used to measure fractional charge and statistics of excitations in fractional quantum Hall systems [53,54]. Quantum noise in photon counting experiments was used to identify non-classical states of light [55].
The philosophy of noise analysis was implemented in several recent types of experiments with cold atoms. These experiments analyzed quantum noise in time of flight images of atoms released from an optical lattice [52,56,57,58,59] (see Ref. [60] for theoretical analysis). Free expansion during the time of flight maps momentum distribution of atoms inside the interacting system into the density distribution after the expansion. Correlation functions ρ(r)ρ(r ′ ) of the images in real space contains information about the n k n k ′ correlation function of the interacting system in momentum space. Note that n k = c † k c k is the occupation number of momentum states and not the density operator at wavevector k. Hence n k n k ′ represents a non-local correlation function for the original interacting system.
What can one learn from n k n k ′ correlation functions? The simple answer is the periodiciy of the original system. Quasimomentum (i.e. lattice momentum) is a good quantum number in a periodic system. Quasimomentum differs from the physical momentum by Bragg reflections by reciprocal lattice vectors. This leads bunching for bosons and peaks in n k n k ′ when the two momenta differ by a reciprocal lattice vector k − k ′ = G. Analogously, fermions exhibit antibunching. Such experiments have been done for spinless bosons [56,57,59] as well as fermions [58] and provided one of the most elegant realizations of Hanbury-Brown-Twiss experiments with neutral atoms. When such experiments are performed with mixtures, we can expect even more exciting results. They should provide an ideal probe for antiferromagnetism or static density wave, since both phenomena increase the size of the unit cell and give rise to additional reciprocal lattice vectors [60]. It is worth pointing out that fluctuating orders, like in Luttinger liquids, can also be observed using this method [61]. Another beautiful application of the noise analysis in time of flight experiments was the measurement of fermion pairing on the BEC side of the Feshbach resonance [58]. These experiments demonstrated correlations in the occupation numbers for fermions with momenta k and −k, which demonstrated the existence of a Cooper pair condensate at zero momentum.
A different method for analyzing correlation functions in interacting many-body systems is based on interference experiments between two independent condensates. It may not be immediately obvious why this measurement corresponds to the noise analysis. In every shot we should find an intereference pattern. However the position of interference fringes is unpredictable from shot to shot. If we average over many experiments we will find that the interference pattern washes out completely. We can put this in a more formal mathematical language: Let A be the quantum mechanical operator that describes the amplitude of interference fringes [85], we have A = 0 but finite |A| 2 . What is interesting for our purposes is the value of |A| 2 itself, or the related quantity of the fringe contrast. When individual condensates are not fluctuating, we expect fringes with perfect contrast. When there are strong thermal or quantum fluctuations we expect suppression of the fringe amplitude. Theoretical analysis performed in Ref. [62] showed that the amplitude of interference fringes can be related to the instantaneous two point correlation function: |A| 2 = Ω Ω G 2 (r), where G(r) = c(r)c † (0) and Ω is the imaging area. By analyzing scaling of |A| 2 with the size of the images region, Hadzibabic et al. recently demonstrated the Berezinskii-Kosterlitz-Thouless transition in two dimensional condensates [63].
The common philosophy of HBT experiments with atoms [52,56,57,58,59,60,61,62,63] and analysis of the average contrast in interference of independent condensates (IIC) experiments is the focus on second order coherence. This means that in the analysis of absorption images, we are not looking for peaks in ρ(r) but in ρ(r)ρ(r ′ ) [86]. An interesting question to ask is whether one can use quantum noise measurements to analyze higher order coherence. In the case of IIC experiments this problem was considered in Refs. [62,64]. The basic idea is that we do not only measure the average value of the amplitude of interference fringes but also their fluctuations. And it is straightforward to show that that higher moments of |A| 2 correspond to higher order correlation functions. In the case of one dimensional Bose systems with interactions it is possible to find the explicit form of the distribution function of |A| 2 . For weakly interacting bosons phase fluctuations within condensates are small and one finds the Gumbel distribution of |A| 2 with a small width. When interactions between atoms become so strong that the system approaches the limit of hard core bosons, the distribution function becomes Poissonian.

C. Quantum magnetism with cold atoms in optical lattices
The main subject of this paper is to extend the approach of quantum noise analysis to the study of higher order correlations in quantum magnetic systems realized using cold atoms. Several approaches have been discussed recently for creating magnetic systems with cold atoms (see Ref. [35] for a review). Our proposal is sufficiently general so we expect that it can be applicable to all of them. The main idea can be understood in the following example. Consider a measurement of the α component of spin magnetization in a one dimensional spin systems (1.1) The average magnetization can be obtained by averaging over many experiments but individual measurements will exhibit quantum noise. Each measurement will give a value of M α that is different from the average M α . It is easy to see that the strength of quantum fluctuations is controlled by correlation functions of the spin operators. For example So the width of the distribution is determined by two point correlation functions. By the same argument higher moments provide a measure of higher order correlation functions The knowledge of high moments is contained in the distribution function of M α . So if we know the distribution function of P (M α ) we obtain information about correlation functions of arbitrarily high order. In real experiments, of course, technical noise and a finite number of experiments will limit the number of moments that can be reliably extracted.
Generalization of this idea to other measurements is straightforward: one should measures not just the average value of some operator but its entire distribution function. The distribution function contains information about high moments which can be related to high order correlation functions. This provides valuable characterization of a many-body state. It is important that measurements are done on a mesoscopic system, hence fluctuations will not be neglibly small compared to the mean value. For macroscopic solid state samples this is not easy to achieve, but for cold atoms systems with tens or hundreds of thousands of atoms this is an experimentally relevant regime. Another advantage of the cold atoms systems is their nearly complete isolation from the environment. We point out that Ref.
[65] considered a related idea of measuring the number of particles in a mesoscopic part of a fermionic system with pairing. They showed that fluctuation in the total number of particles will exhibit interesting evolution as one tunes the system from the BCS to the BEC regime.
In the main part of this paper contained in Sec. II we consider a specific example of the magnetization distribution for the Ising model in a transverse field. We are not the first ones to address this question. Ref. [66] study the distribution functions for both the longitudinal and transverse magnetizations of this model with analytical and numerical techniques. We focus primarily on the transverse magnetization but provide a more detailed analysis of its global structure. In particular, we study how higher order moments encode the different correlations in the ordered, critical, and disordered phases and give a simple physical picture in terms of confinement versus proliferation of domain walls. Sec. III discusses experimental issues involved in applying quantum noise analysis to ultracold atoms including those specific to the results we obtained for the transverse Ising model. We summarize our discussion in the conclusion of Sec. IV. The appendices contain mathematical details of the discussion in Sec. II.

II. MAGNETIZATION DISTRIBUTION IN THE TRANSVERSE ISING MODEL
The one-dimensional Ising model in a transverse field provides a textbook example of a symmetry-breaking quantum phase transition [67]. It is described by the Hamiltonian where S α (i) are spin-1/2 operators at site i. Here g is a dimensionless coupling which sets the strength of the transverse field and J sets the overall energy scale. We take periodic boundary conditions and N is the number of sites. The competition between the transverse field and the Ising term of the Hamiltonian drives a quantum phase transition at g c = 1 from an ordered phase g < g c with spontaneous magnetization along the x direction to a disordered phase g > g c with zero magnetization along x. We are interested in the distribution function for the magnetization M α given by where δ(x) is the delta function and brackets denote ground state expectation values. The global structure of P (M x ) describing the order parameter follows from low-order moments of the longitudinal magnetization M x and symmetry breaking arguments. Away from the critical point g = g c , exact results on S x (i)S x (j) [68,69] show the spin correlation length and the second moment M 2 x are finite. Physically, this means spins fluctuate independently on length scales larger than the correlation length. Thus, for a sum of a large number of spins, the corresponding distribution should be approximately gaussian. In the broken symmetry phase g < g c , the distribution should have two peaks describing the formation of spontaneous magnetization for M x . On the other hand, the unbroken symmetry phase g > g c should only have one peak centered at zero. At the critical point g = g c , the power-law decay of S x (i)S x (j) ∼ |i − j| −1/4 gives a divergent contribution to M 2 x and a non-trivial distribution function.
More analysis is needed for the transverse magnetization M z where symmetry breaking arguments do not apply. Exact results [68,70] for S z (i)S z (j) show the second moment M 2 z is finite even at the critical point. Analogously, P (M z ) is also gaussian. However, additional aspects of the global structure such as the number of peaks and other features are encoded in higher order moments.
Instead of analyzing individual moments, we study the entire distribution P (M z ) directly. We find its global structure reflects the amount of disorder in the system. This is in contrast to P (M x ) which reflects the amount of order. We show P (M z ) is asymptotically gaussian for a large number of spins with a characteristic splitting between even and odd magnetization that differs between the phases. A duality mapping shows this behavior has a physical interpretation in terms of proliferation, criticality, and confinement of dual domain walls when g < g c , g = g c , and g > g c , respectively. The domain wall physics is particularly clear in the small and large g regimes where P (M z ) describes the qualitatively different counting statistcs of pairs and single domain walls, respectively. In section II A we calculate P (M z ) via fermionization. We obtain and discuss its global structure in section II B, focusing on the even versus odd splitting and its interpretation using a duality mapping. Appendix A discusses the use of Toeplitz determinants employed in the calcuation. In addition, appendix B derives the small and large g behavior in more detail. Finally, Appendix C describes the saddle-point analysis used to obtain the asymptotic gaussian behavior.

A. P (Mz) Distribution Function
Jordan-Wigner fermionization (see Ref. [71]) maps spins in one dimension onto spinless fermions where σ α (i) = 2S α (i) are the Pauli matrices and σ ± = σ x ± iσ y are the corresponding raising and lowering operators. The spin Hamiltonian of Eq. 2.1 maps to free fermions where c k are in momentum space. Instead of studying M z directly, the analysis is simpler for the quantity which is simply related to both M z and the fermion number. Recall N is the number of spins in the entire chain and compare this to n which is the number of spins in a block embedded in this system. We will take N → ∞ first to describe the thermodynamic limit and then n → ∞. Physically, we are studying the local magnetization M z of a large block of spins embedded in a much larger system. This has important consequences in the interpretation of the results and experimental detection. We will calculate the generating function for ρ n defined as where the brackets denote ground state expectation values in the thermodynamic limit N → ∞. The corresponding distribution function is given by the fourier transform of the generating function which is related to P (M z ) by the following We first define the following linear combinations of fermionic operators and obtain for the generating function where we have used the operator identity e iλ(1−c † i ci) = A i B i (λ). Using Wick's theorem, the expectation value in Eq. 2.12 can be reduced to products of the following two-point functions where the Bogoliubov angle comes from diagonalizing the Hamiltonian in Eq. 2.5 using a Bogoliubov transformation. Notice A i A j = δ ij . However, the contractions of this form entering χ n (λ) vanish because they always occur with i = j. Moreover, any full contraction containing B i B j also vanishes since it must contain some A i A j contraction. This implies the only full contractions that contribute are those that only pair A i with B j . Taking into account the fermionic signs, these full contractions organize into a determinant of the form where f n (z) are the Fourier coefficients of the function Notice that the coefficients are given by f n (e iλ ) = A i B j (λ) and the function f (z, x) is the integrand in Eq. 2.15 with z = e iλ and x = e ik . The above result gives an explicit expression for the generating function. Here we can make some general comments on the structure of χ n (λ). By diagonalizing the matrix inside the determinant, we see that the generating function is of the form Here µ i are the eigenvalues of the n × n matrix M ij formed from the Fourier coefficients of the Bogoliubov angle This can be rewritten in the more suggestive form as a product of generating functions χ B (λ, p). We recognize these generating functions as that of a Bernoulli random variable with probability p of obtaining one and 1 − p of obtaining zero. This implies ρ n is the sum of independent Bernoulli variables. This is not surprising since ρ n is directly related to the fermion number and each fermionic level can be either occupied or unoccupied. However, even though the matrix M ij is real, the eigenvalues µ i are only required to occur in complex conjugate pairs. In fact, we find the µ i are generically complex and satisfy |µ n | ≤ 1. This means that χ B (λ, p i ) cannot be interpreted as a classical Bernoulli distribution. However, a pair of conjugate eigenvalues gives allowing an interpretation as a classical distribution with the coefficient of e imλ giving the probability of obtaining m.
Notice that the interference of the complex eigenvalues µ i can lead to suppression of odd versus even values. In the fermionic description, this is a result of BCS-type pairing in the Hamiltonian of Eq. 2.5. The anomalous term a k a −k and its hermitian conjugate describe pair creation and annhilation processes which might favor even occupation. We note that the appearance of independent Bernoulli variables for the fermion number distribution has been studied in the mathematical literature (see Ref. [72]) as determinantal point processes. However, it appears the cases that have been studied correspond to states with no BCS-type pairing and purely real eigenvalues µ i . These heuristic considerations take into account only one or two fermionic levels. Determining the global structure of P n (m) is a many-body problem involving a large number of fermionic levels which we address in the next section.

B. Global structure
Before discussing the global structure of P (M z ), we first argue that it gives a measure of the amount of disorder in the system. Physically this is not surprising since S z (i) tends to disorder spins along the x direction by inducing tunneling between the two polarizations. This competes with the Ising term S x (i)S x (i + 1) which favors ordering of spins along the x direction. A duality transformation [73] makes this connection more concrete. This transformation maps the spin Hamiltonain of Eq. 2.1 to a spin Hamiltonian of the same form but with the coupling constant exchanged (2.25) The T α (i) operators describe dual spins that obey the same commutation relations as S α (i) and are given by The center and width of the peak tend to zero for g → ∞ but approach a constant for g → 0.
in terms of the S α (i) operators. Under this transformation the operator ρ n maps to where the term under the summation has eigenvalue 0 for two dual spins aligned on the x direction and eigenvalue 0 for anti-aligned. This implies that ρ n is the number operator for dual domain walls and P n (m) is its probability distribution.
In particular, the fermionic operators c i (c † i ) create (annhilate) dual domain walls. For the BCS-type Hamiltonian of Eq. 2.5, Jg acts as the chemical potential for free fermions with dispersion −J cos(k) with the anomalous terms describing pair creation and annhilation. When g is small, the chemical potential crosses the band and pair fluctuations cost little energy leading to dual domain wall proliferation. When g is large, the chemical potential is below the bottom of the band and dual domain walls are confined due to energy cost. Note this is the exact opposite of the physical domain walls for the S x (i) spins but the underlying physics of confinement versus proliferation [74] is the same.
Analysis of P n (m) shows how the distribution function reflects the above physics Numerically, Eq. 2.17 allows for efficient calculation of the distribution function p n (m). The determinant representation gives χ n (λ) as a polynomial of degree n in e iλ and the coefficient of e imλ can be extracted via a computer algebra system to give P n (m). Analytical insight is also possible in the large-n limit through the use of Toeplitz determinant asymptotics outlined in Appendix A. We find that the generating function has the following asymptotic behavior χ n (λ) ∼ exp c 0 (e iλ )n + c 1 (e iλ ) log n + c 2 (e iλ ) (2.29) where the coefficients c i (z) are explicitly known functions. We first consider the small and large g regimes. We carry out the expansions in Appendix B and find for the generating function where we recognize the form of the generating functions as that of a binomial random variable for g ≪ g c and a sum of two Poissonian random variables for g ≫ g c . However, notice in the g ≫ g c regime that one poissonian variable depends on e 2iλ whereas the other depends on e iλ . This implies the operator ρ n behaves as where B(n, p) is a binomial random variable with n trials and probability of success p and P (λ) is a Poissonian random variable with parameter λ. Explicit expressions for P n (m) for these two regimes are given in Appendix B. We plot the distribution P n (m) for g = 0.5 and g = 2 in Fig. 1. Notice the distinct splitting between even and odd values of m for g < g c that is absent for g > g c . The qualitatively different behavior has a simple physical interpretation in terms of proliferation versus confinement of dual domain walls. We have argued M z simply counts the number of dual domain walls. When they proliferate for small g, individual domain walls are free. Each site contributes independently with probability P occ = 1/2 − g/4 of being occupied and P empty = 1 − P occ of being empty. This leads to the binomial counting statistics that does not distinguish between even and odd total occupation. In contrast, the confined regime for g > g c has dual domain walls tightly bound into pairs that occur with low probability P pair = 1/16g 2 . In the dilute limit, the counting statistics of pairs is Poissonian. However, pairs in the bulk contribute a factor of two to the number of dual domain walls with a probability (n − 1)P pair that scales with n whereas pairs at the boundary contribute a factor of one to the total number with a probability 2P pair since there are two boundaries. The differing contriubtions of pairs in the bulk and boundary leads to the splitting between even and odd total occupation. A schematic illustration of counting dual domain walls in these two regimes is illustrated in Fig. 2.
The above analysis relies on an expansion in the small parameter g or g −1 and identifying the resulting form of χ n (λ). In order to study the evolution from the small to large g regimes, we use a large-n saddle-point integration for P n (m) which works for all g. This comes at the expense of obtaining only the asymptotic gaussian behavior and is carried out in Appendix C. The gaussian nature of the resulting distribution can be established two ways. We have already argued that ρ n is the sum of n independent Bernoulli-like random variables with finite variance which means P n (m) is guassian by the central limit theorem. In addition, we can examine the Taylor series log χ n (λ) = ∞ p=0 λ p p! ρ p n c .
(2.32) whose coefficients ρ p n c give the cumulants or connected moments. We find all cumulants are finite and proportional to n at leading order. By considering the ratio we conclude that the first two cumulants completely characterize the distribution for large-n which implies that P n (m) is gaussian. This is not in contradiction to the small and large g regimes since the resulting distributions are also gaussian for sufficiently large n. The absence or presence of an even versus odd splitting characterizes the small and large g regimes, respectively. This aspect of the global structure of P n (m) reflects higher order moments which encode the underlying correlations of the model and appear in the large-n gaussian distributions. In particular, we find which is of the form of either a single gaussian or the sum of two gaussians. The (−1) m term gives the even versus odd splitting which has a sharp transition and critical behavior at g = g c . This illustrates clearly that the quantum phase transition manifests itself in higher order correlations of M z and is accessible via the distribution function P (M z ). Having disucssed the physical origin of the this splitting, we now show how it also reflects correlations of the order parameter or longitudinal spins S x (i). At first this is suprising since P n (m) describes correlations of the transverse spins S z (i). However, recall the S z (i) correspond to domain walls for T x (i) which are related to S x (i) via duality. Explicitly, we obtain From Eq. 2.28 This gives the two-point function of the dual spins T x (0)T x (n) as a function of ρ n implying that P n (m) can be used to compute its average in the ground state The interpreation of this equation is also illustrated in Fig. 2. An even (odd) number of dual domain walls between two sites gives a spin configuration with spins aligned (anti-aligned) along the x direction. By using the gaussian form of P n (m) in Eq. 2.34 and passing from summation to integration we find where we have dropped the contribution from terms entering with (−1) m in the summation. Here C = 2 1/12 e 1/4 A −13 with A being Glaisher's constant. This agrees with known results on the spin correlation function S x (0)S x (n) [68,69] after applying a duality transformation g → g −1 . Why the agreement is exact is discussed in Appendix C.
We now discuss the mean µ and variance σ 2 of the guassian distributions in Eq. 2.34. They are expressed in terms of the Bogoliubov angle by We plot the mean µ and standard deviation σ in Fig. 3. As g → 0, the distribution has a finite width σ and a finite mean µ. In contrast, as g → ∞, both σ and µ go to zero. This behavior can also be seen in Fig. 1. We plot P n (m) in rescaled variables in Fig. 4. The asymptotic convergence to gaussian distributions is quite good even at the critical point g = g c for sufficiently large n. In addition, the distinctive even versus odd splitting is absent for g < g c , present for g > g c , and decays algebraically for g = g c .
Here we compare our results with that of Ref. [66]. Our new contribution is the demonstration and analysis of the even versus odd splitting in P n (m). Fig. 2 shows that this is fundamentally a boundary effect of a local block of where . . . c denotes the cumulant. The solid lines indicate the asymptotic gaussian form of pn(m). Notice the convergence to guassian behavior even at g = gc and the splitting between even and odd m for g ≥ gc.
spins. Thus it is necessary to separate the system size N from the block size n to see this effect. In contrast, Ref. [66] implicitly set N = n in the beginning of their analysis. At a more technical level, we find an approach based on mapping the problem to Toeplitz determinants allows for separation of the length scales N and n. This allows us to obtain the first finite size correction to the leading gaussian behavior to P n (m) as discussed in Appendix A which turns out to be crucial in obtaining the splitting. In addition, this approach also facilitates study of the small and large g regimes where the underlying physics is most evident. As another technical note, Ref. [66] also obtains the asymptotic gaussian behavior via a saddle-point analysis, but we show that there is an additional saddle-point which gives rise to the splitting and is discussed in Appendix C. Having discussed one application of quantum noise analysis in a specific spin system, we now turn our attention to a discussion of general issues in experimental realizations.

III. EXPERIMENTAL ISSUES
Although the use of ultracold atoms in experimental studies of quantum many-body systems has seen great progress, techniques based on quantum noise analysis are just becoming within reach. The fundamental difficulties arise from the need to collect accurate statistics on noise fluctuations of an observable, ideally to high orders. Here we discuss some proposals for experiments and the relevant issues for quantum noise analysis.
Optical lattices offer one of the most promising avenues for engineering many-body systems in ultracold atoms including spin systems [35] but they may pose problems for noise analysis. This is especially true for low-dimensional systems. They have already been realized in optical lattices [21,22,56] but occur in a large array of up to thousands of one-or two-dimensional systems within each sample. A single measurement such as an absorption image over the entire sample corresponds essentially to averaging over the large array which supresses fluctuations. Low-order statistics such as second-order coherence [56,57,59] are still accessible. However higher order noise fluctuations may be difficult to detect as a result of the intrinsic averaging that occurs within each sample.
Accessing higher order noise may be possible through realizations of single or small numbers of low-dimensional systems where fluctuations are not negligble compared to the average. Such an approach has already been used to study the Berezinskii-Kosterlitz-Thouless transition [63]. In this experiment, the quantum noise in the interference amplitude of just a pair of two-dimensional quasi-condensates was used to detect signatures of the transition. For one-dimensional systems, surface microtraps allow for tight confinement in a waveguide geometry [75,76,77] and optical box traps have recently demonstrated single one-dimensional condensates [78].
In addition to preparation, experiments will also have to advance in measurement techniques. The ideal measurement of quantum noise for magnetization would be state-selective in situ imaging with single atom precision to obtain accurate high order statistics and high resolution to resolve spatial correlations. This has not yet been achieved with optical absorption imaging but may be possible through proposed scanning tunneling microscropy of ultracold atoms [79]. However, single atom counting has been demonstrated in an atom laser setup [51]. However, this is more suited to the study of temporal correlations rather than the spatial correlations describing noise in magnetization. Detection of single atoms has also been used in Hanbury-Brown-Twiss experiments [52] but with imaging after expansion which resolves momentum correlations and not in situ which probes real space correlations. An alternative measurement involves imprinting the quantum noise for magnetization onto a weakly coupled laser probe [80]. In this case, the quantum noise of the laser itself should be less than that of the magnetization and may make features at the single atom level difficult to detect.
The above discussion of general issues for quantum noise experiments also applies to the transverse Ising model we considered. However, resolving the even versus odd splitting in the distribution function requires additional considerations. As shown in Fig. 2, distinguishing between even and odd magnetization requires a sharp boundary between the block of spins undergoing measurement and the rest of the spin chain. This requirement does not pose difficulties for measurement via high resolution in situ imaging. However, for techniques that couple to the entire sample, RF cutting (see Ref. [81]) may be needed to isolate the region of interest beforehand. Sensitivity at the single atom level will also be crucial for detecting the even versus odd splitting. We consider the feasibility of such a measurement in the absence of in situ imaging by analyzing the proposal of Ref. [80] for coupling the magnetization to off-resonant coherent laser light. For our purposes, the main result relates the output quadrature of light X out to the input quadrature X in and the magnetization M z as where κ is the atom-photon interaction and n a is the number of atoms. Here the input quadrature component X i n has been rescaled by the number of photons n p and describes coherent light with noise ∆X in = 1/ √ 2 The quantum noise of the magnetization is thus imprinted on top of the inherent photon shot noise of light. Discerning gross features of the magnetization distribution only requires fluctuations in X out due to quantum noise in the collective magnetization be greater than fluctuations due to photon shot noise. The width of the magnetization distribution gives an estimate ∆M z ∼ √ n a σ where σ is shown in Fig. 3 and yields the condition This enables resolving the smooth distribution shown in Fig. 1 in the regime with no splitting whereas it only gives the average of the two separate distributions of the same figure in the regime with splitting. Detecting this splitting requires X out fluctuations due to single spins rise above photon shot noise and yields a much more stringent condition κ > n a 2 . (3.3) Eq. 3.1 describing the coupling of magnetization to light is valid when decoherence due to spontaneous emission is small and for low photon absorption. This implies the atomic depumping η = n p σ r Γ 2 /A∆ 2 and photon absorption ǫ = n a σ r Γ 2 /A∆ 2 are small. Here n p is the number of photons, σ r is the resonant scattering cross section which is distint from σ the width of the magnetization distribution, Γ is the spontaneous decay rate, A is the illuminated cross sectional area, and ∆ is the detuning. The atom-photon interaction is given by κ = √ n a n p σ r Γ/A∆. For given η, ǫ ≪ 1, the condition for resolving collective magnetization noise in Eq. 3.2 becomes which can be satisfied for sufficiently large detuning. However, Eq. 3.3 for resolving single spin noise becomes which is likely an unphysical regime since since σ r ∼ λ 2 where λ is the laser wavelength and η is a small parameter. This states that the photon shot noise for coherent light overwhelms the single spin noise. Using squeezed states of light with smaller quadrature fluctuations may allow detection of single spin noise but a more detailed analysis would be necessary.

IV. CONCLUSION
Ultracold atoms offer the prospect of experimentally realizing quantum many-body systems in a strongly interacting regime with advancements such as optical lattices, Feshbach resonances, and confinement to low dimensionality. Noise analysis, or the study of fluctuations of an operator potentially to high orders instead of just the average, offers an avenue to access correlation functions and to characterize many-body states. We discussed a proposal for using noise analysis in quantum spin systems realized with ultracold atoms. By studying the entire distribution function for the magnetization of a mesoscopic spin system, the behavior of high order correlation functions can be extracted. We analyzed the Ising spin chain in a transverse field in detail, focusing on the transverse magnetization. The resulting distribution function possesses a distinct splitting between even and odd values of the magnetization. This splitting behaves differently in ordered, critical, and disorderd phases and has a simple interpretation in terms of confinement versus proliferation of domain walls. Finally, we discussed pertinent issues for implementing quantum noise analysis through experiments in the near future.
This replacement corresponds to noticing f i−j (z) is translationally invariant in the sense that it depends only on i − j and evaluating the trace in momentum space. The subleading terms correspond to finite size corrections in n.
The functions c i (z) depend on the analytic structure of log f (z, x) and differ for g < g c , g = g c , and g > g c . It will be useful to define the following functionals acting on functions of x where the contour integrals are along the unit circle and notice h y [f (x)] depends on an additional parameter y. If there are singularities on the unit circle then the integrals are in the principal value sense. First we consider g > g c . For all z, f (z, x) is continuous on |x| = 1 with zero winding number. We obtain (A.8) Next we consider g < g c . For Re[z] ≥ 0 (Re[z] < 0), log f (z, x) is non-singular with winding number zero (one). We obtain where |ξ| < 1 corresponds to the singularity of log f (z, x) closest to the unit circle given by .12) and the square root is positive for positive real arugments. The final case we consider is g = g c . Here, f (z, x) is discontinuous at x = −1 for all z. This singularity can be written in power-law form where the exponent is given by Notice log f (z, x)/g(z, x) is continuous. We obtain where G(x) is the Barnes-G function satisfying G(x + 1) = Γ(x)G(x + 1) with Γ(x) the gamma function. We now turn to the small and large g expansions of the generating function χ n (λ) in Eq. 2.30. For g ≪ g c , we use the determinant representation in Eq. 2.17 and first expand the matrix elements f n (z) in powers of g. We find The determinant det[f i−j (z)] is lower triangular and given by the product of the diagonal matrix elements where Γ(n) is the gamma function. The distribution P n (m) is given by Eq. B.3 with p = (2 − g)/4 which we use to plot Fig. 1.
For g ≫ g c , we analyze the determinant of Eq. 2.30 using the large n asymptotic of Eq. 2.29 as discussed in Appendix A. Instead of the integral representation employed in Eq. A.6 it will be easier to use the equivalent representation in terms of the fourier coefficients of log f (z, x) given by We expand the coefficients h n (z) in powers of g and find Using Eq. B.5, this gives the generating function to O(g −3 ) as which is that of the sum of two Poissonian random variables as in Eq. 2.30. However, notice one term enters with e 2iλ and the other enters with e iλ meaning we identify ρ n ∼ 2P ((n − 1)/16g 2 ) + P (1/8g 2 ) in Eq. 2.31. These sums can be evaluated in terms of the hypergeometric function 1 F 1 (a; b; z) The distribution P n (m) is given by Eq. B.12 with λ 1 = 1/8g 2 , λ 2 = (n − 1)/16g 2 which we use to plot Fig. 1.

APPENDIX C: SADDLE POINT ANALYSIS
In this appendix, we obtain the distribution P n (m) by evaluating the integral of Eq. 2.8 in the saddle-point approximation. We rewrite the integral as a contour integral over the unit circle |z| = 1 with the asymptotic form of χ n (λ) in Eq. 2.29 to obtain P n (m) ∼ 1 2πi dz z m+1 exp [c 0 (z)n + c 1 (z) log n + c 2 (z)] . (C.1) The leading contribution to integral for large n is controlled by the term z −m exp[c 0 (z)n] which has saddle points at where primes denote differentiation with respect to z. For the function c 0 (z) obtained in Appendix A, we find there are two saddle points. One is located on the positive real z axis and the other is on the negative real z axis which we denote as z ± respectively. In addition, the |z| = 1 contours can be deformed to steepest descent contours and we obtain in the saddle-point approximation P n (m) ∼ ± ± z −(m+1) ± exp [c 0 (z ± )n + c 1 (z ± ) log n + c 2 (z ± )] 2π[(m + 1)/z 2 ± + nc ′′ 0 (z ± )] .

(C.3)
We have already argued that the in the large-n limit, P n (m) should be gaussian in the discussion of Eq. 2.33. In general, the gaussian form of P n (m) is good near the peak of the distribution. Thus, we will perform an expansion of Eq. C.3 about the peak located at nµ ± and the corresponding saddle-points at z 0 ± . The two quantities are related by the zeroth order saddle-point condition and we have dropped the ± subscript for clarity. We then substitute m = nµ + nδm (C.5) z = z 0 + δz (C. 6) into the full saddle-point equation of Eq. C.2 and solve order by order to obtain δz = a 1 δm + a 2 δm 2 (C.7) where the coefficients a i are given by a 1 = 1 c ′ 0 (z 0 ) + z 0 c ′′ 0 (z 0 ) (C.8) (C.9) Evaluating Eq. C.3 up to second order in δm for the leading term z −m exp[c 0 (z)n] and zeroth order for the other terms yields P n (m) ∼ ± ±1 z 0 ± e −(m−nµ±) 2 /2nσ 2 ± −m log z 0 ± 2πnσ 2 ± exp c 0 (z 0 ± )n + c 1 (z 0 ± ) log n + c 2 (z 0 ± ) (C. 10) where σ 2 ± is given by and the ± subscripts have been suppressed. Since nµ should correspond to the peak of the distribution, we see that z 0 is determined by the condition that the term linear in m of the exponent in Eq. C.10 should have a real part of zero to ensure that it does not shift the peak. This implies that z 0 ± = ±1. At this point, we have to check whether the z 0 ± saddle-points contribute to the leading term z −m exp[c 0 (z)n] at the same order. The functions c 0 (z) are given in Appendix A and from Eq. A.6 and Eq. A.15 we find that c 0 (±1) = 0 which means that both saddle-points should be included for g ≥ g c . However, from Eq. A.9, we see that c 0 (+1) > c 0 (−1) due to the additional contribution of log(−ξ) of Eq. A.12 with |ξ| < 1. Thus the contribution from z 0 − should be dropped for g < g c .
Evaluating the remaining expressions in Eq. C.10 using the functions c i (z) obtained in Appendix A gives the gaussian form of Eq. 2.34. In addition we find that µ ± = µ given by Eq. C.4 corresponds to the mean and σ 2 ± = σ 2 given by Eq. C.11 corresponds to the variance and are given by the expressions in Eq. 2.38. An enlightening way of writing Eq. C.10 is given by where χ n (λ) with λ = 0, π correspond to contributions from the zeroth order saddle points z = e iλ = +1, −1, respectively. The contribution from the z = +1 saddle point is dictated by normalization of the distribution function is determined by the T x (0)T x (n) correlation function. Since the asymptotic expression we obtained for χ n (λ) is exact for any λ and the even versus odd splitting is given precisely by χ n (π), we see that the splitting gives the exact asymptotic behavior of T x (0)T x (n) .