Higher Order Spin Noise Statistics

The optical spin noise spectroscopy (SNS) is a minimally invasive route towards obtaining dynamical information about electrons and atomic gases by measuring mesoscopic time-dependent spin fluctuations. Recent improvements of the sensitivity of SNS should make it possible to observe higher order spin correlators at thermodynamic equilibrium. We develop theoretical methods to explore higher order (3rd and 4th) cumulants of the spin noise in the frequency domain. We make predictions for the possible functional form of these correlators in single quantum dot experiments and then apply the method of the stochastic path integral to estimate effects of many-body interactions.


INTRODUCTION
The full counting statistics has become a widely discussed topic both in electronics and quantum optics [1][2][3][4][5][6]. Its measurements promise to provide considerably more information about interacting electrons and photons than that could possibly be obtained from standard linear response characteristics.
Spins in a small mesoscopic volume of a semiconductor with N electrons typically experience statistical fluctuations of order √ N , even in a zero magnetic field at the thermodynamic equilibrium. The spin noise spectroscopy (SNS) is an optical technique that provides a viable route to study such fluctuations by directly measuring local spin correlators [7,8]. Experimentally, spin correlators were obtained by measuring the rotation of the linear polarization of a light beam that passes through a mesoscopic region with spins, e.g. as shown in Fig. 1. The polarization rotation angle is proportional to the local instantaneous magnetization of the region and can be traced with sub-nanosecond resolution in time [9]. Characterization of the equilibrium spin noise has been demonstrated and used to determine g-factor, spin coherence and spin relaxation times of electrons in GaAs [10][11][12][13][14][15][16] and atomic gases [8,9,[17][18][19][20]. An application of SNS to the central spin dynamics in InGaAs hole doped quantum dots [21][22][23][24], in particular, revealed the importance of the nuclear quadrupolar coupling for the decoherence and relaxation of a spin qubit [25].
Higher order spin correlators are not the subject of application of the standard fluctuation-dissipation relations. Hence, by studying higher order statistics, one can obtain information about spin dynamics that simply cannot be found in the average linear response characteristics of spin systems [26][27][28]. For example the so called phase transitions at fluctuating level attract lots of interest [29]. Recently, it was realized that such unusual critical phenomena can be effectively studied in physical systems by measuring time-dependent behavior of high order noise cumulants [28,30].
Most of the studies, however, have been focused so far on the spin noise power spectrum [25,[31][32][33][34] and the associated second order spin correlator where S z (t) is the time-dependent spin polarization in a mesoscopic region and z is the measurement axis. The discussion of higher order spin correlators at the thermodynamic equilibrium has been essentially absent. This is partly due to the fact that for a large number N of spins the physical noise is dominated by Gaussian fluctuations which are fully described by (1), as articulated by the law of large numbers. Hence, in order to reveal the additional information about spin dynamics in the form of higher order cumulants one should not only filter the useful signal from the background noise but also filter out the physical Gaussian fluctuations. In noninteracting spin systems, non-Gaussian effects are due to discreteness of spin states. Among the most relevant theoretical and experimental work, which refered to the optical spin noise spectroscopy, we mention here the model of a weak measurement framework for a single spin system, developed in [35], which discusses time interval statistics between detector clicks. The type of time correlators studied in that work, however, is somewhat distant from the one that is currently accessible by SNS. On the experimental side, in [36] and [37], simple examples of measurements of higher order noise statistics of acoustic sound in the frequency domain were provided, claiming that similar methodology may work to study the spin noise in an SNS setup. SNS appears to be particularly promising for characterization of the high order spin statistics. It provides a considerable and continuous stream of data for statistical analysis. Information can arrive with the rate of gigabites per second and the measurement time is formally unrestricted, e.g. it can be several weeks if needed. Moreover, the sensitivity of SNS has been continuously improving. For example, recently introduced Ultrahigh Bandwidth SNS achieved picosecond time-scale resolution [19]. New measurement schemes have also been proposed recently to increase the polarimetric sensitivity by using high polarization exctinction schemes and hence suppress the relative role of the background noise sources [38]. Noise of a single spin of a heavy hole localized in a flat (InGa)As quantum dot has been successfully observed recently by SNS [39]. The latter achievement is an important milestone on the way to experimental studies of the higher order spin statistics because the law of large numbers no longer applies to N = 1, which means that a number of higher order cumulants of spin noise statistics should be of the same order of magnitude as the spin-spin correlator (1). Non-optical methods have also been developed and demonstrated recently to probe correlators of single spins [40,41], which may be considered for studies of higher order spin noise statistics.
It is also helpful that optical SNS is free from many of the problems that were encountered in measurements of higher order statistics of electric currents [42][43][44]. In SNS, since spin fluctuations are probed optically, there is no problem with the often detrimental noise from leads, and one generally does not need a complex nano-lithography. The advancement of the spin noise spectroscopy thus motivates us to reconsider the modern technological capabilities of using higher order noise statistics measurements as an important tool for materials studies and quantum information science.
The plan of our article is as follows. In section II, we define the 3rd and 4th order spin noise cumulants and discuss their properties. We will then consider four different applications (Fig. 2) in which we believe the experimental observation of higher order spin correlations is most likely: In section III, motivated by the observation of the spin noise of a single central spin in a quantum dot, we explore a possible form of 3rd and 4th order correlators of a single spin ( Fig. 2(a)) and make specific predictions for verification in InGaAs quantum dots [39]. In section IV, based on a simple non-interacting spin model, we develop a method based on the stochastic path integral to calculate spin noise statistics and, for illustration, we apply it to an ensemble of a mesoscopic number of noninteracting spins ( Fig. 2(b)). Then in section V we explore the effect of many-body interactions, such as the ones arising near a ferromagnetic phase transition in magnetic semiconductors (Fig. 2(c)). We apply the path integral method to a phenomenological kinetic model based on the Glauber spin dynamics to demonstrate that higher order spin noise statistics becomes particularly insightful to observe near a phase transition point. In section VI, we explore the higher order spin cumulants for conducting electrons (Fig. 2(d)) and explore the effect of the Pauli principle.

II. PROPERTIES OF 3RD AND 4TH ORDER SPIN CORRELATORS IN FREQUENCY DOMAIN
Complete information about an interacting many-body system is contained in the full set of all cumulants of system variables. Let us introduce a normalized spin polarization and consider a random variable where T m is the time of measurement, which is assumed to be much larger than other physical time-scales in the system, including the characteristic relaxation time T of the spin system. Note that since S z (t) is real, we have a * (ω) = a(−ω). For a paramagnetic system, one finds that a(ω) = 0, where the average is considered over repeated measurements under identical experimental conditions. The most accessible physically interesting characteristic is the noise power, which is defined as Its knowledge is equivalent to the knowledge of a spin-spin correlator in the time domain via At steady conditions, only products of a(ω i ) with i ω i = 0 can be nonzero after averaging. Hence, the next nearest nontrivial correlator of a(ω) is of the 3rd order and given by A specific property of C 3 is that it is zero in a system with a time-reversal symmetry. There can be a variety of potentially interesting combinations of 4th powers of a(ω). The one that has the most transparent physical meaning is the so called bispectrum: For ω 1 = ω 2 its definition is modified: The bispectrum (7) indicates how spin noise components at two different frequencies "talk" to each other. For example, if C 4 (ω 1 , ω 2 ) is negative, one can conclude the presence of anti-correlations, i.e. observation of a strong signal at one frequency presumes that another frequency is likely suppressed, etc. The form of, so called, cumulants (6), (7) and (8) is dictated by the requirement that they are zero for Gaussian fluctuations of a(ω), so that they do not duplicate the information that can be obtained from (4). One more important property of cumulants is their "additivity", i.e. independent noise sources additively contribute to their values. This is important because experimental measurements in SNS generally produce a strong background shot noise, which can be filtered out by separately measuring the pure background in a strong magnetic field applied in the direction transverse to the measurement axis. After subtracting the background contribution from the measured cumulants one obtains the physical values of cumulants.
Another consequence of such an additivity of cumulants is that they generally linearly increase with the measurement time T m . Indeed, if T m is much larger than the relaxation time T , one can expect that a measurement during T m is roughly equivalent to T m /T independent measurements which additively contribute to the finite cumulant value. This property is one of the reasons of the difficulty in observing higher order statistics. For example, the dimensionless ratio C 4 (ω)/(C 2 (ω)) 2 should generally be proportional to T /T m ≪ 1. Experimentally, it is important to keep T m large to avoid nonphysical effects of a finite measurement time. On the other hand, one should keep the ratio T /T m not too small in order to be able to filter C 4 (ω) from the Gaussian part of spin noise statistics.
As another word of caution, we would like here to point to a problem that has not appeared in previous measurements of the noise power spectrum (4). Most cumulants are protected only against an additive background noise, i.e. the noise that is uncorrelated from the physical signal. In SNS, fluctuations of spins are deduced from the fluctuations of an optical signal. The probability of light interaction with spins depends on the laser beam intensity, which is often slowly fluctuating. For C 2 and C 3 such fluctuations merely renormalize their amplitude, keeping relative amplitudes at different frequencies the same. In contrast, the expression, e.g. (7), for C 4 is the difference of two terms that may be differently renormalized by the beam intensity fluctuations. Moreover, for a large number of noninteracting spins in the observation region, each of those terms can be much larger than their difference in (7) so that even small instabilities in the beam intensity can lead to an admixture of Gaussian fluctuations and degrade the measurements of C 4 . Hence C 4 should become most accessible when the number of observed spins is relatively small so that the dimensionless ratio, C 4 /(C 2 ) 2 , is as large as possible and the cumulants are averaged over a smaller span of time than the time scale of slow fluctuations of the beam intensity. This situation can be realized in a single spin quantum dot as in [39] or in the case of considerable spin correlations, e.g. when one effectively observes dynamics of a small number of ferromagnetic domains.
Higher order correlators depend on more than one frequency, so obviously, they contain additional and, possibly, considerably larger amount of information about the system than the noise power. In solid state applications, as a proof of principle, C 3 (ω 1 , ω 2 ) was measured as a function of frequencies only recently to describe the shot noise of electric charge currents through an artificial nanostructure [36]. Specially engineered nanoscale systems could allow measurements of charge current cumulants in the time domain up to 15th order and the study of fundamental nonequilibrium fluctuation relations [26][27][28] but such measurements could not be applied to materials characterization. Therefore, we will explore the information that can be obtained about condensed matter systems by measuring higher order spin, rather than current, correlators.

III. SINGLE SPIN NOISE
A. 3rd order correlator of Ising spin dynamics From the point of view of statistical filtering of a useful signal, the 3rd order correlator (Eq. 6) is the next in complexity after the popular 2nd order spin-spin correlator. However, the 3rd order correlator changes sign under time reversal and hence its observation requires a specific breakdown of a time reversal symmetry. In a single spin InGaAs quantum dot, conditions for 3rd order cumulant observation can be created by applying a strong (of the order of 1 Tesla) magnetic field in the out-of-plane (i.e. parallel to the measurement z-axis) direction. At such fields, spin relaxation is dominated by interactions with phonons. At a sufficiently low temperature, Zeeman coupling in such a magnetic field becomes comparable to k B T s , where T s is the system temperature and k B is the Boltzmann constant. Without an in-plane magnetic field component, one can disregard coherence effects and assume that the spin essentially behaves as a classical Ising spin. Spin transitions between up and down states can be then described by kinetic rates, respectively, k 1 and k 2 , which satisfy the detailed balance condition where B z is the z-component of the magnetic field and g z is the corresponding g-factor. Spin polarization dynamics of this system behaves like a telegraph noise, as shown in Fig. 3(a). Let p 1 (t) and p 2 (t) be time-dependent probabilities of a spin being, respectively, in the up and down states. The dynamics of probabilities is governed by the master equation:ṗ which can be solved as P(t) = U(t)P(0), with P(t) = (p 1 , p 2 ) T , where the evolution operator U is given by It can be shown that the average spin polarization is a constant: which is equal to the average spin at t = 0: S z ≡ 1|σ z |P 0 = k2−k1 k1+k2 , where σ z is the Pauli matrix. Here 1| = (1, 1) and the initial spin polarization |P 0 = (p 1 , p 2 ) T with p 1,2 = k 2,1 /(k 1 + k 2 ). Note that, for simplicity, here we normalized the Ising spin polarization values to be ±1 rather than ±1/2.
The second order correlator of spin fluctuations in real time is then given by The spin noise power spectrum, defined as in Eq. (5), reads: with the relaxation rate γ = k 1 + k 2 . The 3rd order spin correlator can be calculated similarly. For t 1 > t 2 > t 3 , we have: For the case of t 2 > t 1 > t 3 , the expression for C 3 can be obtained by exchanging t 1 and t 2 . Finally, by taking the Fourier transform, the 3rd order cumulant in the frequency domain is found: One can conclude from (16) that C 3 is nonzero only when k 1 = k 2 , i.e. when a strong magnetic field changes the relative rates of spin transitions. If k 1 − k 2 changes sign, so does C 3 . Measurements of C 3 in a quantum dot system can be used, e.g. as an independent probe of the g-factor via (9) at high values of the magnetic field. In Fig. 3 Finally, we note that noise in simple kinetic models have been previously studied in different contexts, e.g. in statistics of electric currents coupled with a fluctuating two level degree of freedom [26,42,45]. For example, spin correlator (16) has the same functional form as the third order electric current correlator derived in [45].

B. 4th order cumulant of the Ising spin
The advantage of the 4th order cumulant is that it is nonzero even without time symmetry breaking. For a single spin, it is expected to be nonzero because the discreteness of spin states leads to a binary signal, which is similar to a telegraph noise rather than Gaussian fluctuations. Since strong out-of-plane magnetic fields are not needed for its observation, one can measure C 4 (ω) in a standard setting, e.g. in a zero magnetic field. The spin of a heavy hole in an InGaAs quantum dot at such conditions behaves essentially as an Ising spin with Markovian transitions between up and down states due to the coupling to a quickly fluctuating nuclear spin bath [25].
While, in principle, analytical calculations of C 4 via finding its value in the time domain are possible, we found them very tedious even for the simplest models. In the next section, we will describe an alternative approach, which is based on the stochastic path integral technique, and which provides a simpler framework for the calculation of C 4 , including for independent Ising spins. In this subsection, instead, we will present a numerical approach that simulates the weak measurement setup. Such simulations are close in spirit to the real physics encountered in the experimental measurement process of SNS and allow us, in particular, to study the effect of detector parameters on the observable higher cumulants.
As in the previous subsection, we consider an Ising spin but assume that spin flips between up and down states happen with the same rates k = 1/T , where T is the characteristic life-time that for a hole doped InGaAs quantum dot was estimated to be ∼ 0.4µs [22]. We simulate the weak measurement process with the scheme that was proposed for SNS in [35]. According to it, measurements are performed in discrete time steps τ ≪ T . Suppose that with a probability p D < 1 the coupling to the detector induces the collapse of the state vector at each time step, leading to a detector "click", i.e. the measurement of the spin state. Respectively, with a probability we observe no useful signal at the detector per single measurement. In (17), we introduced the time T D that characterizes the typical time between collapses of the spin vector to the measurement basis.
Our measurement axis is the z-axis, and there are three possible measurement outcomes. Let us call them 1, 0, −1, where ±1 correspond to the collapse of the state vector to one of the spin states |± , which also correspond to specific ±1 outputs of the detector, while 0 corresponds to the absence of a useful detector click.
If the probability per single measurement step to collapse the spin state is small, most of the measurement outcomes would be zeros. However, when the collapse of the spin vector happens, we determine the state of the spin, and the measurement produces +1 or −1. After this, the evolution of the spin is again calculated according to the master equation. A numerical simulation of such a process returns a random sequence that simulates the output of a detector, such as 0000 1 000000000000000 −1000000000000 1 000 . . ., etc.
Evaluating it for a span of time T m ∼ 256τ , we found its Fourier transform and calculated the powers of a(ω). By repeating this process, we produced a 100 million of such random sequences, from which we calculated the cumulants. As in SNS experiments, the fact that most of the time a weak measurement fails to produce the detector response leads to a broad background white noise, which we subtracted to obtain the physical part of the spectrum. Our results for the noise power, Eq. (4), and for the 4th order cumulant, Eq. (8), are shown in Fig. 4. As expected, the noise power spectrum (Fig. 4(a)) is Lorentzian. The 4th correlator appears to be negative and narrower than C 2 . The ratio of the integrated cumulants is found to be η 1 = −1.9(T /T m ), which is close to the theoretically predicted value η 1 = −2(T /T m ) (section IV) in the limit T m ≫ T ≫ τ . The numerical simulations of measurements of C 4 , for a single spin, converge relatively quickly when the time T D is smaller than T , i.e. several nonzero detector clicks happen before a spin flip. In the opposite case, T D < T , convergence of the simulations to the expected values of C 4 (ω) quickly deteriorates, and the result appears to be dependent on the measurement time T m . Hence we predict that the favorable conditions for the reliable observation of C 4 (ω) should be found at strong beam intensities that increase the probability of light interaction with a spin of a quantum dot.

IV. STOCHASTIC PATH INTEGRAL APPROACH APPLIED TO NONINTERACTING SPIN SYSTEMS
Stochastic path integrals are frequently used for the calculation of physical noise statistics. Here, we will use the version of this technique that was developed by Pilgram et al. [46,47], and which is particularly suitable for obtaining cumulants of mesoscopic systems. Originally, this approach was built to study the full counting statistics of electron transport in mesoscopic electric circuits, and then adapted to applications in biochemical reaction networks and explicitly time-dependent systems [48,49]. In this section, we will consider the model of N independent Ising spins as an illustrative example, and in the following two sections we will use this approach to study the effects of ferromagnetic interactions and the Pauli exclusion principle on cumulants of the spin noise.
The basic idea of constructing the path integral is the separation of time scales. If the number of interacting spins N is mesoscopically large, one can identify the time interval ∆t such that the number of spins flipped during this time interval is, on average, much larger than unity but still much smaller than N . One can consider then the total spin polarization of the observed region as a slow variable, which can be approximated by a constant during time ∆t. Suppose that the probability k 1 is to flip from up to down and the probability k 2 is to flip from down to up states of any spin per unit of time. Since the spins do not interact with each other, this model, formally, can be solved without a path integral, as we did in Sec. III because cumulants for noise of independent spins simply add. Our goal, however, is to illustrate how path integrals can be used for calculations of higher order spin noise cumulants, in order to apply them later to more complex interacting spin systems.
Let N ↑ and N ↓ be the numbers of observed spins, respectively, in the up and down states. The average number of spins transferred from up to down during time ∆t is Q 1 = k 1 N ↑ ∆t. Similarly, the transferred number of spins from down to up is: Q 2 = k 2 N ↓ ∆t. Q 1 and Q 2 are random variables obeying the Poisson distribution. Their fluctuations result in the variation of the spin polarization, which we define as M = N ↑ − N ↓ . One can integrate over the fast fluctuating dynamics of Q 1 and Q 2 to obtain the partition function for the slowly changing variable M . During time ∆t, the variation of M is M (t + ∆t) − M (t) = 2(Q 2 − Q 1 ) and Q 1 and Q 2 can be considered uncorrelated to each other.
The Poisson probability distributions of Q 1,2 can be written as integrals P (Q 1,2 ) = dχ 1,2 e −iχ1,2Q1,2+∆tH1,2(χ1,2) , in which: are the cumulant generating functions of corresponding Poisson distributions. Following the standard path integral approach [46,47], we discretize a long measurement time, T m , into ∆t segments around time moments: t n = n∆t with n = 1, ..., N . The sum over all possible system trajectories weighted by their probabilities during this time T m , which is called the partition function, can then be written as: where the δ-function imposes the conservation constraint. One should then express the delta function as an integral over a new variable χ, i.e. as δ(f ) = (2π) −1 dχe iχf . Then one can perform integration over variables Q 1,2 , and χ 1,2 , to end up with : Taking a continuous limit, we obtain: At the steady state, Eq. (24) has a solution One can find that M C coincides with the mean value of the spin polarization in the system. In order to obtain higher order spin correlators, one should consider the action near the saddle point (25) beyond the classical limit and introduce finite values χ and δM that describe the deviation from (25). As it was discussed in [47], in order to calculate the n-th order correlator of variables, it is sufficient to find the Hamiltonian in the action up to the n-th order in total powers of χ and δM . Hence, we write the partition function as where To shorten our notation, we introduce parameters Keeping only the quadratic part of the Lagrangian, L 2 in (26), one can calculate the second order correlation function by writing the action in the frequency domain. By substituting χ(t) = T −1 m ω e iωt χ(ω) and δM (t) = T −1 m ω e iωt δM (ω) into the action, we find with a = (χ, δM ) T and the matrixÂ given by: Similarly, one can write the higher order contributions to the Lagrangian in the frequency domain, e.g., The 2nd order correlators are found from (we recall that Z = 1): Explicitly: Eq. (35) provides the spin noise power in the frequency domain. As expected, its value coincides with the noise power of the single Ising spin, Eq. (14), up to the factor N . Correlators of δM with χ do not describe a measurable characteristic but they are needed for the calculation of higher order cumulants.
To estimate the third order cumulant, we keep the term with L 3 in the action and treat it as a small perturbation: Since the exponent in (38) is quadratic in variables, one can calculate this expression using the Wick rule by summing over all possible products of 2nd order correlators. Terms in such expressions can be represented by the Feynman diagrams, e.g. one of them that contributes a nonzero value to the cumulant (38) is shown in Fig. 5(a). Other nonvanishing diagrams contributing to C 3 are obtained by all possible permutations of arrows entering the node such that the total sum of ingoing and outgoing frequencies is zero. Evaluating these diagrams, we obtain: The fourth order cumulant defined in Eq. (7), at k 1 = k 2 can be found as For the time-reversal symmetric case k 1 = k 2 ≡ k the result is represented by the sum of diagrams shown in Fig. 5 (b) and (c): Generally, there is also a tree diagram contribution to C 4 , which appears when the second power of L 3 is included in the perturbative expansion, as shown in Fig. 5(d). It is also proportional to the measurement time, T m , but it vanishes at identical rates of up-down and down-up transitions (k 1 = k 2 ), so we do not consider it here.
Finally, we discuss the relative strength of C 4 and C 2 . It is useful for this to introduce a dimensionless combination .

V. MODEL WITH A FERROMAGNETIC COUPLING
In order to explore many-body effects on spin noise statistics, here we will study a model with the Glauber dynamics of ferromagnetically coupled Ising spins [30]. We assume that all N observed spins experience an effective magnetic field proportional to the instantaneous spin polarization, B z (t) ∼ M (t), so that the kinetic rates are modified to account for this field and the detailed balance conditions (9). Here, for simplicity, we define M = , which is slightly different from the quantity defined in the noninteracting case.
In comparison to the noninteracting spin model of the previous section, we assume that k 1 and k 2 are no longer constant but rather depend on the local magnetization such that k 1 /k 2 = e −αM , where α is a parameter that characterizes the exchange coupling in the mean field approximation. This parameter also absorbs the inverse temperature α ∝ T −1 s in (9). We choose k 1 = ke −αM/2 and k 2 = ke αM/2 , where k is the characteristic kinetic rate at zero M . Applying the rules for constructing an effective Hamiltonian in the path integral, we find For such a choice of variables M and χ, the partition function is expressed in the form in which it is explicitly clear that the correlators, e.g. M (ω)M (−ω) , depend as 1/N on the total number of observed spins. By recalling that the full polarization is obtained by changing variables M → M N , one can conclude that all cumulants of the total spin polarization will be finally proportional to N . We will assume that parameter α is tuned so that the system is close to a ferromagnetic phase transition, so that the mean magnetization is either zero or small. In such an approximation, the saddle point equations for Eq. (44) have two solutions: one is χ = 0 and M = 0; the other is χ = 0 and M 2 ≈ 4 α 2 ( α−2 1−α/6 ). This means that there is a critical value of the parameter α, namely, α c = 2 that corresponds to the phase transition between paramagnetic and ferromagnetic phases.
A. Correlations in the paramagnetic phase (temperatures above the phase transition) At temperatures slightly above the phase transition point, we would have zero average magnetization, M = 0. We introduce a small parameter, t = (T s − T C )/T C , such that α ≈ 2(1 − t), and expand the Lagrangian up to 4th order in powers of fluctuations from this point: L = L 2 + L 3 + L 4 + . . .: where a = 2k, γ = 2kt.
Evidently, due to the time-reversal symmetry in the paramagnetic phase, the third order cumulant is zero. The second order cumulant has the same Lorentzian form as noninteracting spins in Eq. (35), (36), and (37), but with parameters a and γ defined in (47). At the phase transition, the effective relaxation rate vanishes γ → 0, indicating the critical slowdown. Hence, the amplitude of the 2nd order correlator grows as it approaches the phase transition, as shown in Fig. 6(a). Each term in L 4 corresponds to some contribution to the 4th cumulant of spin noise, so that 4 .
The first two terms in (48) are due to χ 4 and δM χ 3 in (46). We have previously calculated them for the independent Ising spins model. They produce the same expressions as Eq. (41) and Eq. (42) with redefined parameters according to (47). The terms ∼ δM 2 χ 2 and ∼ δM 3 χ in (46) are new. Their contributions to C 4 are given by In Fig. 6(b), we plot C 4 (ω 1 , ω 2 ) in the paramagnetic phase, according to which, the 4th order correlator is still negative but it can have a very large amplitude near the phase transition.
B. Spin noise in the ferromagnetic state (temperature below phase transition) A specific feature of the ferromagnetic state is the absence of the time-reversal symmetry, meaning that the 3rd order cumulant of the spin noise may become nonzero. Below the critical temperature, we would have the saddle point, which corresponds to M = ±M 0 = ± 3|t|. Expanding the Lagrangian up to the 3rd power of fluctuations near M = M 0 , we find In this case, a = 2k, γ = 4k|t|. Note that the relaxation rate is twice the one in the paramagnetic phase at the same distance to the phase transition point. Diagrammatic calculations show that the two 3rd order terms in (51) produce two contributions to the 3rd order cumulant: 3 , where Both terms depend linearly on the average magnetization M 0 , which in turn scales as M 0 ∼ |t| 1/2 as a function of temperature distance to the phase transition. The 4th order cumulant in this phase would be very complicated to show here properly. In Fig. 6(c) we just plot the numerical result for the dimensionless parameter η as a function of t for both t > 0 and t < 0. Since the parameter η characterizes the deviation of the distribution of the spin noise from a Gaussian form, its divergence at the phase transition point indicates a strongly non-Gaussian spin noise statistics. A difference between t > 0 and t < 0 could be traced to the difference of relaxation rates γ in the two cases.

VI. CONDUCTING ELECTRONS WITH PAULI EXCLUSION PRINCIPLE
As another example of nontrivial higher order correlations, consider the Fermi sea of conducting electrons. Previously charge current fluctuations in such systems have been studied extensively. Partial suppression of the shot noise by the Pauli principle has been one of the most important effects in this field. It was studied previously, in particular, by the method of the stochastic path integral [47]. Here we will explore effects of the Pauli exclusion principle on the cumulants of the local spin noise fluctuations. We assume that, due to the phonon coupling in the observation region, the local electron distribution in the momentum space, for each spin species, quickly equilibrates and restores to the Fermi-Dirac distribution at the ambient temperature T s and at local chemical potentials µ ↑ (t) and µ ↓ (t) for, respectively, spin up and spin down electrons. We assume that the spin degree of freedom equilibrates at a much longer time scale, e.g. the spin relaxation due to the Dyakonov-Perel mechanism is at ∼ 100ns, while thermalization of orbital degrees of freedom can be at sub-nanoseconds for conducting electrons in GaAs. Due to spin flipping, chemical potentials, µ ↑ (t) and µ ↓ (t) will fluctuate. We normalize the chemical potentials so that at a zero net spin polarization they are set to zero. The excess of electrons with spin up then can be related to the chemical potential by where D is the density of states in the observation region per spin and per unit of energy near the Fermi surface. Note that D is not an intensive characteristic in the sense that it is not a density per volume of the system. For example it is proportional to the size of the mesoscopic observation region in Fig. 1. We set the Boltzmann constant k B = 1. For a sufficiently large observation region, electroneutrality ensures that N ↑ = −N ↓ . To account for the Pauli exclusion principle, we assume that the average number of spin flips, e.g. from up to down, for electrons with energy ǫ is proportional to the number of electrons in the observation region with spin up and the density of unfilled states with the spin down at energy ǫ. Then the total average rate of transitions from up to down is given by: where k is the characteristic transition rate. Similarly, The local spin polarization is M = N ↑ − N ↓ . Since a single spin flip changes the spin polarization by 2, the Hamiltonian in the path integral action that describes the dynamics of M is written as: where χ is the variable conjugated to M . Explicitly: The saddle point solution corresponds to χ C = 0 and M C = 0. Expanding the Lagrangian up to the fourth order, we obtain: The second order cumulant would have the same Lorentzian shape as it would be for independent Ising spins with parameters a = 4kDT s , γ = 2k. The Pauli exclusion principle changes the expression of a by replacing the total number of spins N with the number of spins near the Fermi surface in the interval of energy of the order of temperature T s . The fourth order cumulant is the sum of three terms that correspond to contributions of each of the three terms in L 4 : where (63) The first two terms have the same form as for noninteracting Ising spins, (one can compare (63), (64) with Eq. (41) and Eq. (42)). The new contribution (65) appears from the term ∼ δM 2 χ 2 in (60). It radically changes properties of the bispectrum, as shown in Fig. 7. In particular, the dimensionless parameter η for the Fermi system becomes positive, η = 2 3a = 1 6kDTsTm . Comparing this result with noninteracting and ferromagnetic interaction cases, for which η is negative, we conclude that the functional form of the fourth order cumulant is sensitive to the details of spin interactions. Note also that at lower temperatures, the Pauli exclusion principle makes η larger, i.e. statistics becomes less Gaussian.

VII. DISCUSSION
Higher order cumulants, by construction, contain additional information to the noise power spectrum. For example, in all models considered in this work, the noise power has a Lorentzian form, from which only a single parameter, i.e. the effective relaxation time can be determined. In contrast, the 3rd order cumulant contains information about the asymmetry of relaxation rates, while measurements of the 4th order cumulant can be used to estimate not only all parameters of the considered models but also distinguish among candidate models if the Hamiltonian of the spin system is not a priory known, because the functional form and the sign of the 4th order cumulant are sensitive to subtle details of the kinetics. We also predict that ratios of cumulants, such as parameters η 1 and η in Eqs. (18) and (43) provide a good estimate of the size of physical spin correlations.
We explored the properties of higher order cumulants of the spin noise in the frequency domain and discussed the conditions for their experimental observation. In a standard framework of most of the SNS experiments, measurements were performed on a mesocopic number, e.g. N ∼ 10 5 , of independent spins. In such a case, the higher cumulants are suppressed in comparison to the noise power, e.g. C 4 /(C 2 ) 2 ∼ 1/N . Our results suggest two strategies that can be used to enhance this ratio. First, one can perform measurements on a smaller number of spins. Since the spin noise characterization of a single spin is now accessible in InGaAs hole doped quantum dots [39], we believe the measurements of higher order noise cumulants in such systems are already possible. Our analytical and numerical calculations predict a negative value of the 4th order cumulant at a zero external magnetic field. The 3rd order cumulant can be also observed in such systems in a strong out-of-plane magnetic field.
The second strategy to observe higher order spin cumulants is to perform measurements on strongly interacting spins. In the case of ferromagnetic interactions, spin fluctuations are strongly enhanced near the paramagnetic/ferromagnetic phase transition, leading to a much stronger signal for all cumulants. In such a case, spins flip not independently but rather as clusters of many correlated spins, and the number N should be interpreted as the typical number of clusters in the observation region, which can be considerably smaller than the total number of observed spins. Hence, we predict that the higher order cumulants can be substantially more important for the characterization of magnetic semiconductors, especially near the phase transition temperatures.
To study spin noise in interacting spin systems, we developed a quantitative theoretical approach, which is based on the stochastic path integral technique. We calculated higher order cumulants in a model of a ferromagnetically coupled interacting spin system. At temperatures below the phase transition point, we found that the 3rd spin cumulant becomes nonzero in the state with a spontaneous symmetry breaking. Approaching the phase transition point from higher temperatures, the higher order cumulants grow not only in the magnitude but also in comparison to the noise power.
Finally, by applying the stochastic path integral technique to conducting electrons we found an enhancement of the relative role of C 4 due to the Pauli exclusion principle. In the metallic phase, we predict that the dimensionless ratio of C 4 and (C 2 ) 2 is inversely proportional to temperature. Hence, by observing the spin noise at the lowest possible temperatures one may achieve the regime with strongly non-Gaussian spin fluctuations.