Zeros of partition function for Continuous Phase Transitions using cumulants
Debjyoti Majumdar Somendra M. Bhattacharjee







debjyoti@iopb.res.in
Institute of Physics, Bhubaneswar, Orissa, 751005, India
Homi Bhabha National Institute, Training School Complex, Anushakti Nagar, Mumbai 400085, India

somen@iopb.res.in
Institute of Physics, Bhubaneswar, Orissa, 751005, India




  This paper explores the use of a cumulant method to determine the zeros of partition functions for continuous phase transitions.  Unlike a first-order transition, with a uniform density of zeros  near the transition point, a continuous transition is expected to  show a power law dependence of the density with a nontrivial slope  for the line of zeros.  Different types of models and methods of  generating cumulants are used as a testing ground for the method.  These include exactly solvable DNA melting problem on hierarchical  lattices, heterogeneous DNA melting with randomness in sequence,  Monte Carlo simulations for the well-known square lattice Ising  model.  The method is applicable for closest zeros near the  imaginary axis, as these are needed for dynamical quantum phase  transitions. In all cases, the method is found to provide the basic
  information about the transition, and most importantly, avoids root  finding methods.


 Zeros of partition function for Continuous Phase Transitions
  using cumulants
    Somendra M. Bhattacharjee
    December 30, 2023
===============================================================================



§ INTRODUCTION

A phase transition is generically defined as a point of non-analyticity of any quantity of interest, like the free energy for an equilibrium system, and developing an understanding of phase transitions is associated with keywords like critical exponents, order parameter, correlation length, etc., and, not the least, the transition point <cit.>.  One of the innovative ways of locating and finding the nature of a transition is the method of finding the zeros of the partition function by allowing the
temperature or similar intensive field like quantities to be complex. This method was first incorporated by Lee and Yang.  The failure of the cluster expansion method by Mayer and coworkers in describing the properties of the liquid phase, motivated Lee and Yang to study the liquid-gas phase transition in the complex fugacity plane and the two-dimensional (2d) Ising model in the complex magnetic field plane <cit.>.  Of late, studies of zeros of partition functions, or its equivalents, have become relevant in many fields, like quantum dynamics, polymers, QCD, and analysis of experimental data, to name a few <cit.>.  The method of Lee-Yang zeros are not only confined to cases of equilibrium systems but have also been applied to systems out of equilibrium successfully <cit.>.  Although it is difficult to realize complex parameters in experiments, recently there have been attempts to explore the Lee-Yang zeros of an
Ising type spin bath <cit.>.

The basic idea of the method of zeros is to find the zeros of the partition function for a finite system.  As the partition function, Z, involves a sum of positive Boltzmann factors, there can be no real positive β=1/k_BT, T being the temperature and k_B the Boltzmann constant, where Z may vanish.  The roots are in complex conjugate pairs in the complex-β plane.  However, as the size of the system increases, the roots may have a limit point on the real axis which then represents the transition point.  The identification follows from the fact that a zero of the partition function translates into a singularity of the free energy F=-k_BTlnZ.  In addition, how the limit point is reached in the complex plane with size, the density of zeros, and the angle at which the line of zeros impinge the real axis, are related to the universality class of the transition <cit.>.  Therefore, locating the zeros in the complex plane near the real axis provides a way of characterizing a phase transition.   There are other methods of finding the transition behaviour, e.g., from the size dependence of response functions like specific heat, susceptibility, etc., with the help of data collapse assuming finite-size scaling <cit.>.  Such methods have worked nicely but often pose problems for unknown systems, with several parameters to be adjusted simultaneously.   In such situations,  the method of zeros with an underlying rigorous basis has been found to be fruitful in different problems <cit.>.


The difficulty of finding the zeros of the partition function using numerical techniques lies in the increasing level of memory requirements in the thermodynamic limit in root-finding methods.  The zeros for large systems may appear as curves, may cover some planar region, or may occur in some other complicated fashion. Although the pattern of
the zeros can be visually appealing, it is only the zeros near the positive real axis (or the imaginary axis in the case of Dynamical Quantum Phase Transition (DQPT) <cit.>) that bear useful information.  Recently, a cumulant method has been proposed to locate the leading zero for a  first order transition, taking  the zipper model
of DNA as an example <cit.>. Although the method  seems to be well behaved in case of first order transition, it calls for more detailed analysis for cases where zeros simply do not distribute uniformly, as for continuous transitions.    Our aim in this paper is to develop an iterative scheme based on the cumulant method that can
successively determine the zeros near the real axis.   

We consider primarily a model of DNA melting that can be studied exactly using real space renormalization group, namely two interacting strands of DNA on hierarchical lattices.  This particular model shows a continuous transition whose nature can be tuned by changing the dimensionality of the lattice. Effectiveness of the cumulant method from different types of data generated such as by exact methods or from Monte Carlo simulation to see to what extent they can provide us with fruitful information.  The outline of the paper is as follows.  In Sec II, we explain the method of finding zeros based on cumulants.  Sec III summarizes the model of DNA melting on hierarchical lattices.  In this section, the melting point is first determined from the closest zeros of a few lower generations (equivalently smaller system sizes).  The second closest zeros are found with some less accuracy and finally, the exponent for diverging length is calculated.  In Sec IV, this procedure is used for hetero-DNA, where the pairing energy has a binary disorder, mimicking the two possible values of base pairings. The melting point is now sample dependent and  forms a distribution. By observing how the width of the distribution varies with the system size, we conclude that there is a unique T_c in the thermodynamic limit.  Secs III and IV deal with data obtained from exact computations.  In
contrast, in Sec V, we apply the method with moments generated from Monte Carlo simulations for a 2-dimensional Ising model. In Sec VI, the method is applied to find imaginary zeros for a one-dimensional Ising model. A summary is given in Sec VII.




§  THE CUMULANT METHOD

For a finite discrete system, the partition function is Z= ∑_j  e^-β E_j, where the sum is over all states j with E_j as the energy, and as already mentioned, cannot have any real, positive β as a zero.  The zeros in the complex-β plane, satisfying Z(β)=0, occur in complex conjugate pairs.  We discuss here a method <cit.> of determining the leading pair closest to the real axis, and then show how the cumulants can also be used to determine the next pair of zeros.  Determination of a few zeros would allow us to compute the critical exponents associated with the transition.

We can define energy cumulants by taking derivatives of the logarithm of the partition function with respect to the inverse temperature β, and study the zeros in the complex-β plane.  The nth order cumulant can be defined as

    U_n = (-1)^n∂^nln Z/∂β^n.

 According to the Weierstrass factorization theorem any entire function can be written in terms of the product of its zeros <cit.>. Thus, the partition function
can be written as <cit.>

    Z(β) = Z(0)e^β c∏_k(1-β/β_k),

where c is a constant independent of β, β_k is the kth zero in the complex-β plane and these zeros occur in complex conjugate pairs. Thus, the free energy becomes a sum over all the zeros

    F =-β^-1[ ∑_kln(1-β/β_k)+lnZ(0)+β c].

The nth order cumulant is completely expressed by the zeros as

    U_n= (-1)^n-1∑_k(n-1)!/(β_k-β)^n.

The cumulants are real quantities since the zeros come in complex conjugate pairs.  Eq. (4) can be further reduced to

    U_n= 2(-1)^n-1(n-1)!∑_k^cos(n ϕ_k)/(ρ_k)^n,

where ρ_k is the distance to the kth zero from the point on the real β axis at which cumulants are being calculated and ϕ_k is the angle with the real β axis from the cumulant calculating point; see Fig. <ref>.  Note that the sum in Eq. (5) is only over the zeros of the upper half of the complex-β plane.  In the higher orders, cumulants will be dominated by the closest zero contribution only, i.e.,

    U_n≈ U_n^(0)= 2(-1)^n-1(n-1)!cos(n ϕ_0)/(ρ_0)^n, n>>1.




Thus, from four such successive cumulants in the higher orders we can solve for the leading pair of zeros, ρ_0 and ϕ_0, using the matrix equation

    [              1     -κ^(+)_n/n;              1 -κ^(+)_n+1/n+1 ][ -2ρ_0cosϕ_0;       ρ_0^2 ]
     =[ (n-1)κ_n^(-);   nκ_n+1^(-) ],

where, κ_n^(±)=U^n±1/U^n <cit.>.

Now, with the closest pair of zeros, we develop an iterative scheme to determine the next set of closest zeros which we call (ρ_1,ϕ_1).  Let us define a truncated cumulant, U_n, such that 
    U_n=U_n - U_n^(0),
 which will now be dominated by the second closest zeros in some region.  We use these U_n to define a new set of cumulant ratios κ_n.  The procedure of Eq. (7) can then be repeated with the new set to obtain the second closest zeros (ρ_1,ϕ_1). 

Implementation of the exact method of the previous paragraph is fraught with the difficulties of propagation of errors, especially when the zeros are determined successively.  In the higher order cumulants the closest zero mostly dominates, while in the lower order other zeros should contribute too. It must be somewhere in between that the closest zero takes over and just before that the largest contribution must come from the closest and the second closest zeros.  Thus, we have to look for this crossover region, which will provide us with an upper bound on the second closest zero. The upper bound, will depend upon the error introduced in calculating the closest pair of zeros and the error propagates into the approximate cumulant calculations.  Let us consider the region where the nth cumulant consists mainly of the first and the second closest zeros. Then, we may write

    U_n≈ (-1)^n-1 (n-1)! [1/(β_0-β)^n+1/(β_1-β)^n] + c.c,

where β is the point on the real axis of the complex-β plane at which the cumulant is being calculated, β_0 and β_1 being the closest and second closest zeros respectively. Consider,

    β_0= β_0^'+δβ_0,

where β^'_0 is the approximate value of the closest zero, and δβ_0 is the error introduced when calculated from the higher order cumulants.  Then, the Eq. (8) can be written as

    U_n≈(-1)^n-1 (n-1)![1/(β^'_0+δβ_0-β)^n
    
    +1/(β_1-β)^n] + c.c



    ≈ (-1)^n-1 (n-1)![1/(β^'_0-β)^n{1-nδβ_0/β^'_0-β}
    
    + 1/(β_1-β)^n]+c.c.

Subtracting the contribution of the closest zeros we obtain

    U_n= 2(-1)^n-1 (n-1)![cos(nϕ_1)/ρ_1^n
    
    -n |δβ_0|cos((n+1) ϕ_0+ψ)/ρ_0^n+1],

where ψ is the phase factor coming from δβ_0=|δβ_0|exp(-iψ).  Now, the first term in Eq. (12) dominates when

    n≪|ln|δβ_0||/ln(ρ_1/ρ_0).


Eq. (13) shows how the error in the first pair of zeros puts a constraint on the order of the cumulants that can be used for the second pair.




§ DNA ON HIERARCHICAL LATTICES

To investigate the effectiveness of the cumulant method for continuous transitions, we first consider to apply it for a model of double stranded DNA (dsDNA) on hierarchical lattice.  Hierarchical lattice with discrete scaling allows exact implementation of the real space renormalization (RG).  Thus one can write down recursion relations for the partition function and for the zeros of the partition function.  This allows us to compare the cumulant method with the exact results from RG.

On untwisting the helical structure of a dsDNA, it takes the shape of a railway track where the ties represent the bonds shared by the two strands. This structure of bond sharing can be mimicked by putting the two strands on a hierarchical lattice with the two endpoints fixed but they can wander at intermediate points.  The lattice is generated iteratively by replacing each bond in the (n-1)th generation with a motif of λ b bonds to get the new nth generation, where λ and b are the bond scaling factor and the bond branching factor respectively; see Fig.  <ref>.  The effective lattice dimension in the thermodynamic limit will be given by <cit.> 
    d=lnλ b/lnλ.


In this paper we shall choose λ=2 and b=4 which corresponds to a three dimensional lattice.  The contact energies are defined to be -ϵ (ϵ > 0) when two strands share a bond.  This constitutes our DNA model <cit.>. 







The partition function of our model can be written in terms of the Boltzmann factor as the variable y=e^βϵ, determined by the basic energy scale of the problem.   For a finite system, the partition function is a polynomial in y of the form 

    Z_m = a_0+a_1y^2+⋯+a_ny^2^m.

  We define cumulant-like quantities by taking derivatives of the logarithm of the partition function with respect to the Boltzmann factor i.e. y instead of β, and study the zeros in the complex-y plane.  The nth order cumulant can be
defined as

    U_n = (-1)^n∂^nln Z/∂ y^n.

  Exact recursion relations followed by the single chain (C_n) and double chain (Z_n) partition functions can be written as  <cit.>

    C_n=b C_n-1^2,


    Z_n=b(b-1)C_n-1^4+b Z_n-1^2,

with the initial conditions taken as

    C_0=1,  Z_0=y.

The zeros of the partition function thus follows the following
recursion relation

    q=±√(bq̃_̃j̃-(b-1)),

where q̃_̃j̃ is one of the zero of the (n-1)th generation, giving rise to two new zeros in the nth generation <cit.>.


It is known that dsDNA undergoes a phase transition from a bound state to an unbound state (dsDNA melting) under variation of temperature, the critical point being known as the duplex melting
point.  From RG relations the critical point is found to be y_c=b-1, with y being the Boltzmann factor.

We used the method of Sec. <ref> to determine the closest zeros.  To do so, the cumulants were calculated at y=2.27 and Eq. (7) used to determine the real and imaginary parts of the closest zeros.  The results obtained for a few smaller generations N=5 to  11 are shown in Fig. <ref>.  As the figure shows, on extrapolation to infinite size, the values are in good agreement with the known exact results.  The iterative scheme of Sec. <ref> is then used to determine the second closest zero.  In Fig. <ref>(a) and Fig. <ref>(b) we plot the real and imaginary parts of the second closest zero calculated at three different temperatures.  Clearly, we did not consider very low order points since all the zeros contribute there.  There appears to be a crossover region after which there is a constant contribution which corresponds
to the first closest zero. In between, there is a flat plateau which is the region of our interest. Thus, we fit our curve such as to get all the points before the crossover region and up to the flat plateau to find our best value of the second closest zero.  The thickness of the shaded region in the flat plateau should serve as an error bar to our value.  While the closest zeros can be determined with a decent accuracy of four decimal places for smaller generations, the next to closest pair of zeros can be calculated with somewhat less accuracy of one place after the decimal.

Next, we find the angle at which the zeros meet the real axis in the thermodynamic limit; see Fig.  <ref>(c). For this, we have taken the average of the angles made by the first and second closest zeros with the critical point (since the zeros do not fall in a smooth straight line, in lower generations) and plotted it with the generation number inverse.  Our estimate for the meeting angle of the zeros with the real y axis θ_estimate=53.87, compares well with the exact result θ_exact=52.646.  The results vary substantially if we
do not include the second closest zero.

This angle contains the information about the order of the phase transition. The angle at which the zeros meet the real y axis is related to the specific heat exponent as

    tan(θν)=-tan(πα)+A_-/A_+(πα),

where θ is the angle between the tangent of the zeros at the limit point and the real axis of y, ν is the critical exponent for diverging length, α is the specific heat exponent and A_± are the amplitudes of the specific heat on the low and high y side of the transition <cit.>. For our problem with double stranded DNA, it is known that A_-/A_+→∞ as A_+=0. Therefore, the angle is given by

    θ=π/2ν.

Once we have θ, we can find the critical exponent ν. But in general we need to know ν separately. To determine ν one can make use of the RG transformation properties of the Hamiltonian to see how the partition function would change under a scaling <cit.>, and obtain the following scaling equation of the distance of the kth zero from the critical point

    Δ_k=k^1/dνL^-1/νf^-1,



where d is the dimension, Δ_k is the distance of the kth zero from the critical point, L is the size of the system and f^-1 is a constant. To determine ν we plot lnΔ_k vs ln L, for the closest zeros (k=1) for generations N=5,6,7,8; see Fig.  <ref>(d).



Our estimate of the length exponent ν_estimate=1.67 compares reasonably well with the independent estimate from  Eq. (22) which also gives ν=1.67.  The exact value of ν_exact=1.70951.  It might be possible to improve the accuracy of the result by iterative inclusion of other zeros.



§ EFFECTS OF  BINARY RANDOM DISORDER

In order to show the usefulness of the cumulant method, we now use it for a heterogeneous DNA.  The non-uniformity in the base sequence of a DNA can be handled by specifying the base pair energies along the chain. A simpler situation is a random sequence where the interaction energy is randomly chosen from a specified distribution <cit.>.

With the introduction of randomness, the obvious question is whether a critical point for melting exists. Since the partition function is different for different realization of randomness, any thermally averaged quantity will have a probability distribution.  One such important distribution is the probability distribution of the closest zero. In the thermodynamic limit, a well-defined
melting point would require a sharp probability distribution.  The width of the distribution should vanish in the large length limit ("self-averaging").  Although, that is not enough for a transition to occur the zeros must touch the real axis.  However, this is not equivalent to the traditional quenched averaging.  The quenched average transition temperature can be determined from the average
specific heat as we discussed below.

As a model system, we take two different types of pairing energy, ϵ_1=1.0 and ϵ_2=0.5  with equal probability.  The randomness depends only on the longitudinal direction and has no dependence on the transverse direction.  For the closest zeros of each sample,  ρ and cosϕ are then determined by the using the method of Sec. <ref>  from the four successive cumulants of order n=19,20,21,22.  The cumulants were calculated at β=ln2.75 and the whole procedure is repeated for 200 random samples.  Fig.  <ref>(a) shows that the average  cosϕ, where the overbar represents an averaging  over samples, approaches 1 on extrapolation.  This  indicates that there is a limit point on the real β axis. see A similar extrapolation of ρ in Fig.  <ref>(b) locates the transition point at  β_c=1.23±0.004.

Fig.  <ref>(c) and Fig.  <ref>(d) show the distributions of ρ and cosϕ for successive generations of N=6,7,8.  The width of the joint distribution is found to go to zero in the thermodynamic limit.  An independent way of determining the transition temperature is to consider the average specific heat or the second energy cumulant U_2 for various sizes of the system and look for the common crossing point.  Fig.  <ref>(e) shows the second energy cumulant U_2 curves (scaled by the length of DNA) for three different generations, N=7,8,9, each one averaged over 500 samples. The intersection of the U_2 curves gives the melting temperature to be β_c=1.35±0.03 (transparent blue circled region), where the uncertainty is based on the sample fluctuations of the intersection point.  For comparison, the melting point for a pure system with ϵ_1 is β_c(ϵ=1)=1.0986, while for the other one, β_c(ϵ=0.5)=2.1972. The zeros seem to be dominated by the tighter pairs.  For small lengths, it is possible to get samples with over representation of the lighter energy bonds that led to a double peaked distribution.  As there are situations like Griffiths' phase <cit.> where zeros of rare regions play a crucial role, this cumulant method of finding zeros will be relevant there.




§ ISING 2D

In the previous sections <ref> and <ref>, we determined the closest zeros as well as the second closest zeros from the higher energy like cumulants of the partition function. There we had the privilege of finding the exact cumulants since the partition function is known exactly. But in simulations, the higher energy cumulants can be calculated more easily, albeit approximately, than finding the partition function.  One such example is the two dimensional Ising model. The Hamiltonian for 2d Ising model with interaction between the nearest neighbours in the absence of an external magnetic field is

    βℋ=-J∑_⟨ ij⟩σ_iσ_j,

where σ_i∈{-1,1} is the spin at site i, J is the coupling constant between neighboring spins and ⟨ ij⟩ indicates nearest neighbors.  It is known that for such a two-dimensional spin system there is a phase transition from paramagnetic to ferromagnetic phase at a finite temperature.  The energy moments were obtained from Monte-Carlo simulation using the Metropolis algorithm <cit.>  and then  MATHEMATICA was used to determine the cumulant <cit.>. Since our aim throughout this work had been to determine quantities in the thermodynamic limit by studying the zeros of small system size, we have calculated the leading zeros for N× N square lattices with N = 5,6,7,8,9,10,11 using the cumulant method.  To find the zeros we have calculated up to 30th moment of energy using quadruple precision data type at β=1/3 with J/k_B=1.







    β_c     β_c     ν 
 3*Estimated
            

     0.445     0.003     1.01 

            


 3*Exact
            

     0.44068     0     1 

            

 
tableComparison with exact values for 2d Ising model.





Fig.  <ref>(a,b) show the real and the imaginary parts of the leading zero.  Fig.  <ref>(b) shows that the imaginary part of the closest zero goes to zero in the thermodynamic limit as required for a phase transition.  The critical temperature is then given by the extrapolated value in Fig.  <ref>(a). Fig.  <ref>(c) shows how the distance between the closest
zeros and the transition point scales with the system size.  The inverse of the slope of the fitting line gives the critical exponent ν for the diverging correlation length according to the Eq (23).  In Table <ref> we have compared our numerical results with the exact results in the thermodynamic limit.




§ IMAGINARY ZEROS IN ISING 1D


Motivated by the fact that a phase transition can happen with time as a parameter in case of DQPT where the zeros meet on the imaginary axis <cit.>, we have found the zeros closest to the imaginary axis, for the 1d Ising model which do not have any transition at any non zero temperature.  The partition function of the 1d Ising model is exactly known to be

    Z = λ_+^N+λ_-^N,

where λ_± are the eigenvalues of the transfer matrix and N is the number of spins.

In absence of any external magnetic field, the eigenvalues become

    λ_±=e^κ± e^-κ,

where κ=J/k_BT and the zeros of the partition function on the complex-κ plane are found to be

    κ_j=1/2log[-itan((2j+1)π/2N)],

with j=0,1,2⋯ N-1. We here consider only even values of N for which the zeros in the complex-κ plane are distributed symmetrically about the imaginary axis.   The zeros meet the
imaginary axis periodically due to the presence of the -i factor inside log in Eq. (27); see Fig.  <ref>(a).


Assuming a symmetric distribution about the imaginary axis, we can define a free energy like quantity as 

    F=ln(Z Z^*),

where Z and Z^* are the partition functions evaluated at a complex temperature and its complex conjugate respectively so that F is real.  For pure imaginary κ, there is a first order transition at |κ|=π/4, where |λ_+|=|λ_-|, and odd multiples of it.  These transition points are the zeros given by Eq. (<ref>), and are related to
DQPT <cit.>.

In terms of zeros of the complex-κ plane

    F=ln[∏_j(κ_j-iκ) ∏_l(κ_l+iκ)].

The n-th cumulant is given by

    U_n   =   (-∂_iκ)^nln (ZZ^*)
       =   (-∂_iκ)^nln[∏_j(κ_j-iκ) ∏_l(κ_l+iκ)].

Retaining only the terms of the closest zeros, the above equation can be written as

    U_n=4(-1)^n-1(n-1)!cos(nϕ_0)/(iρ_0)^n,

which requires a little modification of the matrix equation Eq. <ref> by simply putting iρ_0 instead of ρ_0. Here, ϕ_0 is the angle made with the imaginary axis by the vector
extending from the cumulant calculating point on the imaginary axis.  Because of the periodic pattern of the zeros, care is needed in choosing the cumulant calculating point.  Fig.  <ref>(b) shows the real and imaginary parts of the closest zeros of the branch closest to the real κ axis on the upper half of the complex-κ plane, for different system sizes.  We do recover the closest zero properly.




§ SUMMARY

To summarize, we developed an iterative scheme to determine the zeros of a partition function using cumulant data. The method has been tested for DNA melting transition on a hierarchical lattice for  which the zeros can be determined  exactly.  We  found the leading pair of zeros of the partition function both for pure and hetero-DNA cases.  The next leading pair of zeros are also found which proved to be important in determining the angle at which the zeros approach the real axis and also in determining the critical exponent.  We have chosen such a system to
have a good comparison with the exact results. The results found out to be matching well enough with the theoretical exact results, thus giving us a new way to look into systems in the thermodynamic limit from smaller system sizes. Although the zeros from the lower generations may not always be helpful, still it proved to be quite a powerful method. Moreover, it is shown that the method is well-suited for results from Monte Carlo simulation proved to be equally good to reproduce the known results for the critical exponent and critical point estimation for the case of Ising model in 2d.  An extension of the method has been shown to be useful for the imaginary zeros too.



§ ACKNOWLEDGEMENTS
 
D.M. would like to thank Jaya Maji and Goutam Tripathy for useful discussions.   Computations performed on the SAMKHYA high performance computing facility at Institute Of Physics Bhubaneswar.

99

goldenfeldN. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group, 1st ed. (Westview Press, Boulder,1992).


leeyang01 C. N. Yang and T. D. Lee, Phys. Rev. 87, 404 (1952).


leeyang02T. D. Lee and C. N. Yang,  Phys. Rev. 87, 410 (1952).


heylM. Heyl, A. Polkovnikov, and S. Kehrein, Phys. Rev. Lett. 110, 135704 (2013).


rensburg E J Janse van Rensburg, J. Stat. Mech. (2017) 033208.


nagata K. Nagata, K. Kashiwa, A. Nakamura, and S. M. Nishigaki, Phys. Rev. D 91, 094507 (2015).


evans R. A. Blythe and M. R. Evans, Braz. J. Phys. 33, 464 (2003).


weiB.-B. Wei and R.-B. Liu, Phys. Rev. Lett. 109, 185701 (2012).


pengX. Peng, H. Zhou, B.-B. Wei, J. Cui, J. Du, and R.-B. Liu, "Phys. Rev. Lett. 114, 010601 (2015).


brandner K. Brandner, V. F. Maisi, J. P. Pekola, J. P. Garrahan, and C. Flindt, Phys. Rev. Lett. 118, 180601 (2017).


itzykson C. Itzykson, R.B. Pearson and J.B. Zuber, Nuclear
  Physics B220[FS8] 415 (1983).


stanley H. E. Stanley, Rev. Mod. Phys. 71, S358 (1999).


flindtbig A. Deger, K. Brandner, and C. Flindt, Phys. Rev. E 97, 012115 (2018).


conway J. B.  Conway, Functions of One  Complex  Variable (Springer, New York, 1995).


jmsmb J. Maji and S.M. Bhattacharjee, Phys. Rev E 86, 041147 (2012).



smsmb2 S. M. Bhattacharjee and S. Mukherji, Phys. Rev. Letts. 70, 49 (1993).


smsmb3 S. Mukherji and S. M. Bhattacharjee, Phys. Rev. E48, 3483 (1993).


smsmb1 S. Mukherji and S.M. Bhattacharjee, Phys. Rev E 52, 1930 (1995).


monthus C. Monthus and T. Garel, Eur. Phys. J. B 48, 393 (2005).


griffiths R. B. Griffiths, Phys. Rev. Lett. 23, 17 (1969).

montecarlo  Binder  K  and  Heermann  D  W  1988 Monte  Carlo  Simulation  in  Statistical  Physics.  An  Introduction (Berlin:Springer)


mathematica Mathematica command to calculate cumulants from moments: .


