Critical exponents in mean-field classical spin systems

TL;DR

This paper investigates the critical exponents in mean-field classical spin systems by analyzing their dynamical evolution via the Vlasov equation, revealing significant differences from static analysis and emphasizing the importance of dynamics in understanding phase transitions.

Contribution

The study demonstrates that critical exponents derived from dynamical phase space evolution differ markedly from static analysis, highlighting the necessity of considering dynamics in mean-field systems.

Findings

Dynamical analysis yields different critical exponents than static methods.

Static approaches may be insufficient for accurate critical exponent determination.

Emphasizes the importance of phase space dynamics in phase transition analysis.

Abstract

For a mean-field classical spin system exhibiting a second-order phase transition in the stationary state, we obtain within the corresponding phase space evolution according to the Vlasov equation the values of the critical exponents describing power-law behavior of response to a small external field. The exponent values so obtained significantly differ from the ones obtained on the basis of an analysis of the static phase-space distribution, with no reference to dynamics. This work serves as an illustration that cautions against relying on a static approach, with no reference to the dynamical evolution, to extract critical exponent values for mean-field systems.

Tables

Table 1. Table 1: Critical exponents of the spin model ( 2 ) obtained within the Vlasov dynamics. The critical exponents γ ± superscript 𝛾 plus-or-minus \gamma^{\pm} and δ 𝛿 \delta in Model- 0 0 reflect no response within the Vlasov dynamics. The scaling relation γ − = β ​ ( δ − 1 ) superscript 𝛾 𝛽 𝛿 1 \gamma^{-}=\beta(\delta-1) holds for all cases.

	Model-$n$	$\beta$	${\gamma }^{+}$	${\gamma }^{-}$	$\delta$
Statistical Mechanics	$n\ge 0$	$1/2$	$1$	$1$	$3$
Vlasov dynamics	$n=0$	$1/2$	$0$	$0$	$1$
	$n\ge 1$	$\frac{1}{2}$	$1$	$\frac{\beta }{n+1}$	$\frac{n+2}{n+1}$

Equations

$${\mathbf{S}}_i=(S_{ix},S_{iy},S_{iz});i=1,2,\dots ,N.$$

$$H_N=-\frac{J}{2N}\sum _{i,j=1}^N{\mathbf{S}}_i\cdot {\mathbf{S}}_j+D\sum _{i=1}^N{S}_{iz}^{2n}-\mathbf{h}(t)\cdot \sum _{i=1}^N{\mathbf{S}}_i.$$

$$\mathbf{h}(t)=\Theta (t-t_0)\mathbf{h},$$

$${\dot{\mathbf{S}}}_i=\{{\mathbf{S}}_i,H_N\};i=1,2,\dots ,N,$$

$$\{XY,Z\}=\{X,Z\}Y+X\{Y,Z\}$$

$$\begin{array}{rl}& \{S_{ix},S_{jy}\}={\delta }_{ij}{S}_{iz},\\ & \{S_{iy},S_{jz}\}={\delta }_{ij}{S}_{ix},\\ & \{S_{iz},S_{jx}\}={\delta }_{ij}{S}_{iy}.\end{array}$$

$$\begin{array}{rl}{\dot{S}}_{ix}& =S_{iy}(m_z+h_z)-S_{iz}(m_y+h_y)-2nD{S}_{iy}{S}_{iz}^{2n-1},\\ {\dot{S}}_{iy}& =S_{iz}(m_x+h_x)-S_{ix}(m_z+h_z)+2nD{S}_{ix}{S}_{iz}^{2n-1},\\ {\dot{S}}_{iz}& =S_{ix}(m_y+h_y)-S_{iy}(m_x+h_x),\end{array}$$

$$\mathbf{m}\equiv \frac{1}{N}\sum _{i=1}^N{\mathbf{S}}_i=(m_x,m_y,m_z)$$

$$m_x{h}_y-m_y{h}_x=0$$

$$S_{ix}=\sin {\theta }_i\cos {\varphi }_i,S_{iy}=\sin {\theta }_i\sin {\varphi }_i,S_{iz}=\cos {\theta }_i,$$

$$\begin{array}{rl}{\dot{\theta }}_i& =(m_x+h_x)\sin {\varphi }_i-(m_y+h_y)\cos {\varphi }_i,\\ {\dot{\varphi }}_i& =(m_x+h_x)\cot {\theta }_i\cos {\varphi }_i+(m_y+h_y)\cot {\theta }_i\sin {\varphi }_i\\ & -(m_z+h_z)+2nD{\cos }^{2n-1}{\theta }_i.\end{array}$$

$$S_{ix}=\sqrt{1-p_i^2}\cos {\varphi }_i,S_{iy}=\sqrt{1-p_i^2}\sin {\varphi }_i,S_{iz}=p_i.$$

$$\{X,Y\}=\sum _{i=1}^N\left(\frac{\partial X}{\partial {\varphi }_i}\frac{\partial Y}{\partial p_i}-\frac{\partial X}{\partial p_i}\frac{\partial Y}{\partial {\varphi }_i}\right).$$

$$\frac{\partial f}{\partial t}+\frac{\partial H}{\partial p}\frac{\partial f}{\partial \varphi }-\frac{\partial H}{\partial \varphi }\frac{\partial f}{\partial p}=0,$$

$$H(\varphi ,p,t)=D{p}^{2n}-[\mathbf{m}(t)+\mathbf{h}(t)]\cdot \mathbf{S},$$

$$\mathbf{S}\equiv \left(\sqrt{1-p^2}\cos \varphi ,\sqrt{1-p^2}\sin \varphi ,p\right),$$

$$\mathbf{m}(t)={\iint }_{\mu }\mathbf{S}f(\varphi ,p,t)\mathrm{d}\varphi \mathrm{d}p.$$

$${\iint }_{\mu }X(\varphi ,p)\mathrm{d}\varphi \mathrm{d}p\equiv {\int }_0^{2\pi }\mathrm{d}\varphi {\int }_{-1}^1\mathrm{d}pX(\varphi ,p).$$

$$C[f](t)={\iint }_{\mu }c(f)\mathrm{d}\varphi \mathrm{d}p$$

$$f_0(\varphi ,p)=F(H_0(\varphi ,p))=\frac{G(H_0(\varphi ,p))}{{\iint }_{\mu }G(H_0(\varphi ,p))\mathrm{d}\varphi \mathrm{d}p},$$

$$H_0(\varphi ,p)=D{p}^{2n}-{\mathbf{m}}_0\cdot \mathbf{S},$$

$${\mathbf{m}}_0=(m_{0x},m_{0y},m_{0z})$$

$${\mathbf{m}}_0={\iint }_{\mu }\mathbf{S}{f}_0(\varphi ,p)\mathrm{d}\varphi \mathrm{d}p.$$

$$G(x)=\exp (-x/T).$$

$$G(x)=\frac{1}{\exp [(x-a)/b]+1}.$$

$$H_0(\varphi ,p)=D{p}^{2n}-m_0\sqrt{1-p^2}\cos \varphi ,$$

$$m_0={\iint }_{\mu }\sqrt{1-p^2}\cos \varphi F(H_0(\varphi ,p))\mathrm{d}\varphi \mathrm{d}p.$$

$$H_h(\varphi ,p)=D{p}^{2n}-(m_h+h)\sqrt{1-p^2}\cos \varphi .$$

$$\delta m\equiv m_h-m_0.$$

$$m_0(T)\propto (T_c-T{)}^{\beta };T\to T_c^{-},$$

Full text

I Introduction

Since early days of statistical mechanics, studying phase transitions in physical systems has been a theme of active research in the field. Phase transitions can occur only in the thermodynamic limit. Second-order or continuous phase transitions are characterized by a power-law behavior of macroscopic quantities close to the critical point of transition. Such transitions in different systems may be broadly classified into universality classes identified by different values of critical exponents describing the power-law behavior. For example, for a ferromagnet exhibiting a second-order phase transition as a function of temperature $T$, the magnetization close to and on the lower side of the critical point $T_{c}$ has a power-law dependence on the separation $(T_{c}-T)$ from the critical point, with the corresponding exponent being $\beta$. On applying an external field, the magnetization increases as a function of the field strength, and in the limit of an infinitesimal field, a linear growth for $T\neq T_{c}$ implying a linear response determines the zero-field susceptibility $\chi$. The susceptibility diverges as a power law close to and on both sides of the critical point, with the corresponding exponents denoted by $\gamma^{+}$ and $\gamma^{-}$ on the disordered ($T>T_{c}$) and the magnetized ($T<T_{c}$) phase, respectively. At the critical point, the response becomes nonlinear, being characterized by the critical exponent $\delta$. These critical exponents are known to satisfy the scaling relation $\gamma^{\pm}=\beta(\delta-1)$ Fisher ; Stanley ; nishimori-ortiz-11 .

One representative class of systems exhibiting second-order phase transitions is that of mean-field systems. In thermal equilibrium of such systems, statistical mechanical predictions for the critical exponents, based on an analysis of the thermal equilibrium phase space distribution with no reference to dynamics, yield the values $\beta=1/2,~{}\gamma^{\pm}=1$, $\delta=3$ Stanley . However, owing to the mean-field nature of the time evolution, critical exponents obtained on the basis of dynamics may well have different values. Indeed, dynamics of a mean-field system in the thermodynamic limit is described by the so-called Vlasov equation that allows a vast number of stable stationary states, and thermal equilibrium is just one of them campa-dauxois-ruffo-09 ; levin-etal-14 ; Campa:2014 ; Gupta:2017 . This implies that once the system is in a stable stationary state other than thermal equilibrium, it would not relax to thermal equilibrium. A large but finite system remains trapped in so-called quasistationary states (QSSs) identified as stable stationary solutions of the Vlasov equation, with finite-size effects allowing a slow evolution of the QSSs towards thermal equilibrium over a timescale that diverges with the system size yamaguchi-etal-04 ; binney-tremaine-08 .

Existence of QSSs allows nonequilibrium phase transitions: a generic initial state undergoes a violent relaxation to relax to a QSS barre-etal-06 , and the nonequilibrium phase transition can for example be defined with respect to the value of the order parameter in the QSS. In a given system, these nonequilibrium phase transitions may not necessarily be continuous even when the equilibrium phase transition is continuous, and several discontinuous nonequilibrium phase transitions have been reported in the literature antoniazzi-07 ; Pakter:2011 ; rochafilho-amato-figueiredo-12 ; teles-benetti-pakter-levin-12 . In this article, we, however, focus on families of QSSs that exhibit continuous phase transitions and to which an external field is applied in order to investigate the values of the critical exponents characterizing the response.

The aforementioned trapping scenario holds even when an external field is applied to the system prepared in a thermal equilibrium state: With the field on, a finite system goes from the initial to a new thermal equilibrium state via intermediate QSSs, while a thermodynamic system remains trapped in a QSS and does not relax to thermal equilibrium ogawa-patelli-yamaguchi-14 . The latter fact requires that one invokes an alternative strategy of obtaining susceptibility that is based on the Vlasov dynamics when addressing the issue of response of mean-field systems in thermal equilibrium to an external field. The critical exponents $\gamma^{\pm}$ so obtained may not necessarily coincide with the ones computed within equilibrium statistical mechanics. Indeed, in the so-called Hamiltonian mean-field (HMF) model inagaki-konishi-93 ; antoni-ruffo-95 , a paradigmatic mean-field system exhibiting a second-order phase transition, the critical exponents obtained within the Vlasov dynamics have been shown to be $\gamma^{+}=1,~{}\gamma^{-1}=1/4$ for a family of stable stationary initial states with $\beta=1/2$ ogawa-patelli-yamaguchi-14 . Moreover, at the critical point, the Vlasov dynamics gives $\delta=3/2$. More generally, the critical exponents within the Vlasov dynamics have been obtained as $\gamma^{+}=2\beta,~{}\gamma^{-}=\beta/2,~{}\delta=3/2$ for a class of Hamiltonian particle systems including the HMF model ogawa-yamaguchi-15 . Interestingly, the critical exponents obtained within the two approaches satisfy the same scaling relation, namely, $\gamma^{-}=\beta(\delta-1)$.

The difference in the values of the critical exponents that are obtained based on statistical mechanics and dynamics stems from the existence of an infinite number of so-called Casimir invariants that are constants of motion for the Vlasov dynamics and which make the dynamics non-ergodic. In general, existence of constraints suppresses susceptibility mazur-69 ; suzuki-71 and accordingly the values of the critical exponents $\gamma^{\pm}$. This fact is borne out by the values of the critical exponents obtained in the HMF model by taking into account the Casimir constraints. Another important remark is that divergence of susceptibility is observed even when the dynamics is non-ergodic, as is found in the case of the HMF model ogawa-patelli-yamaguchi-14 .

The HMF model mimics the classical $XY$ model, with an additional kinetic energy term assigned to individual spins. Owing to the latter whose range is the whole real set, the one-particle phase space of the HMF model is a cylinder. In the HMF model, the Poisson bracket between the spin components is taken to vanish identically, In this work, we consider Heisenberg spin systems with mean-field interactions, in which the Poisson brackets between the spin components are strictly nonzero, and the single-particle phase space is the unit sphere. Considering the time evolution of the spin components according to a Hamiltonian with a mean-field interaction and a local anisotropy, we address here several questions of theoretical and practical relevance: Does the universality class for usual Hamiltonian systems defined on a cylinder, e.g., the HMF model, include spin systems defined on the unit sphere? What is the effect of the anisotropy on the critical exponents? Would the scaling relation $\gamma^{-}=\beta(\delta-1)$ still hold even if the spin system is found to be in a different universality class?

This paper is organized as follows. The spin model we study is introduced in Sec. II. Here, the dynamics described by the canonical equations of motion is also discussed, as is the characterization of the dynamics in the thermodynamic limit in terms of the Vlasov equation. Based on the latter, we discuss the setting and the definition of the critical exponents in Sec. III, while our theoretical predictions for the critical exponents are derived in Sec. IV. Detailed numerical checks of our theoretical predictions are pursued in Sec. V. Section VI concludes the paper with discussions.

II The model

II.1 Definition

Our model of study consists of $N$ globally-coupled classical Heisenberg spins of unit length denoted by

[TABLE]

The $N$-body Hamiltonian of the model is given by

[TABLE]

Here, the first term with $J>0$ on the right-hand side models a ferromagnetic mean-field interaction between the spins. The coupling constant $J$ has been scaled down by the system size $N$ in order to make the energy extensive, in accordance with the Kac prescription Kac:1963 . The system (2) is however intrinsically non-additive: it cannot be trivially subdivided into independent macroscopic parts. In the following, we set $J=1$ without loss of generality.

In Eq. (2), the second term with $D>0$ on the right-hand side accounts for local anisotropy; Restricting to the subclass of models that are symmetric under $S_{iz}\to-S_{iz}$, we have made here the choice of even exponent equal to $2n$, with $n$ being a non-negative integer. We refer to the model with exponent $2n$ as Model-$n$. Note that Model-[math] is completely isotropic in the spin space, and there is no preferred direction of orientation of spins. Model-$1$ has been studied previously in the context of QSSs in Refs. Gupta:2011 ; Barre:2014 . Model-$2$ is the special case of a quartic anisotropy; it may be noted that thermodynamic properties of a Heisenberg spin model containing a quartic term have been studied in Ref. kventsel-katriel-84 .

The third term on the right-hand side of Eq. (2) arises due to the application of a time-dependent external magnetic field $\boldsymbol{h}(t)\equiv(h_{x}(t),h_{y}(t),h_{z}(t))$. In this work, we consider the external field to be absent for times previous to instant $t_{0}$, i.e., for times $t<t_{0}$, when the system will be assumed to be existing in a reference state, e.g., a thermal equilibrium state. For times $t\geq t_{0}$, on the other hand, we would put on a constant field in order to measure the response of the reference state to the external field. The explicit form of $\boldsymbol{h}(t)$ is thus given by

[TABLE]

where $\Theta(t)$ is the unit step function, and $\boldsymbol{h}$ is a vector of constant length equal to $h$. The singularity of the unit step function $\Theta(t)$ will have no effect on the values of the critical exponents obtained based on the Vlasov dynamics, and we may replace $\Theta(t)$ with a smooth function ogawa-yamaguchi-12 .

II.2 Spin dynamics

In dimensionless times, the time evolution of system (2) is governed by the set of coupled first-order differential equations

[TABLE]

where the dot denotes derivative with respect to time. The Poisson bracket $\{\cdot,\cdot\}$ is bilinear, skew-symmetric, and satisfies the Leibniz’s rule

[TABLE]

for any functions $X,Y,$ and $Z$ of the spins. The Poisson brackets between two spins are given by

[TABLE]

Using Eqs. (2), (4), and (6), we obtain the time evolution of the spin components as

[TABLE]

where

[TABLE]

is the magnetization vector that serves as the mean field governing 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

[TABLE]

is satisfied. The length of each spin is a constant of motion, and so is the total energy of the system when the field $\boldsymbol{h}$ is time independent.

Writing the spin components in terms of spherical polar angles $\theta_{i}\in[0,\pi]$ and $\phi_{i}\in[0,2\pi)$, as

[TABLE]

we obtain from Eq. (7) the time evolution of the variables $\theta_{i}$ and $\phi_{i}$ as

[TABLE]

For later convenience, we introduce a new variable $p_{i}\equiv\cos\theta_{i}$, in terms of which we have

[TABLE]

In terms of $p_{i}$, which is in fact canonically conjugate to $\phi_{i}$, the Poisson bracket $\{\cdot,\cdot\}$ reads Gupta:2011

[TABLE]

The dynamical variables of the $i$-th spin are thus $\phi_{i}$ and $p_{i}$, while a volume element in the $(\phi_{i},p_{i})$-space is ${\rm d}\phi_{i}{\rm d}p_{i}$.

II.3 Description in the thermodynamic limit

In the thermodynamic limit $N\to\infty$, the dynamics of system (2) is described by the Vlasov equation

[TABLE]

where $f(\phi,p,t)$ is the single-spin distribution function that measures the probability density to find a spin $(\phi,p)$ at time $t$, while the single-spin Hamiltonian $H$ is

[TABLE]

with

[TABLE]

and the magnetization vector $\boldsymbol{m}=(m_{x},m_{y},m_{y})$ given by

[TABLE]

The double integral over any function $X(\phi,p)$ in the single-spin phase space $\mu\equiv(\phi,p)$ is defined as

[TABLE]

Note that the single-spin Hamiltonian (15) depends on time $t$ through the magnetization $\boldsymbol{m}(t)$ and the external field $\boldsymbol{h}(t)$. Normalization of $f(\phi,p,t)$ reads $\iint_{\mu}f(\phi,p,t)=1$ for any time $t$.

Any quantity

[TABLE]

is a constant of motion for any smooth function $c$, as may be seen by considering the time variation of $C$ and using Eq. (14). These invariants of motion are called Casimir invariants, which hold even when the single-spin Hamiltonian depends on time. The Casimir invariants do not allow an initial state with $\boldsymbol{h}=\boldsymbol{0}$ to relax to the thermal equilibrium state with $\boldsymbol{h}\neq\boldsymbol{0}$ when at least one of the Casimir invariants $C[f]$ between the two states is not the same.

III Setting and definition of the critical exponents

III.1 Setting

For $t<t_{0}$, we consider system (2) to be existing in one of a family of stable stationary states with external field $\boldsymbol{h}=\boldsymbol{0}$. In order that we may study the critical exponents associated with the response of the system to an external field that we put on for times $t\geq t_{0}$, we restrict to a family of states that allow a second-order phase transition and consequently a critical point in the stationary state. We refer to such a family of states as our reference states and denote the states by $f_{0}$. From Eq. (14), it is evident that $f_{0}$ of the form

[TABLE]

with $G$ an arbitrary function, is a stationary solution of the Vlasov equation, and we have

[TABLE]

and

[TABLE]

satisfying the self-consistent equation

[TABLE]

The family of functions $G$ may be parametrized by a parameter $T$, which in the case of thermal equilibrium coincides with the temperature note :

[TABLE]

However, the analysis presented in the following applies to other family of functions $G$, such as the Fermi-Dirac-type family

[TABLE]

In this case, the parameter $T$ may be identified with either of the two parameters $a$ and $b$.

Now, from the rotational symmetry of $H_{0}(\phi,p)$ on the $(S_{x},S_{y})$-plane, we may set $m_{0y}=0$ without loss of generality. Moreover, we may assume $m_{0z}=0$, which solves the self-consistent equation for $m_{0z}$. Denoting $m_{0x}$ by $m_{0}$, so that $\boldsymbol{m}_{0}=(m_{0},0,0)$, we have

[TABLE]

while the self-consistent equation (23) reads

[TABLE]

At $t=t_{0}$, we turn on a constant external field $\boldsymbol{h}=(h,0,0)$ pointing in the direction of the reference magnetization $\boldsymbol{m}_{0}=(m_{0},0,0)$. In presence of the external field, the system evolving under the Vlasov dynamics (14) relaxes from the reference state $f_{0}$ to a stationary state $f_{h}$ with magnetization $\boldsymbol{m}_{h}=(m_{h},0,0)$. The single-spin Hamiltonian corresponding to the state $f_{h}$ is

[TABLE]

We stress that $f_{h}$ is not necessarily the thermal equilibrium state proportional to $\exp[-H_{h}(\phi,p)/T]$, and thus could be an out-of-equilibrium state, see Sec. II.3. Within the Vlasov dynamics, the response to the external field is measured by

[TABLE]

In the above setting, we recall the definitions of the critical exponents $\beta,\gamma^{+},\gamma^{-}$ and $\delta$ given in any standard reference on critical phenomena, e.g., Ref. Fisher .

III.2 Definition of the critical exponents

The critical exponent $\beta$ is defined with respect to the reference state, as

[TABLE]

where $T_{c}$ is the critical point. Here, $m_{0}$ is the positive solution of the self-consistent equation (27), and the value of $\beta$ may depend on the choice of the family $F$ of the reference state. The self-consistent equation (27), however, implies quite generally that $\beta=1/2$, see Appendix A.

The critical exponents $\gamma^{\pm}$ are defined in the regime of linear response. The response $\delta m$ depends on $T$ and $h$, and the susceptibility $\chi(T)$ is defined as

[TABLE]

The susceptibility diverges at the critical point $T_{c}$ as

[TABLE]

which defines the exponents $\gamma^{\pm}$.

At the critical point $T_{c}$, one has

[TABLE]

which defines the critical exponent $\delta$. Usually, one has $\delta>1$, since the leading response is nonlinear and is stronger than the linear response.

III.3 Statistical mechanics predictions for the critical exponents

Statistical mechanics analysis that considers studying the equivalent of Eqs. (20), (21) and (23) in presence of a constant field $\boldsymbol{h}=(h,0,0)$, with no reference to dynamics, gives

[TABLE]

irrespective of the value of the exponent $n$, see Appendix A for details. In the next section, we derive the values of the critical exponents within the Vlasov dynamics. We will obtain the response $\delta m$ within the Vlasov dynamics, and hence the exponents $\gamma^{\pm}$ and $\delta$ may very well take values different from the ones in Eq. (34).

IV Theoretical predictions for the critical exponents based on the

Vlasov dynamics

In this section, we derive our results for the critical exponents based on the Vlasov dynamics. As already mentioned above, we have $\beta=1/2$ quite generally for all choices of the family $F$ of the reference state. In the following, we discuss the computation of the critical exponents $\gamma^{\pm}$ and $\delta$ for a given value of the exponent $\beta$.

IV.1 Model-[math]

The single-spin Hamiltonian of Model-[math] is

[TABLE]

The equations of motion are obtained from Eq. (7) as

[TABLE]

Clearly, $S_{ix}$ for any $i$ and consequently $m_{x}$ are constant of motion for any external field $h$ irrespective of its time dependence. Dynamically, each spin rotates on a $S_{x}=$ constant plane. The fact that the variable $S_{x}$ is a constant of motion for both cases of $h=0$ and $h\neq 0$ implies that the reference state $f_{0}=F(H_{0})=F(-m_{0}S_{x})$ is stationary even after the external field is turned on, and we have $m_{h}=m_{0}$. Consequently, no response to the external field is obtained within the Vlasov dynamics. If we have to assign values to the critical exponents, we may say

[TABLE]

owing to the fact that no divergence of the susceptibility is obtained within the Vlasov dynamics. The aforementioned exponent values are quite different from the ones obtained within statistical mechanics, Eq. (34). For the case $G(x)=\exp(-x/T)$, the critical point is $T_{\rm c}=1/3$ (obtained by using the results in Appendix A, in particular, by substituting such a form of $G(x)$ into the function $A(T)$ defined by Eq. (68) and then solving $A(T_{\rm c})=0$).

The isotropic spin model, Model-[math], shows no response to the external field and thus provides a simple and extreme example of dynamical suppression of response, but Model-$n$ with $n\geq 1$ does show non-zero response. In the following subsection, we obtain the values of the critical exponents for $n\geq 1$.

IV.2 Model-$n$ with $n\geq 1$

To compute the values of the critical exponents $\gamma^{\pm}$ and $\delta$, our task is to obtain within the Vlasov dynamics starting from the state $f_{0}$ the asymptotic state $f_{h}$ and hence the response $\delta m$. For this purpose, we remark that the Hamiltonian $H_{h}$ given in Eq. (28) is integrable and has the associated angle-action variables $(w,I)$. In this setting, the response formula

[TABLE]

has been proposed for Hamiltonian systems ogawa-yamaguchi-12 ; ogawa-yamaguchi-15 , where $\langle\cdot\rangle_{h}$ is defined as the average over the angle variable $w$, as

[TABLE]

In Appendix B, we summarize the derivation of Eq. (38). Using $H_{h}=H_{h}(I)$ and

[TABLE]

for any function $\varphi$, and the expansion of the single-spin Hamiltonian as

[TABLE]

with $S_{x}=\sqrt{1-p^{2}}\cos\phi$, we have the expansion of $f_{h}(I)$ as

[TABLE]

Note that, for instance, one has $F^{\prime}(x)=-F(x)/T$ for thermal equilibrium reference state (24). The average $\langle\cdot\rangle_{0}$ is defined as an average over the angle variable associated with the integrable system $H_{0}$.

Multiplying Eq. (42) by $S_{x}=\sqrt{1-p^{2}}\cos\phi$ and then integrating over $\phi$ and $p$, we have the self-consistent equation for the response $\delta m$ as

[TABLE]

where the coefficient $L$ of the linear part is

[TABLE]

while ${\cal N}$ concerns the leading nonlinear part:

[TABLE]

The linear part $L$ gives the values of the critical exponents $\gamma^{\pm}$, while the nonlinear part ${\cal N}$ gives the value of the exponent $\delta$.

Note that $L(T_{c})=0$, so that the contribution of only the nonlinear response appears at the critical point $T_{c}$. Away from the critical point, it is the linear part that gives the dominant contribution in Eq. (43), so that neglecting the nonlinear contribution, we have the linear response

[TABLE]

Then, within linear response, the divergence of the susceptibility,

[TABLE]

is determined by the convergence behavior of $L$ as $T\to T_{c}^{\pm}$.

IV.2.1 Linear response in the disordered phase

In the disordered phase, the angle variable $w$ is nothing but $\phi$, and we have

[TABLE]

This result implies that $L$ has no contribution from the dynamics, and hence may be expanded in a Taylor series in $(T-T_{\rm c})$ around $T_{c}$, resulting in its convergence being proportional to $T-T_{c}$. The critical exponent $\gamma^{+}$ is, therefore, given by

[TABLE]

IV.2.2 Linear response in the ordered phase

In the ordered phase, let us divide $L$ into the two parts:

[TABLE]

with

[TABLE]

and

[TABLE]

The behavior of $L_{1}(T)$ is as in the disordered phase discussed above: $L_{1}(T)=O(T_{c}-T)$. If $L_{2}(T)$ has slower convergence than $L_{1}(T)$, the convergence of $L(T)$ will be dominated by that of $L_{2}(T)$.

In the HMF model, we can construct the angle-action variables explicitly, and the estimation of $L_{2}(T)$ is rather straightforward. In our spin model, such an explicit construction does not seem feasible owing to the form of the single-spin Hamiltonian $H_{0}$, so that we have to invoke some physical observations and assumptions in order to estimate $L_{2}(T)$. Details of the estimation are presented in Appendix C, and one gets

[TABLE]

We remark that this estimation does not depend on the choice of the reference family $G$. From the above equation, it follows that the critical exponent $\gamma^{-}$ is

[TABLE]

provided $\beta\leq n+1$, which is satisfied for $\beta=1/2$ and $n\geq 1$.

IV.2.3 Nonlinear response at the critical point

As mentioned earlier, $L(T_{c})=0$, and Eq. (43) gives to leading order in $h$ the result

[TABLE]

at the critical point $T_{c}$. The first term of ${\cal N}(T)$, see Eq. (45), vanishes on using the fact that $m_{0}=0$ at $T=T_{c}$ gives $\langle S_{x}\rangle_{0}=0$. As a result, ${\cal N}(T_{c})$ becomes

[TABLE]

a form that reduces to the one for $L_{2}(T)$, Eq. (52), on replacing in the latter the reference state $f_{0}$ with the asymptotic state $f_{h}$ in performing the average over $S_{x}$. We thus have an estimation of ${\cal N}(T_{c})$ as

[TABLE]

where we have used Eq. (53) and have replaced $m_{0}$ in it with $m_{h}+h=\delta m+h$. This estimation gives

[TABLE]

and hence, that

[TABLE]

The critical exponent $\delta$ is thus

[TABLE]

IV.2.4 Predicted critical exponents and the scaling relation

The theoretically predicted critical exponents, obtained within the Vlasov dynamics, are thus

[TABLE]

These exponents satisfy the scaling relation

[TABLE]

irrespective of the value of $\beta$.

We remark that Model-[math] corresponds to the limit $n\to\infty$, since this limit eliminates from the Hamiltonian the anisotropic term $Dp^{2n}$ for $|p|<1$. In this limit, we have $\gamma^{-}=0$ and $\delta=1$, which is consistent with (37) and no response in Model-[math]. The obtained critical exponents are displayed in Table 1.

V Numerical tests

In this section, we discuss numerical checks of our theoretical predictions for the critical exponents obtained in the preceding section. As a representative case, we focus on thermal equilibrium states as the reference states, which are represented by Eq. (24) and give $\beta=1/2$. We, however, underline that the theoretical results developed in 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-Kutta algorithm with the timestep $\delta t=0.01$. We choose sufficiently large numbers of spins, namely $N=10^{6}$ or $10^{7}$. We will refer to $T$ as the temperature, but it is just a parameter characterizing the reference state (24) and there is no thermal noise in the dynamics.

V.1 Temporal evolution of magnetization

Considering $n=1$, and preparing the system in the thermal equilibrium state at a temperature $T>T_{\rm c}$, we show in Fig. 1 the behavior of the magnetization $m_{x}$ as a function of time when a constant field of strength $h=0.01$ along the $x$-axis is turned on at $t_{0}=10$. In the figure, 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 $\gamma^{\pm}$ and $\delta$ by varying the temperature $T$ and observing the response that corresponds to the difference between the zero level and the dashed-line level of the magnetization in the figure.

The external field may have induced periodic oscillations in the magnetization Pakter:2013 , but in our case, one may observe from Fig. 1 that the magnetization does not exhibit stable oscillations after the field is turned on. It rather exhibits due to finiteness of the number of spins only fluctuations about the theoretically predicted value valid in the thermodynamic limit, as has also been observed in the HMF model ogawa-patelli-yamaguchi-14 .

V.2 Critical exponents

Turning on a constant external field at $t=10$, we study numerically the response of the system (2) with $N=10^{7}$. Our results, presented in Figs. 2 and 3 for Model-[math], Figs. 4 – 6 for Model-$1$, and in Figs. 7 – 9 for Model-$2$ are all consistent with our theoretical predictions in Table 1. For Model-$1$, we take $D=5$ for which $T_{c}\approx 0.476$ is obtained by using $G(x)=\exp(-x/T)$ in the expression given by Eq. (68) 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}\approx 0.47$. Note that in Figs. 5 and 8, our theoretical results match with our numerical results only for sufficiently small $h$, as expected on the basis of the fact that our theoretical analysis is valid in the linear response regime obtained in the limit $h\to 0$. Let us remark that very close to $T_{c}$, numerical results for finite $N$ shown in Figs. 5 and 8 do not show the divergence predicted by our theory and shown in these figures by red lines, owing to finiteness of the field strength $h$ (the theoretical results are valid in the limit $h\to 0$ and $N\to\infty$ while satisfying the condition $h>1/\sqrt{N}$ that ensures that the response dominates over finite-size fluctuations; in these figures, $N$ is large enough that the condition $h>1/\sqrt{N}$ is satisfied, although not the limit $h\to 0$). Moreover, the convergence to the $h\to 0$ limit is slower for Model-$2$ than for Model-$1$.

VI Conclusions

In this work, we have discussed response to an external field in mean-field systems of classical Heisenberg spins exhibiting a second-order phase transition in the stationary state. The time evolution in the thermodynamic limit of such systems is described by the so-called Vlasov equation for the single-spin phase space distribution function. We have shown that for Vlasov-stationary states that allow a second-order phase transition and when subject to a small external field, the critical exponents characterizing power-law behavior of response close to the critical point and obtained within the Vlasov dynamics may take values different from the ones obtained on the basis of a statistical mechanical analysis, with no reference to the dynamics of the initial state in the presence of the external field. Interestingly, we find that both the sets of values of the critical exponents satisfy the same scaling relation, which is incidentally the same as the one known for mean-field Hamiltonian particle systems that are quite different from the studied spin systems. This work hints on one hand at the universality of critical behavior for mean-field systems evolving under Vlasov dynamics, and cautions on the other hand against relying on a static approach, with no reference to the dynamical evolution, to extract critical exponent values for mean-field systems.

The reason that one has to resort to the dynamics in order to extract the correct critical exponents is the following. Let us first recall the protocol we employ in extracting the critical exponents. We prepare the system in Vlasov-stationary states parametrized by a parameter $T$, and which allow a second-order phase transition and consequently a critical point $T_{\rm c}$. An example of such states is the thermal equilibrium state for which the parameter $T$ is the temperature. We then subject the system to a constant external field. As mentioned in the introduction, mean-field systems like ours when considered in the thermodynamic limit remain trapped in Vlasov-stationary states forever in time. The dynamics in such states is non-ergodic due to existence of the Casimir invariants, so that one may not apply tools from Boltzmann-Gibbs equilibrium statistical mechanics to extract the critical exponents characterizing the response of such states to the external field, for the simple reason that ergodicity lies at the heart of the very foundation of equilibrium statistical mechanics.

It may be noted that one could very well have a class of Vlasov-stationary states that for our model do not allow a second-order but a first-order phase transition. One may mention the case of the HMF model where the Vlasov-stationary thermal equilibrium state allows for a second-order phase transition, but there are other classes of Vlasov-stationary states that allow a first-order phase transition antoniazzi-07 ; Pakter:2011 ; rochafilho-amato-figueiredo-12 ; teles-benetti-pakter-levin-12 . It would be interesting in the context of our model to study the response for other classes of Vlasov-stationary reference states that allow a first-order phase transition and differ from the ones studied here. Another immediate follow up of this work would be to investigate the validity, in the context of the studied spin model, of the Kubo fluctuation-dissipation theorem valid for short-range systems prepared in thermal equilibrium and subject to small external fields. Studies in this direction are underway and comparison with the HMF model results yamaguchi-16 will be reported elsewhere. Collective $1/f$ fluctuation due to the Casimir invariants is also an interesting topic to pursue yamaguchi-kaneko-18 .

Acknowledgements.

The work of D.D. is supported by UGC-NET Research Fellowship Sr. No. 2121450744, Ref. No. 21/12/2014(ii) EU-V. Y.Y.Y. acknowledges the supports of JSPS KAKENHI Grant No. 16K05472. The manuscript was finalized while S.G. was visiting the Quantitative Life Sciences section of the International Centre for Theoretical Physics (ICTP), Trieste, and he would like to acknowledge the support and hospitality of the ICTP. S.G. acknowledges support from the Science and Engineering Research Board (SERB), India, Grant No. TAR/2018/000023.

Appendix A Derivation of the critical exponents (34)

Here we discuss how one may obtain the values of the critical exponents given in Eq. (34). The starting point is the equivalent of Eqs. (20), (21) and (23) in presence of a constant field $\boldsymbol{h}=(h,0,0)$:

[TABLE]

with

[TABLE]

and

[TABLE]

The last equation gives

[TABLE]

with $S_{x}=\sqrt{1-p^{2}}\cos\phi$. We assume $G$ to be a smooth function of its argument. Since we take $G$ to be parametrized by the parameter $T$, this means that $G$ is also smooth with respect to $T$. Let us expand $G$ in a Taylor series around $Dp^{2n}$. Substituting the Taylor series and then integrating over $\phi$, we see that in the numerator on the right hand side of Eq. (66), only odd order terms in $(m_{h}+h)$ survive; In the denominator, on the contrary, only even order terms survive. Consequently, the right hand side of Eq. (66) has only odd order terms of $m_{h}$, and the equation gives

[TABLE]

with

[TABLE]

and we have neglected the higher-order terms in $(m_{h}+h)$ in obtaining Eq. (67). Consider now the latter for $h=0$. Assuming $B(T)>0$, one has a non-zero solution for $m_{0}$ 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 solution, giving the spontaneous magnetization for $T<T_{c}$, is $m_{0}=\sqrt{-A/B}$. The smoothness of $G$ with respect to the parameter $T$ allows $A(T)$ to be expanded in a Taylor series around $T_{c}$, giving $A(T)\propto(T-T_{c})$ close to $T_{c}$, and this gives $\beta=1/2$.

The linear response is obtained by deriving Eq. (67) with respect to $h$, and we have

[TABLE]

From the behavior $A(T)\propto(T-T_{c})$ around $T=T_{c}$, we get $\gamma^{\pm}=1$.

At the critical point, the self-consistent equation (67) reduces to

[TABLE]

The response $m_{h}\propto h^{1/\delta}~{}(\delta>1)$ is larger than $h$ for small $h$, so that we have $m_{h}^{3}\propto h$, implying $\delta=3$.

Appendix B Derivation of the response formula,

Eq. (38)

In this appendix, we summarize the derivation of the response formula (38) by following Ref. ogawa-yamaguchi-14 . Noting that the Vlasov equation is governed by the single-spin Hamiltonian $H$, which depends on $f$ through the magnetization $\boldsymbol{m}$, the idea is to expand the Hamiltonian $H$ as

[TABLE]

where $H_{h}$ defined in (28) is the asymptotic part, characterizing the stationary ($t\to\infty$) state $f_{h}$, while $K$ is the transient part. For our spin model, the explicit form of the transient part is

[TABLE]

where the transient magnetization $m_{\rm T}(t)$ is obtained as

[TABLE]

and the transient state $g$ is defined by $g=f-f_{h}$. The transient quantities, $g,K$ and $m_{\rm T}(t)$, are not known a priori, but they do not appear in the final result of the response formula.

Let us write the Vlasov equation (14) as

[TABLE]

where the linear operator $\mathcal{L}_{H}$ is defined as

[TABLE]

Eq. (74) is still exact. Now, we assume that contribution from the transient part, $\mathcal{L}_{K}f$, is negligible, which is justified under some assumptions for Hamiltonian systems and may be related to the phenomenon of Landau damping, see Ref. ogawa-yamaguchi-14 for details.

Under the aforementioned assumption, the formal solution to the Vlasov equation (74) is

[TABLE]

which represents temporal evolution of $f_{0}$ under Hamiltonian flow associated with the asymptotic Hamiltonian $H_{h}$. Assuming ergodicity, a formula that replaces the time average with a partial phase-space average with respect to a iso-$H_{h}$ surface gives

[TABLE]

Noting that the left hand side is nothing but the asymptotic stationary state $f_{h}$, we obtain the response formula (38).

Appendix C Derivation of the estimation (53) for the quantity $L_{2}$

Here, we derive Eq. (53). For simplicity of notation, we use the same symbols $(w,I)$ for angle-action variables associated with the single-spin Hamiltonian $H_{0}$ as the ones used for $H_{h}$ in Sec. IV.2, but the latter do not appear in this section and no confusion should arise. Here we consider positive $n$.

We start with Eq. (52). Noting that $\langle S_{x}\rangle_{0}^{2}$ and $H_{0}$ depend on only the action variable $I$, we rewrite $L_{2}$ as

[TABLE]

where

[TABLE]

and we have used the fact that canonical transformation from $(\phi,p)$ to $(w,I)$ gives ${\rm d}\phi{\rm d}p={\rm d}w{\rm d}I$. The single-spin Hamiltonian $H_{0}$ has a separatrix which consists of the stable and unstable manifolds of the fixed point $(S_{x},S_{y},S_{z})=(-1,0,0)$ and encloses the point $O=(1,0,0)$ on the phase space of the unit sphere. On the basis of this observation, we make an essential assumption for estimating $L_{2}$, namely, that the main contribution to $L_{2}$ comes from the region around the point $O$. We now change twice the variables of integration in $L_{2}$. First, to divide the phase space into the inside and the outside of the separatrix, the integration variable $I$ is changed to energy $E$ as

[TABLE]

where the frequency $\Omega(I)$ is defined by

[TABLE]

Let $E_{\rm min},E_{\rm sep}$ and $E_{\rm max}$ denote respectively the minimum energy, the separatrix energy, and the maximum energy. Following the essential assumption, we omit the second term of $L_{2}$ in (80). Second, to eliminate the dependence on $m_{0}$ of the integration interval, $E$ is changed to a variable $k$ defined as

[TABLE]

Consequently, one has

[TABLE]

where $\psi(I(k))$ is simply denoted as $\psi(k)$ for instance.

To estimate $\Omega$ around the point $O$, which corresponds to $(\phi,p)=(0,0)$, we approximate $H_{0}$ as

[TABLE]

The action variable is the area enclosed by a periodic orbit, and hence, we have

[TABLE]

with $p_{\rm max}=(E/D)^{1/(2n)}$. In terms of

[TABLE]

we have the action variable as

[TABLE]

which gives

[TABLE]

The frequency $\Omega$ is therefore

[TABLE]

Putting all together, we have the estimation of $L_{2}$ as

[TABLE]

We remark that $\langle S_{x}\rangle_{0}$ is zero in the disordered phase, but it does not vanish in the ordered phase even when the limit $m_{0}\to 0$ is taken, because the iso-$H_{0}$ line is confined to the direction $\phi$ around the point $O$ corresponding to $(\phi,p)=(0,0)$. Equation (90) yields the estimation (53) of the main text.

Appendix D A method to compute averages $\langle\cdot\rangle_{0}$ over angles

In this appendix, we discuss a method to compute the angle average $\langle\cdot\rangle_{0}$. The single-spin Hamiltonian $H_{0}$ has the angle-action variables $(w,I)$ whose temporal evolution is

[TABLE]

where $\Omega(I)={\rm d}H_{0}/{\rm d}I$. Changing the variable from $w$ to $t$, we have the average as

[TABLE]

where $T_{p}$ is the period on the considered iso-$I$ line. We may write down the canonical equations of motion for $H_{0}=Dp^{2n}-m_{0}\sqrt{1-p^{2}}\cos\phi$ as

[TABLE]

and the average $\langle\sqrt{1-p^{2}}\cos\phi\rangle_{0}$ is computed as

[TABLE]

Note that the left hand side depends on the action $I$ only, with the initial condition $(\phi(0),p(0))$ determining the value of $I$.