Non-Gaussian distribution of collective operators in quantum spin chains

We numerically analyse the behavior of the full distribution of collective observables in quantum spin chains. While most of previous studies of quantum critical phenomena are limited to the first moments, here we demonstrate how quantum fluctuations at criticality lead to highly non-Gaussian distributions thus violating the central limit theorem. Interestingly, we show that the distributions for different system sizes collapse after scaling on the same curve for a wide range of transitions: first and second order quantum transitions and transitions of the Berezinskii-Kosterlitz-Thouless type. We propose and carefully analyse the feasibility of an experimental reconstruction of the distribution using light-matter interfaces for atoms in optical lattices or in optical resonators.


Introduction
The understanding of phase transitions and critical phenomena lies at the very heart of condensed-matter physics [1]. Standard quantum phase transitions, i.e. those following Landau's theory, are signalled by the onset of a local order parameter when local symmetries are broken [2]. This theory sets the mean value of the order parameter at the center of stage. However, going beyond the first order moment reveals interesting information about the many-body state without performing a full state tomography. For instance: in experiments with ultracold atoms, noise correlations reveal antiferromagnetic ordering [3,4]; variances of collective operators in the form of structure factors allow to distinguish between quantum phases such as the superfluid and Mott insulator [5], and to detect antiferromagnetic or crystal ordering [6] and collective entanglement [7,8]; the kurtosis, related to higher order moments, provides information about quasiparticle dynamics and their interactions [9]; the full probability distribution of contrast in interference experiments can reveal strongly correlated atomic states [10,11]; and high-order correlation functions are needed to describe properly the physical outcome in the single-shots of several experiments of quantum many-body dynamics [12]. Moreover, in recent experiments with trapped ions, the full-counting statistics of spin fluctuations allows for the detection of entangled, over-squeezed states [13].
In this paper we go beyond the first moments of the order parameter and analyse the full probability distribution function (PDF) of collective operators in spin chains. The PDF of the order parameter plays a central role in statistical mechanics. When the correlation length is finite, deep in an ordered phase, the system can be regarded as the sum of independent subblocks of finite length and the central limit theorem leads to a Gaussian PDF, in agreement with Landau's paradigm. Instead, at criticality, the divergence of the correlation length leads to a non-trivial highly non-Gaussian function, which characterizes the transition. If hyperscaling holds, the PDF shows finite size scaling with the universal critical exponents of the model [14,15]. Moreover, all statistical moments, related to many-body correlation functions, can be extracted from this function. The PDF thus contains non-local information of the system, and it can be connected to non-local order parameters in certain quantum phases.
The PDF of magnetization has been exhaustively studied in classical spin systems undergoing phase transitions [16][17][18][19][20][21][22][23][24][25], and for those models, non-Gaussian distributions exhibiting finite size scaling at criticality have been found. However, much less attention has been paid to its quantum counterpart, and only few models such as the transverse Ising model have been previoulsy investigated [26,27]. Non-Gaussian distributions of light fluctuations have been observed for cold atomic ensembles [28], but nevertheless, a systematic treatment in the case of strongly-correlated atoms is currently missing. In this context, quantum coherent fluctuations of the individual constituents can lead to non-Gaussian PDFs, which are a necessary resource for continuous variables quantum i nformation processing [29].
The quantum-classical correspondence dictates for the quantum distribution the same finite size scaling with the corresponding critical exponents. However, to what extent the exact shape of the PDF is universal or dependent on the microscopic details represents an intriguing question, which we address in this work. Here, using exact solutions and the density matrix renormalisation group (DMRG) [30,31], we obtain the PDF of collective spin variables for different types of quantum phase transitions (first and second order, and BKT type), paying special attention to their behavior at criticality, and compare these results with their classical counterparts. Furthermore, we give an example on how the PDF can be connected to a non-local order parameter.
In particular, we investigate realizations using ultra-cold atoms. Although at present relevant energy scales are still too small compared to the lowest temperatures achieved in current experimental setups, current efforts focused on finding efficient protocols to lower the total entropy of the system [32], for instance, by using Raman transitions [33] or rehaping the trapping potential [34], may well make such investigations feasible in the near future. Here, we propose two possible experimental setups based on optical lattices for the measurement of the PDF, the first employing high resolution microscopy and the second using quantum polarization spectroscopy.

Probability density functions and finite size scaling
The probability to observe an eigenvalue m of an operator M is given by P (m) = µ m µ | ρ |m µ , where ρ is the density matrix describing the system and {|m µ } is an orthonormal set of eigenstates of M compatible with m. If M is the order parameter, this distribution can be related to the free energy functional F[m] appearing in Landau formalism [2]: P (m) ∼ e −F [m] . F[m] can be approximated by a power series of m and, if the correlation length does not diverge, the lowest order terms dominate. To lowest order, this leads to a PDF which is approximately Gaussian, a result which can also be understood by the central limit theorem. In contrast, close to the critical point we expect a highly non-Gaussian PDF.
Close to a continuous (second order) phase transition induced by a parameter g of the Hamiltonian with critical point g c , the correlation length diverges as ξ ∝ τ −ν , where τ = |g − g c | and ν is the critical exponent. Close enough to the critical point, the finite size scaling hypothesis implies that the mean value M scales with the system size L as [35]: where f is an analytic function. It is often more convenient and accurate for determining the critical point, both numerically and experimentally, to compute the Binder cumulant: as its scaling depends only on ν and not on β. That is, U =f (τ L 1/ν ), with a different scaling functionf . U quantifies the Gaussianity of the PDF, being null for a Gaussian distribution centered at zero.
Here we go beyond the first moments and consider the full probability distribution function P L (m, g). The renormalised PDF is expected to be a universal function (although still depending on the boundary conditions): withm = m/L −β/ν and r = L/ξ. Note that in this expression the parameter g driving the second order phase transition needs to be changed when varying L, in order to keep ξ/L fixed. This finite-size scaling of the PDF implies that hyperscaling and finite-size scaling of higher order correlation functions hold [16,17,21]. In fact, and as we will see, when the critical exponents are unknown, or not defined in the standard way as in the BKT transition, one should instead rescale the PDF and the eigenvalue m with σ = M 2 , directly computed from the PDF: withm = m/σ.

Models and Methods
We consider a variety of spin-1/2 chain models encapsulated by the Hamiltonian: where the sum on i extends over the L spins in the chain. The Pauli operators of spin i are denoted by σ i α , while J α and h α are coupling constants and magnetic fields along different directions (α = x, y, z). This model exhibits various transitions [1,36,37], some of which will be discussed below. We will study the behaviour of collective operators such as the total magnetization M α = i σ i α /L and its staggered counterpart Analytical results for the PDF can be obtained for models that can be written with a free fermionic representation after the Jordan-Wigner transformation, i.e. J z = h x,y = 0 in Eq.(5) assuming periodic boundary conditions, and operators that can be written as the sum of single-site operators in this basis, as for example M z . The order parameter M x however, is not separable and contains the string order operator. In principle all its powers could be computed by means of Wick's theorem, but the evaluation becomes involved as the order increases. At criticality and for the quantum Ising spin chain, the PDF can be exactly obtained by exploiting the relation with Kondo physics [27], but in general, reconstructing the PDF analytically is a very challenging task. Here instead, to obtain the PDF, we combine two different numerical methods. We use exact diagonalisation for sizes up to 20 sites, and also the time-dependent density matrix renormalisation group (t-DMRG) [30] with open boundary conditions. The t-DMRG method provides an efficient way to evaluate the characteristic function, defined as the Fourier transform of the PDF. For a pure quantum state |Ψ 0 this can be written as: being M the operator of interest and {m} the corresponding set of eigenvalues. The expression for χ(u) is equivalent to the overlap between the initial state |ψ 0 and its evolution under a fictitious Hamiltonian equal to M at different times u, and thus, it can be computed with t-DMRG. Since there is a one-to-one correspondence between the characteristic function and the PDF, this can be recovered by inverse Fourier transforming χ(u). Moreover, it allows to directly evaluate all central momenta and cumulants.

Second order quantum phase transition
We set J y,z = h x,y = 0, J x = J < 0 and h z = Jg in Eq. (5), corresponding to the ferromagnetic transverse Ising model, which exhibits a second order phase transition at g = ±1 separating a ferromagnetic ordered phase (FM) at low fields and a paramagnetic disordered phase (PM) at high fields. In Fig. 1(a) we show the numerical results for the PDF of the spontaneous magnetization M x , the order parameter in FM. Remarkably, we have obtained data collapse already for relatively small system sizes, assuming the ansatz of Eq.(4) with the predicted values β = 1/8 and ν = 1 (see e.g. [14,15]). This result reinforces the scaling hypothesis of all the statistical moments of the order parameter. Away from the critical point and in the thermodynamic limit (r → ∞) the PDF is Gaussian, in agreement with Landau theory and the central limit theorem.
In the disordered phase, it is centered around m = 0, whereas in the ordered phase the distribution is bimodal with the two peaks corresponding approximately to the broken-symmetry values of the order parameter. In contrast, at criticality (r → 0) the distribution becomes non-Gaussian. The central limit theorem does not necessarily hold anymore due to quasi long-range fluctuations in this regime. We note that, despite the PDF is universal at criticality, it still depends on the boundary conditions that are chosen. For comparison, it is worth noting that for the transverse magnetization M z , which is not the order parameter, the PDF is always Gaussian. For this observable, the moment of n-th order can be directly evaluated from the characteristic function and it always scales as L 1−n . This leads to a Gaussian distribution in the thermodynamic limit.
As previously discussed, it is not clear a priori that the PDF is a universal function (only depending on the universality class), and thus, that it can be retrieved from the classical model analogue. The 1D transverse Ising model maps onto the classical 2D Ising model with spatially anisotropic couplings. The anisotropy depends on the temperature of the quantum system, and therefore, on the system size of the equivalent classical model. Performing finite size scaling in the 2D model simultaneously in all dimensions (as it is usually done for a z = 1 transition) changes the anisotropy at each step, and one might doubt whether this process leads to a universal function. To elucidate this, we compare the previous results with those obtained for the classical model by Monte Carlo in [20,38]. A direct comparison shows that the quantum results for the PDF at the critical point seem to have more fluctuations around m = 0, leading to a flatter distribution with non-zero value, in contrast to the classical case. Gaussianity of the PDF.-We can characterize the gaussianity of the distribution by fitting the data with the sum of two Gaussians as with σ 2 and m 0 as free parameters (see solid lines in Fig. 1(a)). As shown in the figure, we find very good agreement away from the critical point, whereas close to the critical point the fitting yields poor results and one should use instead the exponential of a higher-order polynomial. This can be quantified by the non-linear least square fitting correlation coefficient defined as: where δy 2 is the squared norm of the residuals of the data y and Var[y] the corresponding variance. This quantity is close to one when the fitting works well, whereas it drops to smaller values for a poor fitting. Fig. 2(a) shows c as a function ofr = r sgn(g − 1) = L(g − 1). Clearly, there is a sudden decrease close to the critical point (r = 0), whereas it tends to one away from the critical point (r → ±∞).
An alternative way to quantify the Gaussianity of the PDF is by evaluating the Binder cumulant U , Eq.(2). U vanishes for a Gaussian distribution centered at zero, whereas it tends to U → 2/3 for a distribution close to two symmetric delta functions. Evaluating U is convenient for locating the critical point, since its finite size scaling only depends on the critical exponent ν and not on β, that is, U =f (r), and in this case U =f (|g − g c | · L). Thus, when plotted as a function of the transverse field g, the crossing point for data corresponding to different sizes tends to the actual critical point g c . When instead, plotted as a function ofr, they collapse to the same curve, which depends, however, on the boundary conditions.
In the same figure we also plot the result of U for the transverse Ising model with open boundary conditions (OBC) as a function ofr (b) and g (c) for different system sizes, ranging from L = 20 to L = 100. The result for periodic boundary conditions (PBC) is also shown for comparison.

First order phase transition
We consider the same Hamiltonian as before but with an additional longitudinal field h x = Jh = 0. At fixed value of the transverse field g = ±1 and approaching the Ising critical point by varying only h, the quantum phase transition is characterized by a different set of exponents (β = 1/15 and ν = 8/15 [14]). Our results (not shown) support strong evidence of scaling of the PDF.
If instead, the longitudinal field h is varied across zero, but at fixed value of the transverse field |g| < 1, the system undergoes a first order transition between two ferromagnetic states with opposite magnetization M x . At this transition, neither the correlation length diverges nor the gap closes. Nevertheless, inspired by Ref. [39], we propose a finite size scaling of the PDF.
In a finite size chain, two different energy gaps are present in this model. One is the real energy gap in the thermodynamic limit, which at h = 0 is given by ∆E = ||g| − 1| νz with z = 1, and it closes at the Ising critical point (g = ±1). The other gap δ separates the two lowest eigenstates and it is minimum at h = 0, with value δ 0 = 2|J|(1 − g 2 )g L , decreasing exponentially with L (|g| < 1). In [39], following dimensional arguments, the authors assume that δ and M x only depend on the longitudinal field h and L through the ratio q between the two energy scales: the first associated with the longitudinal field hLM x,0 , being M x,0 = (1 − g 2 ) 1/8 the magnetization at h = 0, and the second δ 0 : Thus, δ ∼ δ 0 (g)f 1 (q) and M x ∼ m 0 (g)f 2 (q), where f 1,2 (q) are analytic functions at q = 0. Here, we conjecture a similar dependence for the correlation length: ξ ∼ ξ 0 (g)f (q) where f is analytic at q = 0, and ξ 0 (g) = ξ(q = 0, g) diverges close to the Ising critical point with the usual power law: ξ 0 (g) ∼ ||g| − 1| −ν . Moreover, any observable M L (g, h) depends only on the ratios r = L/ξ 0 and q when rescaled with the proper critical exponents, e.g. for the magnetizationM x (r, q) = L −β/ν M x,L (g, h). Indeed, at fixed values of r and q (and for large enough chains, or equivalently weak enough external fields), we observe again data collapse for the PDF for different lengths, as shown in Fig. 1(b). We find the distribution to be always bimodal with the relative height of the two peaks ruled by h. We fitted the data with the sum of two Gaussians and found reasonable agreement except at the critical point, due to non-linear effects in the Landau potential.

BKT transition
For completeness, we finish this analysis with a different type of transition, the Berezinski-Kosterlitz-Thouless (BKT) transition in the spin-1/2 XXZ model. We set h x,y,z = 0, J x,y = J > 0 and J z = J∆ in Hamiltonian (5). The phase diagram is well known: the ground state is ferromagnetic for ∆ < −1, critical for |∆| < 1, and with Néel order for ∆ > 1. The BKT transition is at ∆ = 1. This transition is of the same universality class as the classical 2D XY model, for which previous calculations on the PDF [23] and on the Binder cumulant [24,25], show that the order parameter distribution is non-Gaussian at criticality. Here we compute the PDF for the ground state of the XXZ model, for both staggered magnetizations M st x and M st z , corresponding to the order parameters in the critical and Néel phases respectively, and for the full range of values of ∆. In contrast to the previous sections, we rescale the quantities m and P (m) with σ as defined in Eq. (4). For M st z we also fix r = L/ξ, being ξ = e π/ √ |∆−1| [37], and observe data collapse for different system sizes, as shown in Fig. 3(a). As expected, the distribution tends in the thermodynamic limit to a double-and single-peaked Gaussian for ∆ > 1 (ordered phase with respect to M st z ) and ∆ < 1, respectively. As shown in Fig. 3(b), the situation is much more interesting for the PDF of M st x , for which the spin-spin correlations do not decay exponentially in the critical phase (|∆| < 1). Outside this interval the distribution is Gaussian and centered at zero, because the operator M st x is disordered, while in the critical phase it is highly non-Gaussian. We observe again data collapse for different chain lengths, but now when fixing the value of ∆. The scaling with fixed ∆ is expected to occur in a critical phase, where the energy gap has already closed. When crossing the first order transition at ∆ = −1, the PDF shows also a discontinuity, with a sudden jump from a Gaussian to a function that has a singular behavior. In contrast to the usual behavior of the PDF, this function does not tend to zero as m/σ → ±∞, but it increases its value and its derivative when increasing the system size (see blue symbols in Figure 3(b)). Across the critical phase, the function changes continuously from a double-peak structure for ∆ < 0 to a single peak distribution for ∆ > 0. At the BKT point the PDF tends again to a Gaussian in the thermodynamic limit. Note that the results for the ferromagnetic model J < 0 are equivalent when analyzing M x instead, and changing ∆ ↔ −∆.
We emphasise that while the collapse is expected with the correct universal exponents as in the classical case, the specific collapse function can be dependent on the model or on the boundary conditions. This is clear by visual inspection when comparing the previous results at different values of ∆ to those in the 2D classical XY model reported in [23] (Figure 2).  The non-Gaussianity in the critical phase is also evident from the analysis of the Binder cumulant. In Fig. 4

Non-local order parameters
The BKT transition is particularly relevant in the context of 1D optical lattice gases because it rules the superfluid (SF) to Mott insulator (MI) transition in the Bose-Hubbard model. Close to the MI, where density fluctuations are small, and for large enough integer fillings, the model can be approximated to an effective spin-1 model [40,41]. The SF and MI are mapped respectively to the critical ferromagnetic phase, and to a state with local magnetization S i z = 0 perturbed with tightly bound particlehole fluctuations. The nature of density fluctuations is however different in the two phases. In the MI the non-local correlations can be characterized by the parity operator defined as where δn j is the local excess density from the average filling [42][43][44], and corresponds to the local magnetization in the magnetic model. O 2 P is finite in the MI, while it vanishes in the SF, and it has been experimentally reconstructed in [44] using single site microscopy.
In this context, an additional motivation for the study of PDFs is that it captures these non-local correlations. Indeed, it is easy to see that the characteristic function X(u) = m e ium P (m) for the suitable collective variable (in this case, the average magnetization of a subblock of length l) at frequency u = lπ is equal to the parity operator lim l→∞ X(u = lπ) = O 2 P . Thus, a measurement of a collective variable and the reconstruction of its PDF provides an alternative route to evaluate some particular non-local order parameters.

Experimental reconstruction of the PDF
The PDF for the local number parity operator can be reconstructed using single site resolution microscopy [45,46]. Recently, novel schemes which allow to circumvent light-assisted pair loss and resolve the internal atomic degree of freedom have been also demonstrated [47]. An alternative proposal is based on the Faraday effect, and consists in analysing the polarization fluctuations of a strongly polarised laser beam that interacts with the atomic sample [48][49][50]. Since this is based on dispersive lightmatter interaction, it has the advantage of being less destructive, which could potentially pave the way for multiple probing of the same ensemble. Thus, this method could find relevant applications for quantum state control and preparation. In particular, by monitoring the system in the appropriate way, this could be used in principle to lower its total entropy. Already it has been experimentally demonstrated that QPS combined with real-time feedback allows for a dramatic reduction of the atom number fluctuations [51], increasing the purtity of a coherent spin state [52], and inducing spin squeezing [28]. Moreover, QPS can be exploited for quantum state engineering and control [53], inducing non-trivial quantum dynamics [54], and the conversion of atomic correlations and entanglement into the light degree of freedom for quantum information processing [55].
The light-matter interaction with the collective spin operator M , H eff ∼ (κ/τ )LP ph M , is written in terms of the momentum-like light quadrature P ph measuring the photon fluctuations in the circular basis with respect to the strong polarisation axis, and τ is the interaction time. P ph is canonically conjugated to its position-like counterpart: [X ph , P ph ] = i . The coupling constant κ = (ηα) 1/2 can be expressed in terms of the single atom excitation probability or destructivity η, and the resonant optical depth α = N σ cross /A, where A is the overlap area between light and atoms, σ cross is the resonant cross section and for a one-dimensional system N = L/d, being d the interparticle distance. By adjusting the laser intensity, the destructivity parameter η is typically set to values smaller than 0.1 to limit the fraction of excited atoms. The resonant optical depth α should be maximized in order to obtain the largest possible coupling between light and atoms. For ultracold atoms in a one-dimensional optical lattice, considering L = 100 and A = 0.5µm 2 we obtain α ≈ 8 and κ ≈ 1. We definẽ κ = κ/ √ L, which will be approximately independent of L and could be (in the best case scenario) as large asκ ∼ 0.1. We will show that this value of κ is large enough to reconstruct the PDF of the spin operator M by looking at the distribution of the light quadrature X ph . Larger values ofκ could be engineered by coupling atoms with optical cavities [56,57], nanophotonic crystals [58,59] or optical nanofibers [60].
It is possible to show [61] that the light distribution is the sum of vacuum Gaussian distributions of light each displaced by a quantity proportional to the eigenvalue m of the operator M and scaled by the probability P (m) to observe such eigenvalue: where we have considered a Gaussian input light beam with variance σ 2 ph , being 1/2 for the vacuum state.
To evaluate the effectiveness of the method we compare the actual atomic spin distribution with the one of the output light. The distance between both distributions decreases exponentially withκ/σ ph . Thus, one could in principle improve the fidelity by increasingκ or using squeezed light [62]. We show in Fig. 5 the result for the transverse field Ising model for an optimal caseκ = 1 and σ 2 ph = 1/2 and for a more realistic value ofκ = 0.15, but squeezed input light with σ 2 ph = 1/4, recently achieved [62]. For the former, the light distribution faithfully follows the magnetization, whereas for the later, it only agrees qualitatively, but still captures the peaks. Experimentally, one would need to repeat the measurement a number of shots N shots ∼ L 2 to estimate the PDF.

Summary
In conclusion, we have shown that the distribution of collective variables in spin models reveals relevant information of quantum phase transitions. We have shown that for a range of quantum phase transitions a non-Gaussian distribution of the order parameter is a clear signature of criticality, and that the scaling hypothesis holds. Finally we have proposed an experimental method for its measurement using light-matter interfaces, and discussed its feasibility for realistic values.