Critical exponents in mean-field classical spin systems
CCritical exponents in mean-field classical spin systems
Yoshiyuki Y. Yamaguchi
Department of Applied Mathematics and Physics, Graduate School of Informatics, Kyoto University, Kyoto 606-8501, Japan
Debraj Das ∗ and Shamik Gupta Department of Physics, Ramakrishna Mission Vivekananda University, Belur Math, Howrah, 711202, India (Dated: September 24, 2019)For mean-field classical spin systems exhibiting a second-order phase transition in the stationarystate, we obtain within the corresponding phase space evolution according to the Vlasov equationthe values of the critical exponents describing power-law behavior of response to a small externalfield. The exponent values so obtained significantly differ from the ones obtained on the basis of ananalysis of the static phase-space distribution, with no reference to dynamics. This work serves asan illustration that cautions against relying on a static approach, with no reference to the dynamicalevolution, to extract critical exponent values for mean-field systems.
I. INTRODUCTION
Since early days of statistical mechanics, studyingphase transitions in physical systems has been a theme ofactive research in the field. Phase transitions can occuronly in the thermodynamic limit. Second-order or contin-uous phase transitions are characterized by a power-lawbehavior of macroscopic quantities close to the criticalpoint of transition. Such transitions in different systemsmay be broadly classified into universality classes iden-tified by different values of critical exponents describingthe power-law behavior. For example, for a ferromagnetexhibiting a second-order phase transition as a functionof temperature T , the magnetization close to and on thelower side of the critical point T c has a power-law de-pendence on the separation ( T c − T ) from the criticalpoint, with the corresponding exponent being β . On ap-plying an external field, the magnetization increases asa function of the field strength, and in the limit of aninfinitesimal field, a linear growth for T (cid:54) = T c implyinga linear response determines the zero-field susceptibility χ . The susceptibility diverges as a power law close toand on both sides of the critical point, with the corre-sponding exponents denoted by γ + and γ − on the dis-ordered ( T > T c ) and the magnetized ( T < T c ) phase,respectively. At the critical point, the response becomesnonlinear, being characterized by the critical exponent δ .These critical exponents are known to satisfy the scalingrelation γ ± = β ( δ −
1) [1–3].One representative class of systems exhibiting second-order phase transitions is that of mean-field systems. Inthermal equilibrium of such systems, statistical mechan-ical predictions for the critical exponents, based on ananalysis of the thermal equilibrium phase space distri-bution with no reference to dynamics, yield the values β = 1 / , γ ± = 1, δ = 3 [2]. However, owing to themean-field nature of the time evolution, critical expo-nents obtained on the basis of dynamics may well have ∗ Corresponding author: [email protected] different values. Indeed, dynamics of a mean-field sys-tem in the thermodynamic limit is described by the so-called Vlasov equation that allows a vast number of sta-ble stationary states, and thermal equilibrium is just oneof them [4–7]. This implies that once the system is in astable stationary state other than thermal equilibrium,it would not relax to thermal equilibrium. A large butfinite system remains trapped in so-called quasistation-ary states (QSSs) identified as stable stationary solutionsof the Vlasov equation, with finite-size effects allowing aslow evolution of the QSSs towards thermal equilibriumover a timescale that diverges with the system size [8, 9].Existence of QSSs allows nonequilibrium phase transi-tions: a generic initial state undergoes a violent relax-ation to relax to a QSS [10], and the nonequilibriumphase transition can for example be defined with re-spect to the value of the order parameter in the QSS. Ina given system, these nonequilibrium phase transitionsmay not necessarily be continuous even when the equi-librium phase transition is continuous, and several dis-continuous nonequilibrium phase transitions have beenreported in the literature [11–14]. In this article, we,however, focus on families of QSSs that exhibit contin-uous phase transitions and to which an external field isapplied in order to investigate the values of the criticalexponents characterizing the response.The aforementioned trapping scenario holds even whenan external field is applied to the system prepared in athermal equilibrium state: With the field on, a finite sys-tem goes from the initial to a new thermal equilibriumstate via intermediate QSSs, while a thermodynamic sys-tem remains trapped in a QSS and does not relax tothermal equilibrium [15]. The latter fact requires thatone invokes an alternative strategy of obtaining suscep-tibility that is based on the Vlasov dynamics when ad-dressing the issue of response of mean-field systems inthermal equilibrium to an external field. The criticalexponents γ ± so obtained may not necessarily coincidewith the ones computed within equilibrium statisticalmechanics. Indeed, in the so-called Hamiltonian mean-field (HMF) model [16, 17], a paradigmatic mean-fieldsystem exhibiting a second-order phase transition, the a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p critical exponents obtained within the Vlasov dynam-ics have been shown to be γ + = 1 , γ − = 1 / β = 1 / δ = 3 /
2. More generally, the critical expo-nents within the Vlasov dynamics have been obtainedas γ + = 2 β, γ − = β/ , δ = 3 / γ − = β ( δ − γ ± . This fact is borne out bythe values of the critical exponents obtained in the HMFmodel by taking into account the Casimir constraints.Another important remark is that divergence of suscepti-bility is observed even when the dynamics is non-ergodic,as is found in the case of the HMF model [15].The HMF model mimics the classical XY model, withan additional kinetic energy term assigned to individualspins. Owing to the latter whose range is the whole realset, the one-particle phase space of the HMF model isa cylinder. In the HMF model, the Poisson bracket be-tween the spin components is taken to vanish identically,In this work, we consider Heisenberg spin systems withmean-field interactions, in which the Poisson brackets be-tween the spin components are strictly nonzero, and thesingle-particle phase space is the unit sphere. Consid-ering the time evolution of the spin components accord-ing to a Hamiltonian with a mean-field interaction anda local anisotropy, we address here several questions oftheoretical and practical relevance: Does the universalityclass for usual Hamiltonian systems defined on a cylin-der, e.g., the HMF model, include spin systems definedon the unit sphere? What is the effect of the anisotropyon the critical exponents? Would the scaling relation γ − = β ( δ −
1) still hold even if the spin system is foundto be in a different universality class?This paper is organized as follows. The spin modelwe study is introduced in Sec. II. Here, the dynamicsdescribed by the canonical equations of motion is alsodiscussed, as is the characterization of the dynamics inthe thermodynamic limit in terms of the Vlasov equa-tion. Based on the latter, we discuss the setting andthe definition of the critical exponents in Sec. III, whileour theoretical predictions for the critical exponents arederived in Sec. IV. Detailed numerical checks of our the-oretical predictions are pursued in Sec. V. Section VIconcludes the paper with discussions.
II. THE MODELA. Definition
Our model of study consists of N globally-coupled clas-sical Heisenberg spins of unit length denoted by S i = ( S ix , S iy , S iz ); i = 1 , , . . . , N. (1)The N -body Hamiltonian of the model is given by H N = − J N N (cid:88) i,j =1 S i · S j + D N (cid:88) i =1 S niz − h ( t ) · N (cid:88) i =1 S i . (2)Here, the first term with J > J has been scaled downby the system size N in order to make the energy exten-sive, in accordance with the Kac prescription [21]. Thesystem (2) is however intrinsically non-additive: it can-not be trivially subdivided into independent macroscopicparts. In the following, we set J = 1 without loss of gen-erality.In Eq. (2), the second term with D > S iz → − S iz ,we have made here the choice of even exponent equal to2 n , with n being a non-negative integer. We refer to themodel with exponent 2 n as Model- n . Note that Model-0 is completely isotropic in the spin space, and there isno preferred direction of orientation of spins. Model-1has been studied previously in the context of QSSs inRefs. [22, 23]. Model-2 is the special case of a quarticanisotropy; it may be noted that thermodynamic prop-erties of a Heisenberg spin model containing a quarticterm have been studied in Ref. [24].The third term on the right-hand side of Eq. (2) arisesdue to the application of a time-dependent external mag-netic field h ( t ) ≡ ( h x ( t ) , h y ( t ) , h z ( t )). In this work, weconsider the external field to be absent for times previ-ous to instant t , i.e., for times t < t , when the systemwill be assumed to be existing in a reference state, e.g., athermal equilibrium state. For times t ≥ t , on the otherhand, we would put on a constant field in order to mea-sure the response of the reference state to the externalfield. The explicit form of h ( t ) is thus given by h ( t ) = Θ( t − t ) h , (3)where Θ( t ) is the unit step function, and h is a vector ofconstant length equal to h . The singularity of the unitstep function Θ( t ) will have no effect on the values ofthe critical exponents obtained based on the Vlasov dy-namics, and we may replace Θ( t ) with a smooth function[25]. B. Spin dynamics
In dimensionless times, the time evolution of system(2) is governed by the set of coupled first-order differen-tial equations˙ S i = { S i , H N } ; i = 1 , , . . . , N, (4)where the dot denotes derivative with respect to time.The Poisson bracket {· , ·} is bilinear, skew-symmetric,and satisfies the Leibniz’s rule { XY, Z } = { X, Z } Y + X { Y, Z } (5)for any functions X, Y, and Z of the spins. The Poissonbrackets between two spins are given by { S ix , S jy } = δ ij S iz , { S iy , S jz } = δ ij S ix , { S iz , S jx } = δ ij S iy . (6)Using Eqs. (2), (4), and (6), we obtain the time evolutionof the spin components as˙ S ix = S iy ( m z + h z ) − S iz ( m y + h y ) − nDS iy S n − iz , ˙ S iy = S iz ( m x + h x ) − S ix ( m z + h z ) + 2 nDS ix S n − iz , ˙ S iz = S ix ( m y + h y ) − S iy ( m x + h x ) , (7)where m ≡ N N (cid:88) i =1 S i = ( m x , m y , m z ) (8)is the magnetization vector that serves as the mean fieldgoverning the time evolution of the individual spins.Summing the third equation of (7) over i , we find that m z is a constant of motion if the condition m x h y − m y h x = 0 (9)is satisfied. The length of each spin is a constant of mo-tion, and so is the total energy of the system when thefield h is time independent.Writing the spin components in terms of spherical po-lar angles θ i ∈ [0 , π ] and φ i ∈ [0 , π ), as S ix = sin θ i cos φ i , S iy = sin θ i sin φ i , S iz = cos θ i , (10)we obtain from Eq. (7) the time evolution of the variables θ i and φ i as˙ θ i = ( m x + h x ) sin φ i − ( m y + h y ) cos φ i , ˙ φ i = ( m x + h x ) cot θ i cos φ i + ( m y + h y ) cot θ i sin φ i − ( m z + h z ) + 2 nD cos n − θ i . (11) For later convenience, we introduce a new variable p i ≡ cos θ i , in terms of which we have S ix = (cid:113) − p i cos φ i , S iy = (cid:113) − p i sin φ i , S iz = p i . (12)In terms of p i , which is in fact canonically conjugate to φ i , the Poisson bracket {· , ·} reads [22] { X, Y } = N (cid:88) i =1 (cid:18) ∂X∂φ i ∂Y∂p i − ∂X∂p i ∂Y∂φ i (cid:19) . (13)The dynamical variables of the i -th spin are thus φ i and p i , while a volume element in the ( φ i , p i )-space is d φ i d p i . C. Description in the thermodynamic limit
In the thermodynamic limit N → ∞ , the dynamics ofsystem (2) is described by the Vlasov equation ∂f∂t + ∂H∂p ∂f∂φ − ∂H∂φ ∂f∂p = 0 , (14)where f ( φ, p, t ) is the single-spin distribution functionthat measures the probability density to find a spin ( φ, p )at time t , while the single-spin Hamiltonian H is H ( φ, p, t ) = Dp n − [ m ( t ) + h ( t )] · S , (15)with S ≡ (cid:16)(cid:112) − p cos φ, (cid:112) − p sin φ, p (cid:17) , (16)and the magnetization vector m = ( m x , m y , m y ) givenby m ( t ) = (cid:90) (cid:90) µ S f ( φ, p, t )d φ d p. (17)The double integral over any function X ( φ, p ) in thesingle-spin phase space µ ≡ ( φ, p ) is defined as (cid:90) (cid:90) µ X ( φ, p )d φ d p ≡ (cid:90) π d φ (cid:90) − d p X ( φ, p ) . (18)Note that the single-spin Hamiltonian (15) dependson time t through the magnetization m ( t ) and theexternal field h ( t ). Normalization of f ( φ, p, t ) reads (cid:82)(cid:82) µ f ( φ, p, t ) = 1 for any time t .Any quantity C [ f ]( t ) = (cid:90) (cid:90) µ c ( f )d φ d p (19)is a constant of motion for any smooth function c , asmay be seen by considering the time variation of C andusing Eq. (14). These invariants of motion are calledCasimir invariants, which hold even when the single-spinHamiltonian depends on time. The Casimir invariantsdo not allow an initial state with h = to relax to thethermal equilibrium state with h (cid:54) = when at least oneof the Casimir invariants C [ f ] between the two states isnot the same. III. SETTING AND DEFINITION OF THECRITICAL EXPONENTSA. Setting
For t < t , we consider system (2) to be existing inone of a family of stable stationary states with externalfield h = . In order that we may study the criticalexponents associated with the response of the system toan external field that we put on for times t ≥ t , werestrict to a family of states that allow a second-orderphase transition and consequently a critical point in thestationary state. We refer to such a family of states asour reference states and denote the states by f . FromEq. (14), it is evident that f of the form f ( φ, p ) = F ( H ( φ, p )) = G ( H ( φ, p )) (cid:82)(cid:82) µ G ( H ( φ, p ))d φ d p , (20)with G an arbitrary function, is a stationary solution ofthe Vlasov equation, and we have H ( φ, p ) = Dp n − m · S , (21)and m = ( m x , m y , m z ) (22)satisfying the self-consistent equation m = (cid:90) (cid:90) µ S f ( φ, p )d φ d p. (23)The family of functions G may be parametrized by aparameter T , which in the case of thermal equilibriumcoincides with the temperature [26]: G ( x ) = exp( − x/T ) . (24)However, the analysis presented in the following appliesto other family of functions G , such as the Fermi-Dirac-type family G ( x ) = 1exp[( x − a ) /b ] + 1 . (25)In this case, the parameter T may be identified with ei-ther of the two parameters a and b .Now, from the rotational symmetry of H ( φ, p ) on the( S x , S y )-plane, we may set m y = 0 without loss of gen-erality. Moreover, we may assume m z = 0, which solvesthe self-consistent equation for m z . Denoting m x by m , so that m = ( m , , H ( φ, p ) = Dp n − m (cid:112) − p cos φ, (26)while the self-consistent equation (23) reads m = (cid:90) (cid:90) µ (cid:112) − p cos φF ( H ( φ, p ))d φ d p. (27) At t = t , we turn on a constant external field h =( h, ,
0) pointing in the direction of the reference mag-netization m = ( m , , f to a stationary state f h with magnetization m h = ( m h , , f h is H h ( φ, p ) = Dp n − ( m h + h ) (cid:112) − p cos φ. (28)We stress that f h is not necessarily the thermal equilib-rium state proportional to exp[ − H h ( φ, p ) /T ], and thuscould be an out-of-equilibrium state, see Sec. II C.Within the Vlasov dynamics, the response to the externalfield is measured by δm ≡ m h − m . (29)In the above setting, we recall the definitions of the crit-ical exponents β, γ + , γ − and δ given in any standard ref-erence on critical phenomena, e.g., Ref. [1]. B. Definition of the critical exponents
The critical exponent β is defined with respect to thereference state, as m ( T ) ∝ ( T c − T ) β ; T → T − c , (30)where T c is the critical point. Here, m is the positive so-lution of the self-consistent equation (27), and the valueof β may depend on the choice of the family F of the ref-erence state. The self-consistent equation (27), however,implies quite generally that β = 1 /
2, see Appendix A.The critical exponents γ ± are defined in the regime oflinear response. The response δm depends on T and h ,and the susceptibility χ ( T ) is defined as χ ( T ) ≡ ∂ ( δm ) ∂h (cid:12)(cid:12)(cid:12)(cid:12) h → . (31)The susceptibility diverges at the critical point T c as χ ( T ) ∝ (cid:40) ( T − T c ) − γ + ( T → T + c ) , ( T c − T ) − γ − ( T → T − c ) , (32)which defines the exponents γ ± .At the critical point T c , one has δm ∝ h /δ , ( T = T c ) , (33)which defines the critical exponent δ . Usually, one has δ >
1, since the leading response is nonlinear and isstronger than the linear response.
C. Statistical mechanics predictions for the criticalexponents
Statistical mechanics analysis that considers studyingthe equivalent of Eqs. (20), (21) and (23) in presenceof a constant field h = ( h, , β = 12 , γ ± = 1 , δ = 3 , (34)irrespective of the value of the exponent n , see AppendixA for details. In the next section, we derive the values ofthe critical exponents within the Vlasov dynamics. Wewill obtain the response δm within the Vlasov dynamics,and hence the exponents γ ± and δ may very well takevalues different from the ones in Eq. (34). IV. THEORETICAL PREDICTIONS FOR THECRITICAL EXPONENTS BASED ON THEVLASOV DYNAMICS
In this section, we derive our results for the criticalexponents based on the Vlasov dynamics. As alreadymentioned above, we have β = 1 / F of the reference state. Inthe following, we discuss the computation of the criticalexponents γ ± and δ for a given value of the exponent β . A. Model- The single-spin Hamiltonian of Model-0 is H h = − ( m h + h ) S x . (35)The equations of motion are obtained from Eq. (7) as˙ S ix = 0 , ˙ S iy = S iz ( m h + h ) , ˙ S iz = − S iy ( m h + h ) . (36)Clearly, S ix for any i and consequently m x are constantof motion for any external field h irrespective of its timedependence. Dynamically, each spin rotates on a S x =constant plane. The fact that the variable S x is a con-stant of motion for both cases of h = 0 and h (cid:54) = 0 impliesthat the reference state f = F ( H ) = F ( − m S x ) is sta-tionary even after the external field is turned on, and wehave m h = m . Consequently, no response to the exter-nal field is obtained within the Vlasov dynamics. If wehave to assign values to the critical exponents, we maysay β = 12 , γ ± = 0 , δ = 1 , (37)owing to the fact that no divergence of the susceptibilityis obtained within the Vlasov dynamics. The aforemen-tioned exponent values are quite different from the onesobtained within statistical mechanics, Eq. (34). For the case G ( x ) = exp( − x/T ), the critical point is T c = 1 / G ( x ) into the function A ( T ) defined by Eq. (A6) and then solving A ( T c ) = 0).The isotropic spin model, Model-0, shows no responseto the external field and thus provides a simple and ex-treme example of dynamical suppression of response, butModel- n with n ≥ n ≥ B. Model- n with n ≥ To compute the values of the critical exponents γ ± and δ , our task is to obtain within the Vlasov dynamicsstarting from the state f the asymptotic state f h andhence the response δm . For this purpose, we remarkthat the Hamiltonian H h given in Eq. (28) is integrableand has the associated angle-action variables ( w, I ). Inthis setting, the response formula f h ( I ) = (cid:104) f ( φ, p ) (cid:105) h (38)has been proposed for Hamiltonian systems [18, 25],where (cid:104)·(cid:105) h is defined as the average over the angle vari-able w , as (cid:104) A (cid:105) h ≡ π (cid:90) π A ( φ ( w, I ) , p ( w, I )) d w. (39)In Appendix B, we summarize the derivation of Eq. (38).Using H h = H h ( I ) and (cid:104) ϕ ( I ) (cid:105) h = ϕ ( I ) (40)for any function ϕ , and the expansion of the single-spinHamiltonian as H h = H + δH, δH = − ( δm + h ) S x , (41)with S x = (cid:112) − p cos φ , we have the expansion of f h ( I )as f h = f − ( δm + h ) [ S x − (cid:104) S x (cid:105) ] F (cid:48) ( H ) − ( δm + h ) [ (cid:104) S x (cid:105) F (cid:48) ( H ) − (cid:104) S x (cid:105) h F (cid:48) ( H h )] . (42)Note that, for instance, one has F (cid:48) ( x ) = − F ( x ) /T forthermal equilibrium reference state (24). The average (cid:104)·(cid:105) is defined as an average over the angle variable asso-ciated with the integrable system H .Multiplying Eq. (42) by S x = (cid:112) − p cos φ and thenintegrating over φ and p , we have the self-consistent equa-tion for the response δm as L ( δm + h ) + N ( δm + h ) − h = (higher order terms in h ) , (43)where the coefficient L of the linear part is L ( T ) = 1 + (cid:90) (cid:90) µ (cid:2) S x − (cid:104) S x (cid:105) (cid:3) F (cid:48) ( H ( φ, p ))d φ d p, (44)while N concerns the leading nonlinear part: N ( T ) = (cid:90) (cid:90) µ (cid:2) (cid:104) S x (cid:105) F (cid:48) ( H ) − (cid:104) S x (cid:105) h F (cid:48) ( H h ) (cid:3) d φ d p. (45)The linear part L gives the values of the critical expo-nents γ ± , while the nonlinear part N gives the value ofthe exponent δ .Note that L ( T c ) = 0, so that the contribution of onlythe nonlinear response appears at the critical point T c .Away from the critical point, it is the linear part thatgives the dominant contribution in Eq. (43), so that ne-glecting the nonlinear contribution, we have the linearresponse δm = 1 − LL h. (46)Then, within linear response, the divergence of the sus-ceptibility, χ = 1 − LL , (47)is determined by the convergence behavior of L as T → T ± c .
1. Linear response in the disordered phase
In the disordered phase, the angle variable w is nothingbut φ , and we have (cid:104) S x (cid:105) = (cid:104) (cid:112) − p cos φ (cid:105) = 0 . (48)This result implies that L has no contribution from thedynamics, and hence may be expanded in a Taylor seriesin ( T − T c ) around T c , resulting in its convergence be-ing proportional to T − T c . The critical exponent γ + is,therefore, given by γ + = 1 . (49)
2. Linear response in the ordered phase
In the ordered phase, let us divide L into the two parts: L ( T ) = L ( T ) + L ( T ) , (50)with L ( T ) ≡ (cid:90) (cid:90) µ S x F (cid:48) ( H ( φ, p ))d φ d p (51)and L ( T ) ≡ − (cid:90) (cid:90) µ (cid:104) S x (cid:105) F (cid:48) ( H ( φ, p ))d φ d p. (52)The behavior of L ( T ) is as in the disordered phase dis-cussed above: L ( T ) = O ( T c − T ). If L ( T ) has slower convergence than L ( T ), the convergence of L ( T ) will bedominated by that of L ( T ).In the HMF model, we can construct the angle-actionvariables explicitly, and the estimation of L ( T ) is ratherstraightforward. In our spin model, such an explicit con-struction does not seem feasible owing to the form of thesingle-spin Hamiltonian H , so that we have to invokesome physical observations and assumptions in order toestimate L ( T ). Details of the estimation are presentedin Appendix C, and one gets L ( T ) = O (( m ) / ( n +1) ) = O (( T c − T ) β/ ( n +1) ) . (53)We remark that this estimation does not depend on thechoice of the reference family G . From the above equa-tion, it follows that the critical exponent γ − is γ − = βn + 1 , (54)provided β ≤ n + 1, which is satisfied for β = 1 / n ≥
3. Nonlinear response at the critical point
As mentioned earlier, L ( T c ) = 0, and Eq. (43) gives toleading order in h the result N ( T c )( δm + h ) − h = 0 (55)at the critical point T c . The first term of N ( T ), seeEq. (45), vanishes on using the fact that m = 0 at T = T c gives (cid:104) S x (cid:105) = 0. As a result, N ( T c ) becomes N ( T c ) = − (cid:90) (cid:90) µ (cid:104) S x (cid:105) h F (cid:48) ( H h ( φ, p ))d φ d p, (56)a form that reduces to the one for L ( T ), Eq. (52), onreplacing in the latter the reference state f with theasymptotic state f h in performing the average over S x .We thus have an estimation of N ( T c ) as N ( T c ) = O (( δm + h ) / ( n +1) ) , (57)where we have used Eq. (53) and have replaced m in itwith m h + h = δm + h . This estimation gives( δm + h ) ( n +2) / ( n +1) ∝ h (58)and hence, that δm ∝ h ( n +1) / ( n +2) . (59)The critical exponent δ is thus δ = n + 2 n + 1 . (60) Model- n β γ + γ − δ Statistical Mechanics n ≥ / n = 0 1 / n ≥ βn + 1 n + 2 n + 1TABLE I. Critical exponents of the spin model (2) obtainedwithin the Vlasov dynamics. The critical exponents γ ± and δ in Model-0 reflect no response within the Vlasov dynamics.The scaling relation γ − = β ( δ −
1) holds for all cases.
4. Predicted critical exponents and the scaling relation
The theoretically predicted critical exponents, ob-tained within the Vlasov dynamics, are thus β = 12 , γ + = 1 , γ − = βn + 1 , δ = n + 2 n + 1 , ( n ≥ . (61)These exponents satisfy the scaling relation γ − = β ( δ − , (62)irrespective of the value of β .We remark that Model-0 corresponds to the limit n →∞ , since this limit eliminates from the Hamiltonian theanisotropic term Dp n for | p | <
1. In this limit, we have γ − = 0 and δ = 1, which is consistent with (37) and noresponse in Model-0. The obtained critical exponents aredisplayed in Table I. V. NUMERICAL TESTS
In this section, we discuss numerical checks of our the-oretical predictions for the critical exponents obtained inthe preceding section. As a representative case, we fo-cus on thermal equilibrium states as the reference states,which are represented by Eq. (24) and give β = 1 /
2. We,however, underline that the theoretical results developedin Sec. IV hold for other families of reference states.The numerical simulations of the equations of motion(7) are performed by using a fourth-order Runge-Kuttaalgorithm with the timestep δt = 0 .
01. We choose suffi-ciently large numbers of spins, namely N = 10 or 10 .We will refer to T as the temperature, but it is just a pa-rameter characterizing the reference state (24) and thereis no thermal noise in the dynamics. A. Temporal evolution of magnetization
Considering n = 1, and preparing the system in thethermal equilibrium state at a temperature T > T c , weshow in Fig. 1 the behavior of the magnetization m x asa function of time when a constant field of strength h =0 .
01 along the x -axis is turned on at t = 10. In thefigure, we also show by the dashed line the value of the magnetization induced by the field and given by Eq. (46).We can numerically examine the critical exponents γ ± and δ by varying the temperature T and observing theresponse that corresponds to the difference between thezero level and the dashed-line level of the magnetizationin the figure.The external field may have induced periodic oscilla-tions in the magnetization [27], but in our case, one mayobserve from Fig. 1 that the magnetization does not ex-hibit stable oscillations after the field is turned on. Itrather exhibits due to finiteness of the number of spinsonly fluctuations about the theoretically predicted valuevalid in the thermodynamic limit, as has also been ob-served in the HMF model [15]. t − . − . . . . . . m x FIG. 1.
Model- : Considering D = 5 and N = 10 , thefigure shows the temporal evolution of m x . The initial stateis the thermal equilibrium state at temperature T ( > T c ≈ . .
8. A constant field of strength h = 0 .
01 alongthe x -axis is turned on at time t = 10. The dashed linegives the value of the magnetization induced by the field andgiven by Eq. (46) to be ≈ . B. Critical exponents
Turning on a constant external field at t = 10, westudy numerically the response of the system (2) with N = 10 . Our results, presented in Figs. 2 and 3 forModel-0, Figs. 4 – 6 for Model-1, and in Figs. 7 – 9 forModel-2 are all consistent with our theoretical predic-tions in Table I. For Model-1, we take D = 5 for which T c ≈ .
476 is obtained by using G ( x ) = exp( − x/T ) inthe expression given by Eq. (A6) for the quantity A ( T )and then solving A ( T c ) = 0, while for Model-2, we take D = 15 for which one has T c ≈ .
47. Note that in Figs. 5and 8, our theoretical results match with our numericalresults only for sufficiently small h , as expected on thebasis of the fact that our theoretical analysis is valid inthe linear response regime obtained in the limit h → T c , numerical results forfinite N shown in Figs. 5 and 8 do not show the diver-gence predicted by our theory and shown in these figuresby red lines, owing to finiteness of the field strength h (the theoretical results are valid in the limit h → N → ∞ while satisfying the condition h > / √ N thatensures that the response dominates over finite-size fluc-tuations; in these figures, N is large enough that thecondition h > / √ N is satisfied, although not the limit h → h → − − − T c − T . . . . . m FIG. 2.
Model- : Spontaneous magnetization m (points),obtained by solving the self-consistent equation (27) withthermal equilibrium as the reference state. The critical pointis T c = 1 /
3. The line corresponds to the behavior (30), with β given by our theoretical analysis as β = 1 /
2, see Table I.
VI. CONCLUSIONS
In this work, we have discussed response to an externalfield in mean-field systems of classical Heisenberg spinsexhibiting a second-order phase transition in the station-ary state. The time evolution in the thermodynamic limitof such systems is described by the so-called Vlasov equa-tion for the single-spin phase space distribution function.We have shown that for Vlasov-stationary states that al-low a second-order phase transition and when subject toa small external field, the critical exponents characteriz-ing power-law behavior of response close to the criticalpoint and obtained within the Vlasov dynamics may takevalues different from the ones obtained on the basis of astatistical mechanical analysis, with no reference to thedynamics of the initial state in the presence of the exter-nal field. Interestingly, we find that both the sets of val-ues of the critical exponents satisfy the same scaling rela-tion, which is incidentally the same as the one known formean-field Hamiltonian particle systems that are quitedifferent from the studied spin systems. This work hintson one hand at the universality of critical behavior formean-field systems evolving under Vlasov dynamics, andcautions on the other hand against relying on a static ap-proach, with no reference to the dynamical evolution, toextract critical exponent values for mean-field systems. − − − h − × − × − × − × − m h T = 0 . . . − − − h − − m h T = 0 . . . FIG. 3.
Model- : For
T > T c = 1 / T < T c (lower panel), the figure shows the magnetization m h obtained by numerically integrating the equations of mo-tion (36) with N = 10 and performing a time average of theinstantaneous magnetization over an interval of length 20,which is further averaged over 5 realizations of the dynamics.The magnetization m h is independent of h in both cases, thuslending support to the behavior (32), with γ ± given by ourtheoretical analysis as γ + = γ − = 0, see Table I. − − − T c − T . . . . . . m FIG. 4.
Model- : For D = 5, the figure shows the spon-taneous magnetization m (points), obtained by solving theself-consistent equation (27) with thermal equilibrium as thereference state. The critical point is T c ≈ . β given by our theoreticalanalysis as β = 1 /
2, see Table I.
The reason that one has to resort to the dynamics in − − − T − T c . . . . . . . χ h = 0 . . . . . − − − − − T c − T χ h = 0 . . . . . FIG. 5.
Model- : For D = 5, the figure shows the sus-ceptibility χ ( T ) denoted by points; here, the magnetization m h is obtained by numerically integrating the equations ofmotion (7) with N = 10 and performing a time average ofthe instantaneous magnetization over an interval of length20, which is further averaged over 2 realizations of the dy-namics. The upper panel corresponds to the disordered phase T > T c ≈ . T < T c . In either case, the black dashed line corresponds tothe behavior close to T c , Eq. (32), with γ + = 1 and γ − = 1 / L given by Eq. (44) computed numerically by using themethod detailed in Appendix D. order to extract the correct critical exponents is the fol-lowing. Let us first recall the protocol we employ in ex-tracting the critical exponents. We prepare the system inVlasov-stationary states parametrized by a parameter T ,and which allow a second-order phase transition and con-sequently a critical point T c . An example of such states isthe thermal equilibrium state for which the parameter T is the temperature. We then subject the system to a con-stant external field. As mentioned in the introduction,mean-field systems like ours when considered in the ther-modynamic limit remain trapped in Vlasov-stationarystates forever in time. The dynamics in such states isnon-ergodic due to existence of the Casimir invariants,so that one may not apply tools from Boltzmann-Gibbsequilibrium statistical mechanics to extract the criticalexponents characterizing the response of such states to − − h − m h FIG. 6.
Model- : For D = 5, the figure shows the nonlinearresponse at the critical point T c ≈ . m h , denoted by points, is obtained by numerically integratingthe equations of motion (7) with N = 10 and performing atime average of the instantaneous magnetization over an inter-val of length 20, which is further averaged over 2 realizationsof the dynamics. The black line corresponds to the behav-ior (33), with δ given by our theoretical analysis as δ = 3 / h does our theorymatch with numerical results. − − T c − T . . . . . m FIG. 7.
Model- : For D = 15, the figure shows the spon-taneous magnetization m (points), obtained by solving theself-consistent equation (27) with thermal equilibrium as thereference state. The critical point is T c ≈ .
47. The line cor-responds to the behavior (30), with β = 1 /
2, as predicted byour theory, see Table I. the external field, for the simple reason that ergodicitylies at the heart of the very foundation of equilibriumstatistical mechanics.It may be noted that one could very well have a classof Vlasov-stationary states that for our model do notallow a second-order but a first-order phase transition.One may mention the case of the HMF model where theVlasov-stationary thermal equilibrium state allows for asecond-order phase transition, but there are other classesof Vlasov-stationary states that allow a first-order phasetransition [11–14]. It would be interesting in the con-text of our model to study the response for other classesof Vlasov-stationary reference states that allow a first-0 − − − T − T c . . . . . . χ h = 0 . . . . . − − − − − − T c − T − χ h = 0 . . . . . FIG. 8.
Model- : For D = 15, the figure shows the sus-ceptibility χ ( T ), denoted by points; here, the field-inducedmagnetization m h is obtained by numerically integrating theequations of motion (7) with N = 10 and performing a timeaverage of the instantaneous magnetization over an intervalof length 20, which is further averaged over 2 realizations ofthe dynamics. The upper panel corresponds to the disorderedphase T > T c ≈ .
47, while the lower panel is for the orderedphase
T < T c . In either case, the black dashed line corre-sponds to the behavior close to T c , Eq. (32), with γ + = 1 and γ − = 1 / L given by Eq. (44) computed numerically byusing the method given in Appendix D. order phase transition and differ from the ones studiedhere. Another immediate follow up of this work would beto investigate the validity, in the context of the studiedspin model, of the Kubo fluctuation-dissipation theoremvalid for short-range systems prepared in thermal equilib-rium and subject to small external fields. Studies in thisdirection are underway and comparison with the HMFmodel results [29] will be reported elsewhere. Collective1 /f fluctuation due to the Casimir invariants is also aninteresting topic to pursue [30]. ACKNOWLEDGMENTS
The work of D.D. is supported by UGC-NET ResearchFellowship Sr. No. 2121450744, Ref. No. 21/12/2014(ii) − − h − − m h FIG. 9.
Model- : For D = 15, the figure shows the nonlinearresponse at the critical point T c ≈ .
47. The magnetization m h , denoted by red points, is obtained by numerically inte-grating the equations of motion (7) with N = 10 and per-forming a time average of the instantaneous magnetizationover an interval of length 20, which is further averaged over2 realizations of the dynamics. The black line corresponds tothe behavior (33), with δ given by our theoretical analysis as δ = 4 /
3, see Table I. As expected, only for small h does ourtheory match with numerical results. EU-V. Y.Y.Y. acknowledges the supports of JSPS KAK-ENHI Grant No. 16K05472. The manuscript was fi-nalized while S.G. was visiting the Quantitative LifeSciences section of the International Centre for Theo-retical Physics (ICTP), Trieste, and he would like toacknowledge the support and hospitality of the ICTP.S.G. acknowledges support from the Science and En-gineering Research Board (SERB), India, Grant No.TAR/2018/000023.
Appendix A: Derivation of the criticalexponents (34)
Here we discuss how one may obtain the values of thecritical exponents given in Eq. (34). The starting pointis the equivalent of Eqs. (20), (21) and (23) in presenceof a constant field h = ( h, , f h ( φ, p ) = F ( H h ( φ, p )) = G ( H h ( φ, p )) (cid:82)(cid:82) µ G ( H h ( φ, p ))d φ d p , (A1)with H h ( φ, p ) = Dp n − ( m h + h ) · S , (A2)and m h = (cid:90) (cid:90) µ S f h ( φ, p )d φ d p. (A3)The last equation gives m h = (cid:82)(cid:82) µ S x G ( Dp n − ( m h + h ) S x )d φ d p (cid:82)(cid:82) µ G ( Dp n − ( m h + h ) S x )d φ d p , (A4)1with S x = (cid:112) − p cos φ . We assume G to be asmooth function of its argument. Since we take G tobe parametrized by the parameter T , this means that G is also smooth with respect to T . Let us expand G in aTaylor series around Dp n . Substituting the Taylor seriesand then integrating over φ , we see that in the numeratoron the right hand side of Eq. (A4), only odd order termsin ( m h + h ) survive; In the denominator, on the contrary,only even order terms survive. Consequently, the righthand side of Eq. (A4) has only odd order terms of m h ,and the equation gives A ( T )( m h + h ) + B ( T )( m h + h ) − h = 0 , (A5)with A ( T ) ≡ (cid:82) − (1 − p ) G (cid:48) ( Dp n )d p (cid:82) − G ( Dp n )d p , (A6)and we have neglected the higher-order terms in ( m h + h )in obtaining Eq. (A5). Consider now the latter for h = 0. Assuming B ( T ) >
0, one has a non-zero so-lution for m for T < T c and only a zero solution for T > T c , where the critical point T = T c is where we have A ( T c ) = 0, while A ( T ) is positive (respectively, negative)for T > T c (respectively, T < T c ). The nonzero solu-tion, giving the spontaneous magnetization for T < T c ,is m = (cid:112) − A/B . The smoothness of G with respect tothe parameter T allows A ( T ) to be expanded in a Taylorseries around T c , giving A ( T ) ∝ ( T − T c ) close to T c , andthis gives β = 1 / h , and we have χ ( T ) = d m h d h (cid:12)(cid:12)(cid:12)(cid:12) h =0 = (cid:40) (1 − A ) /A ( T > T c ) , (1 + 2 A ) / ( − A ) ( T < T c ) . (A7)From the behavior A ( T ) ∝ ( T − T c ) around T = T c , weget γ ± = 1.At the critical point, the self-consistent equation (A5)reduces to B ( T c )( m h + h ) = h. (A8)The response m h ∝ h /δ ( δ >
1) is larger than h for small h , so that we have m h ∝ h , implying δ = 3. Appendix B: Derivation of the response formula,Eq. (38)
In this appendix, we summarize the derivation of theresponse formula (38) by following Ref. [28]. Notingthat the Vlasov equation is governed by the single-spinHamiltonian H , which depends on f through the mag-netization m , the idea is to expand the Hamiltonian H as H = H h + K, (B1) where H h defined in (28) is the asymptotic part, charac-terizing the stationary ( t → ∞ ) state f h , while K is thetransient part. For our spin model, the explicit form ofthe transient part is K ( φ, p, t ) = − m T ( t ) (cid:112) − p cos φ, (B2)where the transient magnetization m T ( t ) is obtained as m T ( t ) = (cid:90) (cid:90) µ (cid:112) − p cos φ g ( φ, p )d φ d p, (B3)and the transient state g is defined by g = f − f h . Thetransient quantities, g, K and m T ( t ), are not known apriori, but they do not appear in the final result of theresponse formula.Let us write the Vlasov equation (14) as ∂f∂t = L H f = L H h f + L K f, (B4)where the linear operator L H is defined as L H f ≡ ∂H∂φ ∂f∂p − ∂H∂p ∂f∂φ . (B5)Eq. (B4) is still exact. Now, we assume that contribu-tion from the transient part, L K f , is negligible, whichis justified under some assumptions for Hamiltonian sys-tems and may be related to the phenomenon of Landaudamping, see Ref. [28] for details.Under the aforementioned assumption, the formal so-lution to the Vlasov equation (B4) is f ( φ, p, t ) = exp[ t L h ] f ( φ, p ) , (B6)which represents temporal evolution of f under Hamil-tonian flow associated with the asymptotic Hamiltonian H h . Assuming ergodicity, a formula that replaces thetime average with a partial phase-space average with re-spect to a iso- H h surface giveslim t →∞ t (cid:90) t e s L H f ( φ, p )d s = (cid:104) f (cid:105) h . (B7)Noting that the left hand side is nothing but the asymp-totic stationary state f h , we obtain the response formula(38). Appendix C: Derivation of the estimation (53) forthe quantity L Here, we derive Eq. (53). For simplicity of notation,we use the same symbols ( w, I ) for angle-action variablesassociated with the single-spin Hamiltonian H as theones used for H h in Sec. IV B, but the latter do notappear in this section and no confusion should arise. Herewe consider positive n .2We start with Eq. (52). Noting that (cid:104) S x (cid:105) and H depend on only the action variable I , we rewrite L as L = (cid:90) (cid:90) µ ψ ( I )2 π d w d I = (cid:90) ψ ( I )d I, (C1)where ψ ( I ) = − π (cid:104) S x (cid:105) F (cid:48) ( H ( I )) , (C2)and we have used the fact that canonical transforma-tion from ( φ, p ) to ( w, I ) gives d φ d p = d w d I . Thesingle-spin Hamiltonian H has a separatrix which con-sists of the stable and unstable manifolds of the fixedpoint ( S x , S y , S z ) = ( − , ,
0) and encloses the point O = (1 , ,
0) on the phase space of the unit sphere. Onthe basis of this observation, we make an essential as-sumption for estimating L , namely, that the main con-tribution to L comes from the region around the point O . We now change twice the variables of integration in L . First, to divide the phase space into the inside andthe outside of the separatrix, the integration variable I is changed to energy E as L = (cid:90) E sep E min ψ ( I )Ω( I ) d E + (cid:90) E max E sep ψ ( I )Ω( I ) d E, (C3)where the frequency Ω( I ) is defined byΩ( I ) ≡ d H d I ( I ) . (C4)Let E min , E sep and E max denote respectively the mini-mum energy, the separatrix energy, and the maximumenergy. Following the essential assumption, we omit thesecond term of L in (C3). Second, to eliminate the de-pendence on m of the integration interval, E is changedto a variable k defined as k ≡ E − E min E sep − E min . (C5)Consequently, one has L (cid:39) (cid:90) ψ ( I ( k )) E sep − E min Ω( I ( k )) d k = 2 m (cid:90) ψ ( k )Ω( k ) d k, (C6)where ψ ( I ( k )) is simply denoted as ψ ( k ) for instance.To estimate Ω around the point O , which correspondsto ( φ, p ) = (0 , H as H ( φ, p ) (cid:39) Dp n + m φ . (C7)The action variable is the area enclosed by a periodicorbit, and hence, we have I = 12 π (cid:73) φ d p = 2 π (cid:114) m (cid:90) p max (cid:112) E − Dp n d p, (C8)with p max = ( E/D ) / (2 n ) . In terms of p ≡ (cid:18) ED (cid:19) / (2 n ) u, (C9) we have the action variable as I = I E ( n +1) / n √ m , I = 2 √ πD / (2 n ) (cid:90) (cid:112) − u n d u, (C10)which gives E = m n/ ( n +1)0 (cid:18) II (cid:19) n/ ( n +1) . (C11)The frequency Ω is thereforeΩ = m n/ ( n +1)0 (cid:101) Ω , (cid:101) Ω = 2 nn + 1 1 I (cid:18) II (cid:19) ( n − / ( n +1) . (C12)Putting all together, we have the estimation of L as L (cid:39) m / ( n +1)0 (cid:90) ψ ( k ) (cid:101) Ω( k ) d k. (C13)We remark that (cid:104) S x (cid:105) is zero in the disordered phase,but it does not vanish in the ordered phase even when thelimit m → H line is confinedto the direction φ around the point O corresponding to( φ, p ) = (0 , Appendix D: A method to compute averages (cid:104)·(cid:105) over angles In this appendix, we discuss a method to compute theangle average (cid:104)·(cid:105) . The single-spin Hamiltonian H hasthe angle-action variables ( w, I ) whose temporal evolu-tion is w ( t ) = w (0) + Ω( I ) t, I ( t ) = I (0) , (D1)where Ω( I ) = d H / d I . Changing the variable from w to t , we have the average as (cid:104) B (cid:105) = (cid:82) π B ( w, I )d w (cid:82) π d w = (cid:82) T p B ( w ( t ) , I )Ω( I )d t (cid:82) T p Ω( I )d t = (cid:82) T p B ( w ( t ) , I )d t (cid:82) T p d t , (D2)where T p is the period on the considered iso- I line. Wemay write down the canonical equations of motion for H = Dp n − m (cid:112) − p cos φ asd φ d t = 2 nDp n − + m p (cid:112) − p cos φ, d p d t = − m (cid:112) − p sin φ, (D3)and the average (cid:104) (cid:112) − p cos φ (cid:105) is computed as (cid:104) (cid:112) − p cos φ (cid:105) = (cid:82) T p (cid:112) − p ( t ) cos φ ( t ) d t (cid:82) T p d t . (D4)3Note that the left hand side depends on the action I only, with the initial condition ( φ (0) , p (0)) determiningthe value of I . [1] M. E. Fisher, The theory of equilibrium critical phenom-ena, Rep. Prog. Phys. , 615 (1967).[2] H. E. Stanley, Introduction to Phase Transitions andCritical Phenomena (Oxford University Press, UK,1987).[3] H. Nishimori and G. Ortiz,
Elements of Phase Transi-tions and Critical Phenomena (Oxford University Press,Oxford, 2011).[4] A. Campa, T. Dauxois and S. Ruffo, Statistical mechan-ics and dynamics of solvable models with long-range in-teractions, Phys. Rep. , 57 (2009).[5] Y. Levin, R. Pakter, F. B. Rizzato, T. N. Teles and F.P. C. Benetti, Nonequilibrium statistical mechanics ofsystems with long-range interactions, Phys. Rep. , 1(2014).[6] A. Campa, T. Dauxois, D. Fanelli and S. Ruffo,
Physicsof Long-range Interacting Systems (Oxford UniversityPress, Oxford, 2014).[7] S. Gupta and S. Ruffo, The world of long-range interac-tions: A birds eye view, Int. J. Mod. Phys. A , 1741018(2017).[8] Y. Y. Yamaguchi, J. Barr´e, F. Bouchet, T. Dauxois andS. Ruffo, Stability criteria of the Vlasov equation andquasi-stationary states of the HMF model, Physica A , 36 (2004).[9] J. Binney and S. Tremaine, Galactic Dynamics, 2nd ed. (Princeton University Press, Princeton, NJ, 2008).[10] J. Barr´e, F. Bouchet, T. Dauxois, S. Ruffo, and Y. Y.Yamaguchi, The Vlasov equation and the Hamiltonianmean-field model, Physica A , 177 (2006).[11] A. Antoniazzi, D. Fanelli, S. Ruffo, and Y. Y. Yamaguchi,Nonequilibrium Tricritical Point in a System with Long-Range Interactions, Phys. Rev. Lett. , 040601 (2007).[12] R. Pakter and Y. Levin, Core halo distribution in theHamiltonian mean-field model, Phys. Rev. Lett. ,200603 (2011).[13] T. M. Rocha Filho, M. A. Amato, and A. Figueiredo,Nonequilibrium phase transitions and violent relaxationin the Hamiltonian mean-field model, Phys. Rev. E ,062103 (2012).[14] T. N. Teles, F. P. da C. Benetti, R. Pakter, and Y. Levin,Nonequilibrium Phase Transitions in Systems with Long-Range Interactions, Phys. Rev. Lett. , 230601 (2012).[15] S. Ogawa, A. Patelli and Y. Y. Yamaguchi, Non-mean-field critical exponent in a mean-field model: Dynamicsversus statistical mechanics, Phys. Rev. E , 032131 (2014).[16] S. Inagaki and T. Konishi, Dynamical Stability of a Sim-ple Model Similar to Self-Gravitating Systems, Publ. As-tron. Soc. Japan , 733 (1993).[17] M. Antoni and S. Ruffo, Clustering and relaxation inHamiltonian long-range dynamics, Phys. Rev. E , 2361(1995).[18] S. Ogawa and Y. Y. Yamaguchi, Landau-like theory foruniversality of critical exponents in quasistationary statesof isolated mean-field systems, Phys. Rev. E , 062108(2015).[19] P. Mazur, Non-ergodicity of phase functions in certainsystems, Physica (Amsterdam) , 533 (1969).[20] M. Suzuki, Ergodicity, constants of motion, and boundsfor susceptibilities, Physica (Amsterdam) , 277 (1971).[21] M. Kac, G. E. Uhlenbeck and P. C. Hemmer, On thevan der Waals Theory of the Vapor-Liquid Equilibrium.I. Discussion of a one-dimensional model, J. Math. Phys. , 216 (1963).[22] S. Gupta and D. Mukamel, Quasistationarity in a modelof classical spins with long-range interactions, J. Stat.Mech.: Theory Exp., P03015 (2011).[23] J. Barr´e and S. Gupta, Classical Heisenberg spins withlong-range interactions: relaxation to equilibrium for fi-nite systems, J. Stat. Mech.: Theory Exp. P02017 (2014).[24] G. F. Kventsel and J. Katriel, Static properties ofinfinite-range spin Hamiltonians, Phys. Rev. B , 2828(1984).[25] S. Ogawa and Y. Y. Yamaguchi, Linear response theoryin the Vlasov equation for homogeneous and for inhomo-geneous quasistationary states, Phys. Rev. E , 061115(2012).[26] In this work, we measure temperature in units of theBoltzmann constant.[27] R. Pakter and Y. Levin, Nonequilibrium dynamics of aninfinite range XY model in an external field, J. Stat.Phys. , 531 (2013).[28] S. Ogawa and Y. Y. Yamaguchi, Nonlinear response forexternal field and perturbation in the Vlasov system,Phys. Rev. E , 052114 (2014).[29] Y. Y. Yamaguchi, Strange scaling and relaxation of finite-size fluctuation in thermal equilibrium, Phys. Rev. E ,012133 (2016).[30] Y. Y. Yamaguchi and K. Kaneko, Collective 1 /f fluc-tuation by pseudo-Casimir-invariants, Phys. Rev. E98