Nondivergent and negative susceptibilities around critical points of a long-range Hamiltonian system with two order parameters

TL;DR

This paper investigates the linear response of a long-range Hamiltonian system with two order parameters near critical points, revealing suppressed responses, non-diverging susceptibilities, and negative off-diagonal elements, confirmed by numerical simulations.

Contribution

It develops a Vlasov-based linear response theory for systems with two order parameters, uncovering non-divergent susceptibilities and negative responses at critical points.

Findings

Some susceptibility matrix elements do not diverge at critical points.

Negative off-diagonal susceptibility elements can occur, indicating inverse responses.

Theoretical predictions are validated by numerical Vlasov simulations.

Abstract

The linear response is investigated in a long-range Hamiltonian system from the view point of dynamics, which is described by the Vlasov equation in the limit of large population. Due to existence of the Casimir invariants of the Vlasov dynamics, an external field does not drive the system to the forced thermal equilibrium in general, and the linear response is suppressed. With the aid of a linear response theory based on the Vlasov dynamics, we compute the suppressed linear response in a system having two order parameters, which introduce the conjugate two external fields and the susceptibility matrix of size two accordingly. Moreover, the two order parameters bring three phases and the three types of second-order phase transitions between two of them. For each type of the phase transitions, all the critical exponents for elements of the susceptibility matrix are computed. The critical exponents reveal that some elements of the matrices do not diverge even at critical points, while the mean-field theory predicts divergences. The linear response theory also suggests appearance of negative off-diagonal elements, in other words, an applied external field decreases the value of an order parameter. These theoretical predictions are confirmed by direct numerical simulations of the Vlasov equation.

Tables

Table 1. Table 1: Critical temperature and critical energy of second-order phase transitions. The critical energy of the Nematic-Ferro phase transition is ϵ NF = ( 2 + Λ 1 ) / 4 + m 20 ​ ( Λ 1 − 2 ​ Λ 2 ​ m 20 ) / 4 subscript italic-ϵ NF 2 subscript Λ 1 4 subscript 𝑚 20 subscript Λ 1 2 subscript Λ 2 subscript 𝑚 20 4 \epsilon_{\rm NF}=(2+\Lambda_{1})/4+m_{20}(\Lambda_{1}-2\Lambda_{2}m_{20})/4 .

	Critical temperature	Critical energy
Para-Ferro	$T={\Lambda }_1/2$	${ϵ}_0=\left(2+{\Lambda }_1\right)/4$
Para-Nematic	$T={\Lambda }_2/2$	${ϵ}_0=\left(2+{\Lambda }_2\right)/4$
Nematic-Ferro	$T={\Lambda }_1\left(1+m_{20}\right)/2$	${ϵ}_0={ϵ}_{\mathrm{NF}}$

Table 2. Table 2: τ 𝜏 \tau dependence of the spontaneous order parameters m 10 , m 20 , m 30 subscript 𝑚 10 subscript 𝑚 20 subscript 𝑚 30 m_{10},m_{20},m_{30} and m 40 subscript 𝑚 40 m_{40} , where τ 𝜏 \tau is the parameter distance from the critical point. PF and FP represent, for instance, the Para side and the Ferro side of the Para-Ferro phase transition, respectively.

	$m_{10}$	$m_{20}$	$m_{30}$	$m_{40}$
PF	0	0	0	0
FP	$O\left({\tau }^{1/2}\right)$	$O\left(\tau \right)$	$O\left({\tau }^{3/2}\right)$	$O\left({\tau }^2\right)$
PN	0	0	0	0
NP	0	$O\left({\tau }^{1/2}\right)$	0	$O\left(\tau \right)$
NF	0	$O\left(1\right)$	0	$O\left(1\right)$
FN	$O\left({\tau }^{1/2}\right)$	$O\left(1\right)$	$O\left({\tau }^{1/2}\right)$	$O\left(1\right)$

Equations

$$\begin{array}{rl}{H}_N& =\sum _{j=1}^N\frac{p_j^2}{2}+\frac{1}{2N}\sum _{j,k=1}^N\Phi (q_j-q_k)-\sum _{a=1}^2{h}_a\sum _{j=1}^N\cos a{q}_j.\end{array}$$

$$\Phi (q)=1-{\Lambda }_1\cos q-{\Lambda }_2\cos 2q$$

$$\frac{1}{N}\sum _{j=1}^N(\cos q_j,\sin q_j),\frac{1}{N}\sum _{j=1}^N(\cos 2q_j,\sin 2q_j),$$

$$M_1=\frac{1}{N}\sum _{j=1}^N\cos q_j,M_2=\frac{1}{N}\sum _{j=1}^N\cos 2q_j,$$

$$\frac{\partial f}{\partial t}+\frac{\partial \mathcal{H}[f]}{\partial p}\frac{\partial f}{\partial q}-\frac{\partial \mathcal{H}[f]}{\partial q}\frac{\partial f}{\partial p}=0,$$

$$\mathcal{H}[f](q,p,t)=\frac{p^2}{2}+\mathcal{V}[f](q,t),$$

$$\mathcal{V}[f](q,t)=\int \Phi (q-q^{\prime })f(q^{\prime },p^{\prime },t)d{q}^{\prime }d{p}^{\prime }.$$

$${\mathcal{M}}_a[f](t)={\iint }_{\mu }\cos aqf(q,p,t)dqdp,(a=1,2)$$

$$\mathcal{V}[f](q,t)=1-\sum _{a=1}^2({\Lambda }_a{\mathcal{M}}_a[f](t)+h_a)\cos aq.$$

$$f_0(q,p)=\frac{G(H_0,\beta )}{{\iint }_{\mu }G(H_0,\beta )dqdp},G(E,\beta )=e^{-\beta E},$$

$$H_0=\mathcal{H}[f_0]=\frac{p^2}{2}+V_0(q),$$

$$V_0(q)=\mathcal{V}[f_0]=1-{\Lambda }_1{m}_{10}\cos q-{\Lambda }_2{m}_{20}\cos 2q.$$

$$m_{10}={\mathcal{M}}_1[f_0],m_{20}={\mathcal{M}}_2[f_0].$$

$$\begin{array}{rl}\text{Para:}& m_{10}=0,m_{20}=0,\\ \text{Nematic:}& m_{10}=0,m_{20}>0,\\ \text{Ferro:}& m_{10}>0,m_{20}>0.\end{array}$$

$$f^{\mathrm{c}}(q,p)=\frac{G(H^{\mathrm{c}},\beta )}{{\iint }_{\mu }G(H^{\mathrm{c}},\beta )dqdp}$$

$$m_a^{\mathrm{c}}={\mathcal{M}}_a[f^{\mathrm{c}}],(a=1,2).$$

$$m_a^{\mathrm{c}}=m_{a0}+\delta m_a^{\mathrm{c}},(a=1,2),$$

$$H^{\mathrm{c}}=H_0+\delta V^{\mathrm{c}}$$

$$\delta V^{\mathrm{c}}=-({\Lambda }_1\delta m_1^{\mathrm{c}}+h_1)\cos q-({\Lambda }_2\delta m_2^{\mathrm{c}}+h_2)\cos 2q.$$

$$f^{\mathrm{c}}\simeq f_0-\beta f_0\left(\delta V^{\mathrm{c}}-{⟨\delta V^{\mathrm{c}}⟩}_0\right)+O({\left|\left|\mathbf{h}\right|\right|}^2).$$

$${⟨X⟩}_0={\iint }_{\mu }X(q,p)f_0(q,p)dqdp.$$

$$\delta {\mathbf{m}}^{\mathrm{c}}=(D^{\mathrm{c}}{)}^{-1}\beta C^{\mathrm{c}}\mathbf{h},$$

$$D^{\mathrm{c}}=1_2-\beta C^{\mathrm{c}}\Lambda ,$$

$$C_{ab}^{\mathrm{c}}={⟨\cos aq\cos bq⟩}_0-{⟨\cos aq⟩}_0{⟨\cos bq⟩}_0,$$

$$\delta {\mathbf{m}}^{\mathrm{c}}=\left(\begin{array}{c}\delta m_1^{\mathrm{c}}\\ \delta m_2^{\mathrm{c}}\end{array}\right),1_2=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\Lambda =\left(\begin{array}{cc}{\Lambda }_1& 0\\ 0& {\Lambda }_2\end{array}\right).$$

$${\chi }_{ab}^{\mathrm{c}}=\underset{||\mathbf{h}||\to 0}{\lim }\frac{\partial \delta m_a^{\mathrm{c}}}{\partial h_b},$$

$${\chi }^{\mathrm{c}}=(D^{\mathrm{c}}{)}^{-1}\beta C^{\mathrm{c}}=(D^{\mathrm{c}}{)}^{-1}(1_2-D^{\mathrm{c}}){\Lambda }^{-1}.$$

$$m_{a0}={\iint }_{\mu }\cos aq{f}_0(q,p)dqdp(a=3,4),$$

$$C^{\mathrm{c}}=\left(\begin{array}{cc}\frac{1+m_{20}}{2}-m_{10}^2& \frac{m_{10}+m_{30}}{2}-m_{10}{m}_{20}\\ \frac{m_{10}+m_{30}}{2}-m_{10}{m}_{20}& \frac{1+m_{40}}{2}-m_{20}^2\end{array}\right).$$

$$D^{\mathrm{c}}=\left(\begin{array}{cc}1-\beta {\Lambda }_1\frac{1+m_{20}}{2}& 0\\ 0& 1-\beta {\Lambda }_2\left(\frac{1+m_{40}}{2}-m_{20}^2\right)\end{array}\right).$$

Full text

I Introduction

The phase transition is one of the central issues in the field of many-body systems. It is classified into some universality classes, and in particular, the mean-field universality class is easily understood by the Landau’s phenomenological theory landau-37 . Nevertheless, a new aspect of the mean-field universality class is recently revealed by considering dynamics.

Dynamics of the mean-field class, including the systems having long-range interaction campa-giansanti-moroni-00 ; mori-10 ; campa-dauxois-ruffo-09 , is described by the VLF’s equation, or the collisionless Boltzmann equation, in the limit of large population braun-hepp-77 ; dobrushin-79 ; neunzert-84 ; spohn-91 . The Vlasov equation has the infinite number of Casimir invariants, and these invariants may prevent the system from relaxing to thermal equilibrium. Indeed, when the initial state has different values of the Casimir invariants from ones in thermal equilibrium, then the relaxation is impossible. We note that, with finite population, the finite-size fluctuation plays the role of collision and drives the system to thermal equilibrium, while the relaxation time gets longer as the population increases zanette-montemurro-03 ; yamaguchi-04 ; barre-06 ; moyano-anteneodo-06 ; jain-bouchet-mukamel-07 ; debuyl-mukamel-ruffo-11 ; chavanis-12 .

The Casimir invariants hold even when an external field is applied, and the invariants suppress the response mazur-69 ; suzuki-71 . This suppression may induce reduction of the critical exponent for the linear response in the Vlasov dynamics. In a ferromagnetic body, the critical exponents $\gamma_{\pm}$ of susceptibility $\chi$ are defined as $\chi\propto\tau^{-\gamma_{\pm}}$ around the second order phase transition. Here $\tau$ is the parameter distance from the critical point like $|T-T_{\rm c}|$ with temperature $T$ and its critical value $T_{\rm c}$, and $\gamma_{+}$ ($\gamma_{-})$ is defined in the paramagnetic (ferromagnetic) phase. The classical values of $\gamma_{\pm}$ in the mean-field universality class are $\gamma_{\pm}=1$. However, in the Vlasov dynamics of the Hamiltonian mean-field (HMF) model inagaki-konishi-93 ; antoni-ruffo-95 , which is a paradigmatic toy model of a ferromagnetic body in the mean-field class, the linear response theory gives $\gamma_{+}=1$ patelli-gupta-nardini-ruffo-12 ; ogawa-yamaguchi-12 but $\gamma_{-}=1/4$ ogawa-patelli-yamaguchi-14 . The nonclassical critical exponent is not restricted in the HMF model, and the universality is discussed for spatially periodic one-dimensional systems ogawa-yamaguchi-15 .

In the HMF model, detection of nonclassical critical exponents is extended to the nonlinear response at the critical point ogawa-yamaguchi-14 and to the correlation length yamaguchi-16 , which is generalized to the infinite-range models by introducing the coherent number of particles botet-jullien-pfeuty-82 ; botet-jullien-83 . Interestingly, the nonclassical critical exponents share some scaling relations with the classical critical exponents.

Another direction of detecting nonclassical critical exponents is to consider the linear response in extended models. In this article, we consider the so-called generalized Hamiltonian mean-field (GHMF) model teles-12 . In the HMF model, particles are confined on the unit circle, and interaction potential consists of the spatial first Fourier mode only. Introducing the second Fourier mode, the GHMF model acquires the Nematic phase in addition to the paramagnetic (Para) and the ferromagnetic (Ferro) phases. As a result, the GHMF model has the new two phase transitions: the Para-Nematic and the Nematic-Ferro phase transitions. As observed in the HMF model, the critical exponents in the Vlasov dynamics may differ between the two sides of a phase transition, and hence we need to consider six sides for the three phase transitions. Moreover, the susceptibility in one side is described by a $2\times 2$ matrix, since the three phases are characterized by the two order parameters corresponding to the two Fourier modes and each order parameter has the conjugate external field. Consequently, we must consider $6$ critical exponent matrices of the size $2\times 2$ and the total number of $\gamma$ is $24$ accordingly.

Appearance of the Nematic phase and the matrix form of the susceptibility give natural questions: Does the appearance of the Nematic phase drastically change the critical exponents from the HMF model? Are there any differences in the off-diagonal elements of the critical exponent matrix between the mean-field theory and the Vlasov dynamics? The purpose of this paper is to answer to these questions. We compute the $24$ critical exponents theoretically by using a response theory based on the Vlasov dynamics ogawa-yamaguchi-15 , and confirm theoretical predictions by performing direct numerical simulations of the Vlasov dynamics. In the HMF model the reduction of the critical exponent is observed, but we show a stronger result in the GHMF model that some elements of the susceptibility matrices do not diverge at the critical point, even they diverge in the mean-field theory. Further, close to the critical point, we demonstrate that the off-diagonal elements of susceptibility matrix become negative in the Ferro side of the Para-Ferro phase transition.

This article is constructed as follows. The GHMF model and the three phases are introduced in Sec.II with the corresponding Vlasov equation. Responses in statistical mechanics and in the Vlasov dynamics are derived in Secs.III and IV, respectively. Theoretical predictions are examined numerically in Sec.V. The last section VI is devoted to a summary and discussions.

II Generalized Hamiltonian mean-field model

II.1 The model

The GHMF model represents particles confined on the unit circle and is described by the Hamiltonian

[TABLE]

The position of $j$-th particle is $q_{j}\in(-\pi,\pi]$, and $p_{j}\in\mathbb{R}$ is the conjugate momentum. $h_{1}$ and $h_{2}$ represent strength of the external fields. The interaction potential $\Phi$ is

[TABLE]

where $\Lambda_{1}$ and $\Lambda_{2}$ are non-negative constants. Setting $\Lambda_{1}=1$ and $\Lambda_{2}=0$, and restricting $h_{2}=0$, the GHMF model results to the HMF model with the external field $h_{1}$. The coefficients are originally defined as $\Lambda_{1}=\Delta$ and $\Lambda_{2}=1-\Delta$ with $\Delta\in[0,1]$ to ensure the attractive interaction, but we slightly restrict the parameter interval as $\Lambda_{1}\in(0,1)$ and $\Lambda_{2}=1-\Lambda_{1}\in(0,1)$ for later convenience. Corresponding to the two Fourier modes in $\Phi(q)$, the two order parameter vectors are defined as

[TABLE]

but we may set the sine parts to be zeros from the rotational symmetry of the system and omit them accordingly. The remaining parts,

[TABLE]

are conjugate to the external field $h_{1}$ and $h_{2}$ respectively.

In the limit $N\to\infty$, dynamics is described by the Vlasov equation

[TABLE]

where the one-particle distribution function $f(q,p,t)$ is defined on the two-dimensional phase space $\mu=(-\pi,\pi]\times\mathbb{R}$. The one-particle Hamiltonian functional $\mathcal{H}[f]$ is defined by

[TABLE]

where the potential functional is

[TABLE]

Omitting the sine part in $\mathcal{V}[f]$ again from the rotational symmetry of the system, and introducing the order parameter functionals defined by

[TABLE]

the potential functional is rewritten as

[TABLE]

II.2 Three phases in unforced equilibrium state

The canonical thermal equilibrium states with zero external fields, $h_{1}=h_{2}=0$, are written as

[TABLE]

where $\beta$ is the inverse temperature,

[TABLE]

and

[TABLE]

The values of $m_{10}$ and $m_{20}$ are determined by solving the simultaneous self-consistent equations

[TABLE]

Note that the right-hand-sides depend on $m_{10}$ and $m_{20}$ through $H_{0}$. We may set $m_{10}$ and $m_{20}$ as non-negative from the rotational symmetry of the system.

The three phases of the GHMF model are characterized as

[TABLE]

On the two-dimensional phase space $\mu$, the separatrix is an iso-$H_{0}$ contour, and forms the skeleton of the phase space. The Para phase has no separatrix, since the iso-$H_{0}$ contours in the Para phase coincide with the iso-$p$ contours. On the other hand, the Ferro and the Nematic phases have separatrices as schematically shown in Fig.1.

III Response in statistical mechanics

Before going to the linear response in the Vlasov dynamics, we revisit the linear response in statistical mechanics for comparison. The Vlasov dynamics corresponds to the microcanonical ensemble, but the microcanonical one gives the equivalent phase diagram with the canonical one except for the parameter region where the first-order phase transition exists pikovsky-14 . We are interested in the susceptibility around the second-order phase transition, and hence we discuss on the response in the canonical ensemble for simplicity.

The susceptibility is defined in Sec. III.1 by applying constant external fields, $h_{1}$ and $h_{2}$. The critical lines are discussed in Sec. III.2 based on the divergence of the susceptibility. The critical exponent matrices are obtained in Sec. III.3.

III.1 Susceptibility

Let the system be in the canonical thermal equilibrium state $f_{0}$ in the time interval $t<0$. We apply constant external fields $h_{1}$ and $h_{2}$ at the time $t=0$ and wait a long time. Then the system is expected to relax to the forced canonical thermal equilibrium state

[TABLE]

with the forced equilibrium Hamiltonian $H^{\rm c}=\mathcal{H}[f^{\rm c}]$. The relaxation is not always true in the Vlasov dynamics, but we consider the state $f^{\rm c}(q,p)$ in this section. The order parameters in $f^{\rm c}$ are denoted by

[TABLE]

Expanding the order parameters as

[TABLE]

the Hamiltonian $H^{\rm c}$ is also expanded as

[TABLE]

with the discrepancy of potential

[TABLE]

When the external field $\boldsymbol{h}=(h_{1},h_{2})^{\rm T}$ is small, where the superscript T represents the transposition, the discrepancy $\delta V^{\rm c}$ is also small and $f^{\rm c}$ is expanded as

[TABLE]

Here the symbol $\left\langle X\right\rangle_{0}$ represents the average of the observable $X(q,p)$ over $f_{0}$ as

[TABLE]

Multiplying (20) by $\cos aq$ and integrating over $\mu$, we have the self-consistent equations and their formal solution

[TABLE]

where the matrix $D^{\rm c}$ is defined by

[TABLE]

the $(a,b)$-element of the matrix $C^{\rm c}$ is defined by

[TABLE]

and

[TABLE]

The susceptibility matrix $\chi^{\rm c}=(\chi^{\rm c}_{ab})$ is defined by

[TABLE]

and the response formula (22) gives

[TABLE]

III.2 Critical lines

Extending the number of order parameters as

[TABLE]

the matrix $C^{\rm c}$ is expressed as

[TABLE]

On the three critical lines, the order parameter $m_{10}$ is always zero, which induces $m_{30}=0$ by the parity of the mode numbers, and the matrix $D^{\rm c}$ can be reduced to

[TABLE]

The critical point has $\det D^{\rm c}=0$, which determines the critical inverse temperature $\beta$ for fixed $\Lambda_{1}$ and $\Lambda_{2}$, or the critical parameter $\Lambda_{1}$ ($\Lambda_{2}$) for fixed $\beta$ and $\Lambda_{2}$ ($\Lambda_{1}$).

The Para-Ferro and the Nematic-Ferro phase transitions are ruled by the order parameter $m_{10}$, and the Para-Nematic phase transition by $m_{20}$. Therefore, the critical lines are obtained as

[TABLE]

where we used the fact that $m_{20}=m_{40}=0$ on the critical lines of the Para-Ferro and Para-Nematic phase transitions. The value of $m_{20}$ in the Nematic-Ferro phase transition are determined for a given set of $\beta,\Lambda_{1}$ and $\Lambda_{2}$ by solving the self-consistent equations (13) with $m_{10}=0$.

Temperature in the canonical ensemble can be transformed to energy in the microcanonical ensemble. The energy functional is defined by

[TABLE]

where the potential is divided by $2$ to avoid the double counting of pair interactions. The value of $\mathcal{E}[f]$ is conserved in the Vlasov dynamics. The unforced equilibrium value of energy $\epsilon_{0}=\mathcal{E}[f_{0}]$ is related to the temperature $T=1/\beta$ as

[TABLE]

The critical temperature and the critical energy for a given set of $\Lambda_{1}$ and $\Lambda_{2}$ are arranged in Table 1.

III.3 Critical exponent matrix $\gamma^{\rm c}$

The critical exponent matrix $\gamma^{\rm c}=(\gamma^{\rm c}_{ab})$ is defined by

[TABLE]

where $\tau$ is the parameter distance from the critical point. Looking back (27), we find that the divergences of the susceptibility comes from the inverse matrix $(D^{\rm c})^{-1}$, and hence we have to compute $\tau$ dependence of the matrix $D^{\rm c}$.

For later convenience, we decompose the matrix $D^{\rm c}$ into the two parts as

[TABLE]

where

[TABLE]

and

[TABLE]

As shown later, the $A$ part is common to the Vlasov dynamics, but the $B$ part is modified. The estimations of $m_{a0}~{}(a=1,2,3,4)$ are obtained from the self-consistent equations for $m_{10}$ and $m_{20}$, (13), and from the definitions of $m_{30}$ and $m_{40}$, (28). The analyses are given in the Appendix A, and the estimated orders are arranged in Table 2.

We may assume, around the critical lines, the left-hand-sides of (31) are of $O(\tau)$ in general. This assumption and Table 2 give estimations of the matrices $D^{\rm c}$’s as

[TABLE]

where NF and FN represent, for instance, the Nematic side and the Ferro side of the Nematic-Ferro phase transition, respectively. We remark that the orders of elements of the matrix $B^{\rm c}$ are equal to or higher than the matrix $A^{\rm c}$, and the matrix $B^{\rm c}$ is negligible for computing the critical exponent matrices in thermal equilibrium.

Coming back to the formula (27), we have the critical exponent matrices as

[TABLE]

Here we assigned the critical exponent [math] if no divergence appears. These critical exponent matrices are reported in Fig.2 with the critical lines obtained in Sec.III.2.

The microcanonical susceptibility does not exceed the canonical susceptibility due to existence of energy conservation mazur-69 ; suzuki-71 , but the two types of susceptibility share the critical exponents as shown in the HMF model ogawa-patelli-yamaguchi-14 . See Appendix B for the critical exponents in the microcanonical ensemble.

IV Response in Vlasov dynamics

A nonlinear response theory is recently proposed for the Vlasov systems ogawa-yamaguchi-14 . Based on it, a simple response formula has been provided ogawa-yamaguchi-15 , which unifies the nonlinear response theory with the linear response theory patelli-gupta-nardini-ruffo-12 ; ogawa-yamaguchi-12 . The formula is valid under some conditions, but they are satisfied for computing the critical exponents yamaguchi-ogawa-15 . We first review the formula in Sec.IV.1, and the critical exponent matrices are obtained in Sec.IV.2. We further discuss on negative elements of a susceptibility matrix in Sec.IV.3.

IV.1 Response formula

The setting is the same with the case of thermal equilibrium discussed in Sec.III. The initial state is the unforced thermal equilibrium state $f_{0}$ (10) and then a small external field $\boldsymbol{h}$ is applied at the time $t=0$. Under the Vlasov dynamics, the state, however, does not go to the forced thermal equilibrium state $f^{\rm c}$: it is trapped at a nonequilibrium state denoted by $f^{\rm V}$ due to the Casimir invariants of the form $\iint_{\mu}s(f)dqdp$, where $s$ is an arbitrary differentiable function. The response formula predicts $f^{\rm V}$ from $f_{0}$.

The associated one-particle Hamiltonian $H^{\rm V}=\mathcal{H}[f^{\rm V}]$ is integrable since $f^{\rm V}$ is stationary and $H^{\rm V}$ has one degree of freedom accordingly. The integrability introduces the angle-action variables $(\theta,J)$ associated with $H^{\rm V}$, and $H^{\rm V}$ can be written as a function of $J$ only.

Roughly speaking, the response formula is expressed as ogawa-yamaguchi-14 ; ogawa-yamaguchi-15

[TABLE]

where the bracket means the average over the angle variable as

[TABLE]

In other words, $f^{\rm V}$ is obtained by taking time average of the initial state $f_{0}$ under the Hamiltonian flow associated with $H^{\rm V}$.

For obtaining the response, as done in Sec.III.1, we expand the right-hand-side of (40) with the expansion

[TABLE]

where

[TABLE]

The key idea for expanding the right-hand-side of the formula (40) is to use the equality

[TABLE]

which holds from the definition of the partial average $\left\langle\cdot\right\rangle_{H^{\rm V}}$. Substituting $H_{0}=H^{\rm V}-\delta V^{\rm V}$ into the explicit expression

[TABLE]

and expanding the right-hand-side with respect to the small $\delta V^{\rm V}$, we have

[TABLE]

In the way we performed the expansion again by using $H^{\rm V}=H_{0}+\delta V^{\rm V}$. The bracket $\left\langle\cdot\right\rangle_{H_{0}}$ is the average over the angle variable associated with $H_{0}$. We omitted a higher order contribution coming from the replacement of $\left\langle\cdot\right\rangle_{H^{\rm V}}$ with $\left\langle\cdot\right\rangle_{H_{0}}$.

Multiplying (46) by $\cos aq$ and integrating over $\mu$, we have a similar formula for susceptibility with thermal equilibrium (27) as

[TABLE]

but with the different matrix $D^{\rm V}$

[TABLE]

where the $(a,b)$-element of the matrix $C^{\rm V}$ is

[TABLE]

Here we used the equality

[TABLE]

The matrix $D^{\rm V}$ is decomposed into the two parts as

[TABLE]

where the $(a,b)$-element of the matrix $B^{\rm V}$ is

[TABLE]

See the Appendix C for a definite integral formula of each element in a reduced case. The matrix $B^{\rm V}$ results to the matrix $B^{\rm c}$ if we replace the partial average over the angle variable, $\left\langle\cos bq\right\rangle_{H_{0}}$, with the average over $f_{0}$, $\left\langle\cos bq\right\rangle_{0}$. However, existence of the partial average modifies the critical exponents.

IV.2 Critical exponent matrix $\gamma^{\rm V}$

According to the Appendix D, the matrices $B^{\rm V}$’s are estimated as

[TABLE]

The constants $c_{1}$ and $c_{2}$ in $B^{\rm V}({\rm FN})$ are positive and we skip to compute their precise values since they do not contribute to the critical exponents as shown later.

As contrasted with thermal equilibrium case, the matrix $B^{\rm V}$ can partially dominate the matrix $D^{\rm V}$. This domination modifies the estimations of $D^{\rm V}$’s from $D^{\rm c}$’s as

[TABLE]

where $\bar{c}_{j}=\min\{1/2,c_{j}\}~{}(j=1,2)$. Recalling the susceptibility formula (47), we have the critical exponent matrices $\gamma^{\rm V}$ as

[TABLE]

where we assigned the critical exponents [math] when no divergences appear even if $\chi_{ab}$’s go to zeros in the limit $\tau\to 0$. The obtained critical exponent matrices are shown in Fig.3 with stressing the different values from the thermal equilibrium case.

We remark that existence of invariants suppress the response mazur-69 ; suzuki-71 , and hence $\chi^{\rm V}\leq\chi^{\rm c}$. This fact implies that no critical lines exist on the parameter plane except for the ones obtained in canonical statistical mechanics. Consequently, there are no shifts of the critical lines and the new zero critical exponents correctly capture the dynamical obstacle to divergences of susceptibility at the critical point.

IV.3 Negative susceptibility in the Para-Ferro phase transition

We note that the susceptibility matrix $\chi$ is proportional to $D^{-1}(\mathbb{1}_{2}-D)=D^{-1}-\mathbb{1}_{2}$, which is written as

[TABLE]

We consider if the signs of the off-diagonal elements can change around the critical point. The sign of $\det D$ must be fixed in one side of a phase transition around the critical line, since $\det D=0$ appears only on the critical lines. Thus, we concentrate on the off-diagonal elements of the matrix $D$ by focusing on the Ferro side of the Para-Ferro phase transition.

The off-diagonal elements of $D^{\rm c}({\rm FP})$ are proportional to

[TABLE]

and are dominated by the negative first term. Thus, no change of sign is possible around the critical point in thermal equilibrium. Indeed, $\chi_{21}^{\rm c}({\rm FP})\to+\infty$ in the limit $\tau\to 0$ reflecting the positive critical exponent $\gamma_{12}^{\rm c}({\rm FP})=1/2$.

However, in the Vlasov dynamics, the off-diagonal elements of $D^{\rm V}({\rm FP})$, which are proportional to

[TABLE]

may change the signs around the critical point. The second term is of $O(\tau^{1/4})$ [see $B^{\rm V}({\rm FP})$ in (53)], and is positive by considering iso-$H_{0}$ contour and $f_{0}$, while the first term is negative and is of $O(\tau^{1/2})$. Thus, the second term can dominate close to the critical point and change the signs of susceptibility. We will numerically demonstrate the negative susceptibility, i.e. $\chi_{21}^{\rm V}({\rm FP})<0$, in Sec.V.3.

V Numerical tests

The critical exponents $1/4$ of $\gamma^{\rm V}_{11}({\rm FP})$ and $\gamma^{\rm V}_{22}({\rm NP})$ are direct extensions of the HMF model, whose Ferro phase also has the same critical exponent ogawa-patelli-yamaguchi-14 . $\gamma^{\rm V}_{12}({\rm FN})=\gamma^{\rm V}_{21}({\rm FN})=0$ may associate with $\gamma^{\rm V}_{11}({\rm FN})=0$. Therefore, interesting exponents are $\gamma^{\rm V}_{11}({\rm NF})=\gamma^{\rm V}_{11}({\rm FN})=0$ and $\gamma^{\rm V}_{12}({\rm FP})=\gamma^{\rm V}_{21}({\rm FP})=0$. We confirm $\gamma^{\rm V}_{11}({\rm NF})=\gamma^{\rm V}_{11}({\rm FN})=0$ in Sec.V.1 and $\gamma^{\rm V}_{21}({\rm FP})=0$ in Sec.V.2 by direct numerical simulations of the Vlasov equation (5). The negative susceptibility discussed in Sec.IV.3 is also examined in Sec.V.3.

We perform the semi-Lagrangian simulations debuyl-10 with the fixed time step $\Delta t=0.05$. The phase space $\mu$ is truncated into $(-\pi,\pi]\times[-4,4]$, and is divided into the grid size $G\times G$. The initial state is the unforced thermal equilibrium state $f_{0}(q,p)$, (10), and a small external field $h_{1}$ is applied with keeping $h_{2}=0$. We compute temporal evolution of the order parameter values $m_{1}=\mathcal{M}_{1}[f]$ and $m_{2}=\mathcal{M}_{2}[f]$ both for $h_{1}=0$ and for $h_{1}>0$, which are denoted by $m_{a}(t,0)$ and $m_{a}(t,h_{1})$ $(a=1,2)$ respectively at the time $t$. Then, we observe the normalized discrepancy between the two to exclude numerical errors for $h_{1}=0$. That is, we observe the quantities

[TABLE]

The family of initial states is characterized by the inverse temperature $\beta=1/T$, but $\beta$ is just a parameter and is interpreted as energy by the relation (33). All the simulations are performed for the deterministic Vlasov equation (5), and no thermal noise is introduced.

V.1 $\gamma^{\rm V}_{11}({\rm NF})=\gamma^{\rm V}_{11}({\rm FN})=0$

Following the previous work teles-12 , we fix energy as $\epsilon=0.55$ to avoid the first-order phase transition region, and vary $\Lambda_{1}$. The parameter $\tau$ is, therefore, $\tau=|\Lambda_{1}-\Lambda_{1{\rm c}}|$, where the critical value $V_{1{\rm c}}$ and the value of $m_{20}$ at the critical point are computed as

[TABLE]

for $h_{1}=h_{2}=0$. We concentrate on the nondivergence of $\chi_{11}$ at the critical point of the Nematic-Ferro phase transition, which implies $\gamma^{\rm V}_{11}({\rm NF})=\gamma^{\rm V}_{11}({\rm FN})=0$.

The $(1,1)$-element of the matrix $B^{\rm V}({\rm NF})$ is expressed in the integral form as

[TABLE]

where $I_{0}$ is the zeroth order modified Bessel function, and $K(k)$ is the complete elliptic integral of the first kind. $k=0$ and $k=1$ correspond to the minimum energy and the separatrix energy respectively. We used the fact that $\left\langle\cos q\right\rangle_{H_{0}}=0$ in the separatrix outside. See the Appendix C.4 for the derivation of (61). Computing the integral numerically, we have the values of $D^{\rm V}_{11}({\rm NF})$ and $\chi_{11}^{\rm V}({\rm NF})$ as

[TABLE]

at the critical point.

The grid size dependence of $\chi_{11}(t)$ is reported in Fig.4, and the numerical results approach to the theoretical level as the grid gets finer. We also computed $h_{1}$ dependence of $\chi_{11}(t)$ with the grid size $G=512$, and the three numerical curves for $h_{1}=10^{-4},10^{-5}$ and $10^{-6}$ almost collapse within the symbol size of Fig.4 (not shown). We, therefore, conclude that the nondivergence of susceptibility and the finite theoretical level (62) are successfully confirmed at the critical point.

V.2 $\gamma^{\rm V}_{21}({\rm FP})=0$

To avoid the first-order phase transition region again, we set $\Lambda_{1}=0.8$ and $\Lambda_{2}=0.2$, which gives the critical energy $\epsilon_{\rm c}=0.7$, and vary $\epsilon$ below the critical point $\epsilon_{\rm c}$. Thus, the parameter $\tau$ is $\tau=\epsilon_{\rm c}-\epsilon$, since we are in the Ferro, low energy phase. We compute the time averages of $\chi_{11}(t)$ and $\chi_{21}(t)$ in the time window $t\in[200,1000]$. The averaged susceptibilities are reported in Fig.5 as functions of $\epsilon_{\rm c}-\epsilon$ for the three Grid sizes, $G=128,~{}256$ and $512$. The numerical results are in good agreements with an approximate theory, in which $m_{a0}~{}(a=2,3,4)$ are neglected (see the Appendix C for the integral form of each element of the matrix $B^{\rm c}({\rm FP})$ in this approximated case). The critical exponent $\gamma^{\rm V}_{11}({\rm FP})=1/4$ is successfully reproduced as in the HMF model ogawa-patelli-yamaguchi-14 , and no divergence of $\chi_{21}$ to $+\infty$ is also confirmed.

V.3 $\chi^{\rm V}_{21}({\rm FP})<0$ close to the critical point

The susceptibility $\chi_{21}$ in Fig.5 is hidden close to the critical point, since $\chi_{21}$ becomes negative. The negative susceptibility is confirmed as shown in Fig.6 by taking the linear scale for the vertical axis. This observation in the Vlasov dynamics gives a sharp contrast with in thermal equilibrium, in which the susceptibility $\chi_{21}^{\rm c}$ is positive and diverges at the critical point.

VI Summary and discussions

We considered responses to the external fields in the GHMF model. This model has the two order parameters, which characterize the Para, the Ferro and the Nematic phases. In each of thermal equilibrium and of the Vlasov dynamics, we derived $6$ critical exponent matrices corresponding to the two sides of the three phase transitions, where each critical exponent matrix is of $2\times 2$ associated with the two order parameters and their conjugate external fields.

As in the HMF model, the Para phase in the Vlasov dynamics has the same critical exponent matrices with thermal equilibrium. This agreement comes from the fact that both $\left\langle\cos aq\right\rangle_{0}$ and $\left\langle\cos aq\right\rangle_{H_{0}}$ vanish in the Para phase and no dynamical effects appear in the matrix $D^{\rm V}$. In the Ferro side of the Para-Ferro phase transition, and the Nematic side of the Para-Nematic phase transition, we obtained the suppressed critical exponent $\gamma=1/4$ as the straightforward extension from the HMF model, where $\gamma=1$ in statistical mechanics. However, in the Ferro and the Nematic phases, all the other exponents are zeros, and no divergences of susceptibility appear at the critical points. The vanishing critical exponents in the Vlasov dynamics are stronger suppression than the reduced value of the previously mentioned $\gamma=1/4$. These theoretical predictions of no divergences are successfully confirmed by direct numerical simulations of the Vlasov dynamics.

We found two types of nondivergences of susceptibilities: one appears in $\chi_{11}^{\rm V}({\rm NF})$ and $\chi_{11}^{\rm V}({\rm FN})$, and the other in the off-diagonal elements of $\chi^{\rm V}({\rm FP})$. The former type might be understood by the potential barrier formed spontaneously by $m_{20}$. Around the critical point, the potential is $V_{0}\simeq-\Lambda_{2}m_{20}\cos 2q$, and has the two wells centered at $q=0$ (well-1) and $q=\pi$ (well-2). Applying the external field to the direction of $q=0$, particles in the well-2 tend to move to the well-1, but the potential barrier may prevent them from moving, since each particle must conserve energy in the Vlasov dynamics. On the other hand, in the latter type, the nondivergences comes from the domination of $O(\sqrt{m_{10}})$ in the off-diagonal elements of the matrix $D^{\rm V}$, but the thermal equilibrium case also have the leading $O(m_{10})$ terms in the off-diagonal elements of the matrix $D^{\rm c}$. Thus, the mechanism might not be straightforward comparing with the former type.

We remark that a non-divergent susceptibility is reported in the HMF model with a family of spatially homogeneous but asymmetric momentum distributions yamaguchi-15 at the point of stability change. The thermal equilibrium states, discussed in the present article, are symmetric and accept non-homogeneous distributions, and hence the reported non-divergences might have a larger impact than the asymmetric case.

We also revealed that negative elements appear in the susceptibility matrix for the Ferro side of the Para-Ferro phase transition. The negative susceptibility has been reported in the HMF model chavanis-11 ; deninno-fanelli-12 ; debuyl-fanelli-ruffo-12 and in the $\phi^{4}$ model campa-ruffo-touchette-07 ; campa-dauxois-ruffo-09 , but they appear under the energy conservation between with and without the external field chavanis-11 (see also Appendix B), the fixed value of magnetization campa-ruffo-touchette-07 ; campa-dauxois-ruffo-09 , or the nonstationary initial states deninno-fanelli-12 ; debuyl-fanelli-ruffo-12 . The negative susceptibility reported in this article is observed for the initial thermal equilibrium states by applying an external field without any additional constraints, and therefore, it might be rather easy to compare with experiments.

Finally, it might be worth noting that the nonclassical critical exponents of the HMF model are also observed in a model of coupled oscillators by setting the so-called natural frequencies deterministically hong-chate-tang-park-15 . In the model, the oscillators are confined on the unit circle and the interaction is realized only through the first Fourier mode as the HMF model. Thus, it might be interesting to consider a similar extension in the coupled oscillator system.

Acknowledgements.

Y.Y.Y. thanks S.Ogawa for valuable discussions. He acknowledges the supports of JSPS KAKENHI Grant Numbers 23560069 and 16K05472.

Appendix A Estimations of $m_{a0}$

Around the critical point, we estimate the spontaneous order parameters $m_{a0}$’s, which are written as

[TABLE]

where the denominator is of $O(1)$. The key idea of this section is to use the orthogonality of $\{\cos aq\}_{a\in\mathbb{Z}}$, which gives

[TABLE]

It is, therefore, enough to estimate $m_{10}$ and $m_{20}$. We first consider the Para-Ferro and the Para-Nematic phase transitions, around which $m_{10}$ and $m_{20}$ are small enough, and then go to the Nematic-Ferro transition.

A.1 Para-Ferro and Para-Nematic transitions

All the order parameters are zeros in the Para phase, and we focus on estimating the order parameters in the Ferro and the Nematic sides.

The normalization factor, the numerator of the right-hand-side of (63), is expanded as

[TABLE]

Thus, the self-consistent equations, which are $a=1$ and $2$ in (63), are reduced to

[TABLE]

In the Ferro side of the Para-Ferro phase transition, the ordering is $O(m_{20})\leq O(m_{10}^{2})$ ogawa-yamaguchi-15 , and $m_{10}$ is determined by the leading two terms as

[TABLE]

As assumed at the head of Sec.III.3, the coefficient of the first term is of $O(\tau)$. Thus, we obtain $m_{10}=O(\tau^{1/2})$ and $m_{20}=O(m_{10}^{2})=O(\tau)$ from the second equation of (66). Further, the estimations (64) give $m_{30}=O(\tau^{3/2})$ and $m_{40}=O(\tau^{2})$.

In the Nematic side of the Para-Nematic transition $m_{10}$ is always zero, and $m_{20}$ is determined by the equation

[TABLE]

As discussed above, we have $m_{20}=O(\tau^{1/2})$. Further, $m_{30}=0$ and $m_{40}=O(\tau)$ from (64).

A.2 Nematic-Ferro transition

In the Nematic-Ferro transition $m_{20}$ is of $O(1)$, and we need to estimate $m_{10}$. Smallness of $m_{10}$ reduces the normalization factor as

[TABLE]

The self-consistent equation for $m_{10}$ is expanded as

[TABLE]

and the definition of $m_{20}$ gives

[TABLE]

where we used the fact

[TABLE]

Recalling the critical line (31), from (71), we conclude that $m_{10}=O(\tau^{1/2})$. From (64) we also estimate $m_{30}=O(\tau^{1/2})$ and $m_{40}=O(1)$.

Appendix B Critical exponents in microcanonical ensemble

In the microcanonical ensemble, temperature may be modified by applying an external field $\boldsymbol{h}$ at the time $t=0$ due to the energy conservation. Denoting the modified temperature by $T^{\rm mc}$, we consider the energy conservation relation between $t=0^{+}$ and $t\to\infty$ as

[TABLE]

where $m_{1}^{\rm mc}$ and $m_{2}^{\rm mc}$ are the values of order parameters in the microcanonical ensemble. Introducing the response

[TABLE]

which will be determined later, the above relation gives the shift of inverse temperature from $\beta$ to $\beta+\delta\beta$, where

[TABLE]

up to the leading order.

Let us introduce the vectors

[TABLE]

and the discrepancy of potential

[TABLE]

The self-consistent equation in the microcanonical ensemble is obtained by replacing $\beta\delta V^{\rm c}$ with $\beta\delta V^{\rm mc}+H_{0}\delta\beta$ in (20),

[TABLE]

Noting that $p^{2}/2$ term of $H_{0}$ is canceled between the two terms, and substituting (75), (77) and

[TABLE]

into the self-consistent equation (78), we have

[TABLE]

where $C^{\rm c}=\left\langle\boldsymbol{c}\boldsymbol{c}^{\rm T}\right\rangle_{0}$ by the definition. The response $\delta\boldsymbol{m}^{\rm mc}$ in the microcanonical ensemble is, therefore,

[TABLE]

where the matrix $D^{\rm mc}$ is defined by

[TABLE]

The matrix $D^{\rm mc}$ is expressed as

[TABLE]

and the second term of the right-hand-side does not change the dominating $\tau$ dependence of $D^{\rm mc}$ from $D^{\rm c}$. This concludes that the critical exponents are shared between the canonical and the microcanonical ensembles.

We give a remark on usage of the energy conservation. If we require the energy conservation between $t=0^{-}$ and $t\to\infty$, the energy conservation relation is modified from (73) to

[TABLE]

where all the superscripts are replaced to represent the considering situation. Then, $\delta\beta$ is modified to

[TABLE]

In the previous setting, the last term was not proportional to $\boldsymbol{m}_{0}^{\rm T}\boldsymbol{h}$ but to $\delta\boldsymbol{m}^{\rm T}\boldsymbol{h}$ and was omitted since it was of higher order. With the modified $\delta\beta^{\rm ene}$, we have the linear response $\delta\boldsymbol{m}^{\rm ene}$ as

[TABLE]

Divergences of the linear response come from $(D^{\rm mc})^{-1}$ again, and hence the critical exponents are not modified. On the other hand, the response $\delta\boldsymbol{m}^{\rm ene}$ may be negative in the off-diagonal elements due to the factor $\mathbb{1}_{2}-2\beta\Lambda\boldsymbol{m}_{0}\boldsymbol{m}_{0}^{\rm T}$, and even in the diagonal elements for large $\beta$, which implies large $\boldsymbol{m}_{0}$.

Appendix C Integral formula for elements of the matrix $B^{\rm V}$ in a reduced case

We give a useful formula of the matrix $B^{\rm V}$ in the case that the one-particle Hamiltonian is written in the form

[TABLE]

This form includes the Nematic phase by replacing $q$ with $2q$ and setting $M=\Lambda_{2}m_{20}$, and the approximate theory used in Sec.V.2, which is for the Ferro side of the Para-Ferro phase transition, by setting $m_{20}=0$ and $M=\Lambda_{1}m_{10}$. We note that this Hamiltonian is symmetric with respect to both $q$ and $p$, namely $H_{0}(-q,p)=H_{0}(q,-p)=H_{0}(q,p)$. We will use this symmetry for reducing computations.

C.1 Angle-action variables and elliptic integrals

The Hamiltonian system $H_{0}$ (87) is integrable, and we can introduce the associated angle-action variables $(\theta,J)$. They are written in the use of the Legendre elliptic integrals, defined by

[TABLE]

and

[TABLE]

The complete elliptic integrals of the first and the second kinds are defined respectively as

[TABLE]

The Hamiltonian system $H_{0}$ (87) has two hyperbolic fixed points, $(q,p)=(-\pi,0)$ and $(\pi,0)$, and they are connected by the separatrix. The phase space $\mu$ is divided into outside and inside of the separatrix. See Fig.1(a) for a schematic picture of the phase space. Based on this knowledge, the angle-action variables $(\theta,J)$ are introduced as barre-olivetti-yamaguchi-10

[TABLE]

and

[TABLE]

where

[TABLE]

and $Q$ is defined by

[TABLE]

The energy minimum corresponds to $k=0$, and the separatrix energy to $k=1$.

For later convenience, we introduce the integrals

[TABLE]

This integrals have the recursive formula

[TABLE]

and hence

[TABLE]

C.2 Computations of $\left\langle\cos bq\right\rangle_{H_{0}}$

Let us compute the averages

[TABLE]

Using the elliptic function ${\rm sn}$, which is the inverse function of $F(\phi,k)$ and is defined by

[TABLE]

we can write $\cos q$ as ogawa-yamaguchi-14

[TABLE]

Changing the variable as

[TABLE]

and using the symmetry of $H_{0}(q,p)$, we have the expressions of $\left\langle\cos^{n}q\right\rangle_{H_{0}}$ as

[TABLE]

Thus, the integrals $N_{n}(k)$, (97), derive $\left\langle\cos q\right\rangle_{H_{0}}$ and $\left\langle\cos^{2}q\right\rangle_{H_{0}}$. The concrete expressions of $\left\langle\cos q\right\rangle_{H_{0}}$ and $\left\langle\cos 2q\right\rangle_{H_{0}}$, which can be computed from $\left\langle\cos^{2}q\right\rangle_{H_{0}}$, are

[TABLE]

and

[TABLE]

C.3 Computations of $\left\langle\left\langle\cos aq\right\rangle_{H_{0}}\left\langle\cos bq\right\rangle_{H_{0}}\right\rangle_{0}$

We first show the equality (50). Noting that $f_{0}$ and $\left\langle\psi_{2}\right\rangle_{H_{0}}$ depends on $J$ only, and using $dqdp=d\theta dJ$, we can show

[TABLE]

We then consider the average

[TABLE]

As shown previously, the average $\left\langle\cos aq\right\rangle_{H_{0}}$ is obtained as a function of $k$, and accordingly, we change the integral variable from $J$ to $k$ by using the Jacobian

[TABLE]

where we used the derivatives of $K(k)$ and $E(k)$

[TABLE]

Denoting the initial stationary state as $f_{0}(q,p)=G_{0}(k)$, and recalling that the separatrix outside has two contributions from the upper and the lower of the separatrix, we have

[TABLE]

C.4 The $(1,1)$-element in the Nematic phase

We give the $(1,1)$-element in the Nematic phase. We derive it via replacing $q$ with $2q$ in the obtained results. We note that the same formula is also derived by starting from the Hamiltonian

[TABLE]

The Nematic phase has two hyperbolic fixed points of $(q,p)=(-\pi/2,0)$ and $(\pi/2,0)$ and the separatrix connects them by forming two “eyes” centered at $(q,p)=(0,0)$ (eye-1) and $(\pi,0)$ (eye-2). See Fig.1(b).

From symmetry, the average $\left\langle\cos q\right\rangle_{H_{0}}$ is canceled in separatrix outside. Indeed, the average is modified as

[TABLE]

In the eye inside, we have the transform

[TABLE]

and, referring to (100), $\cos 2q$ is expressed as

[TABLE]

Therefore, we totally have

[TABLE]

We have two contributions from the eye-1 and the eye-2, but the factor $2$ is canceled with the factor $1/2$ from the Jacobian $dJ/dk$. Indeed, the action variable defined as

[TABLE]

becomes half since the traveling distance of a periodic orbit becomes half in the separatrix inside of Fig.1(b) comparing with Fig.1(a). We remark that the action in the separatrix outside does not change since the traveling distance does not change.

Putting all together with the thermal equilibrium state

[TABLE]

we have

[TABLE]

This expression results to $B^{\rm V}_{11}({\rm NF})$ (61) by setting $M=\Lambda_{2}m_{20}$.

Appendix D Estimations of $B^{\rm V}$ matrix

We give estimations of the matrix $B^{\rm V}$ (52) by using the formula (109), which implies $B^{\rm V}=O(\sqrt{M})$ by appropriately setting $M$, since the integral part does not vanish even in the limit $M\to 0$. Keeping this ordering in mind, we separately discuss on the Para, the Nematic and the Ferro phases.

D.1 Para phase

All the order parameters are zeros in the Para phase, and the angle variable is nothing but $q$. Thus, we have $\left\langle\cos aq\right\rangle_{H_{0}}=0$ and hence

[TABLE]

This is consistent with setting $M=0$ in the formula (109).

D.2 Nematic phase

The parameter $M$ is regarded as $\Lambda_{2}m_{20}$, and the matrices $B^{\rm V}$ is of $O(m_{20})$. However, the off-diagonal elements vanish due to cancellation. The cancellation can be found as follows. In separatrix outside, we recall $\left\langle\cos q\right\rangle_{H_{0}}=0$ and there is no contribution from the separatrix outside to the off-diagonal elements. In separatrix inside, we have contributions from two eyes (see the Appendix C.4). The contribution from the eye-2 is obtained by shifting the variable $q$ with $\pi$ in the contribution from the eye-1, and is multiplied by $\cos a\pi$. Thus, the total contribution from the two eyes has the prefactor $1+\cos a\pi\cos b\pi$, and vanishes for $(a,b)=(1,2)$ and $(2,1)$. The ordering of $m_{20}$ is $m_{20}=O(\tau^{1/2})$ for the Para-Nematic phase transition, and $m_{20}=O(1)$ for the Nematic-Ferro phase transition. These estimations give

[TABLE]

D.3 Ferro phase

For the Para-Ferro phase transition, we may approximate the potential as

[TABLE]

and hence the parameter $M$ is regarded as $\Lambda_{1}m_{10}$ and is of $O(\tau^{1/2})$. There is no reason of cancellation which occurs in the Nematic phase, and hence we have

[TABLE]

For the Nematic-Ferro phase transition, we may approximate the potential as

[TABLE]

and hence the parameter $M$ is regarded as $\Lambda_{2}m_{20}$ and is of $O(1)$. The approximated potential is the same with one in the Nematic phase, but the cancellation does not exactly occur by symmetry breaking due to non-zero $m_{10}$. The off-diagonal elements may be non-zeros and tend to vanish as approaching to the critical line. However, the off-diagonal elements are not important to observe nondivergence of susceptibility at the critical point of the Nematic-Ferro phase transition, and we skip the precise computations. Consequently, we have

[TABLE]

with $c_{1},c_{2}>0$.