text:
# Thermodynamic limit of the canonical partition function with respect to the quark number in QCD

## I INTRODUCTION

Extensive efforts have been shedding light upon the phase structure of Quantum Chromodynamics (QCD) at finite temperature, though there remain many subtleties yet. The center symmetry is a prosperous implement to characterize the deconfinement transition in pure non-abelian gauge theories (gauge theories without dynamical quarks) sve82 . The order parameter to examine whether the symmetry is broken or not is the expectation value of the Polyakov loop pol78 , that is, the Wilson line winding around the Euclidean thermal torus. The Polyakov loop vanishes in the confined phase, while it takes a finite value in the deconfined phase. In pure gauge theories the critical properties associated with the deconfinement transition have been ascertained by the lattice Monte-Carlo simulations. The results are in accord with those anticipated from the center symmetry and the universality argument mcl81 .

In contrast to pure gauge theories any definite indicator in order to distinguish the deconfined phase from confined one is not established so far yet for systems including dynamical quarks in the fundamental representation. Once the thermal excitation of light quarks is allowed, the center symmetry is broken explicitly. If we construct a 3-d effective model in terms of the order parameter, namely the Polyakov loop, in the presence of dynamical quarks, those quark contributions bring about an external magnetic-like field acting onto the Polyakov loop, which breaks the center symmetry ban83 . Then it is obvious that the Polyakov loop is no longer a proper indicator for confinement because it remains finite even in the confined phase as well as in the deconfined phase due to the absence of the center symmetry.

We note, however, that the notion of confinement should be still articulate even when dynamical quarks are present simply because it is an experimental fact. To find out an appropriate indicator to identify the confined phase in the presence of dynamical quarks it is essential to clarify where the dynamical screening against the Polyakov loop should originate from. In fact, DeTar and McLerran det82  argued that some fractional excitations of dynamical quarks are responsible for the explicit breaking of the center symmetry and proposed a new order parameter for the deconfinement transition based on the canonical description with respect to the quark triality or quark number. This idea had been deserted since Meyer-Ortmanns mey84  showed that DeTar–McLerran’s order parameter behaves unexpectedly and always indicates the confining character for the Z(2) Higgs model on the lattice. Although Meyer-Ortmanns has explicitly demonstrated a failure in the proposed order parameter in the model study, it is not apparent in principle what would cause DeTar–McLerran’s formulation not to fulfill the naive expectation as a proper order parameter. Similar ideas for the deconfinement order parameter based on the canonical ensemble have been proposed and discussed by several authors also wei87 ; ole92 ; fab95  but those arguments are essentially integrated into DeTar–McLerran’s first insight.

Interestingly enough, we can find a clear-sighted comment in Meyer-Ortmanns’ paper as follows: “The states at zero temperature have a finite particle number hence zero particle density, while the relevant states at finite temperature have finite particle density. Therefore one cannot check whether there exist states at zero temperature with fractional baryon number by looking at what states contribute to the partition function at finite temperature.” In her analysis on the Z(2) Higgs model, we can discern no clear embodiment for the statement quoted above. Nevertheless, we would emphasize that this point is actually the most critical in constructing the order parameter for the deconfinement transition. The purpose of the present paper is to give a simple and transparent demonstration to disclose what becomes of the partition function with the particle number kept fixed under the thermodynamic limit.

The problem here can be also stated in more general grounds. The question is whether the grand canonical description with zero quark chemical potential is equivalent to the canonical description with zero quark number or not. We will find that the canonical description with zero quark number would amount to the canonical description with zero quark density whenever the thermodynamic limit is taken. As a result, DeTar–McLerran’s idea to project out the fractional excitations of dynamical quarks does not work because the states with zero particle density may have any fractional excitation of particles which becomes irrelevant eventually in the thermodynamic limit. Then by all means can we reach the desired canonical ensemble in which the particle number is fixed at zero? We will propose one possibility in the present work.

This paper is organized as follows: In Sec. II we review DeTar–McLerran’s idea from the point of view of the canonical ensemble with respect to the quark number. Sec. III is devoted to investigating the thermodynamic limit of the canonical partition function. Using the simplest example we elucidate analytically and numerically that the canonical partition function would be reduced to the conventional grand canonical one under the thermodynamic limit. We reveal the cause of the failure in defining an order parameter for the deconfinement transition. In Sec. IV, we propose an idea to overcome the problem in connection with taking the thermodynamic limit. Performing the Monte-Carlo simulation, we put the idea to the test to confirm that our order parameter would be prosperous. We also find that the problem of constructing an order parameter for the deconfinement transition is deeply related to the sign problem of the Dirac determinant. The results of the numerical simulations suffer from the severe statistical errors due to the sign problem. The concluding remarks are in Sec. V.

## II IDEA OF THE CANONICAL ENSEMBLE

### II.1 Canonical ensemble with respect to the quark number

We will briefly review the idea of the canonical ensemble with emphasis on the possibility to define an order parameter for the deconfinement transition.

First of all, we should notice that the conventional QCD partition function at finite temperature and zero density is implicitly formulated in the grand canonical ensemble with zero quark chemical potential, that is equivalent to the canonical ensemble with zero quark density under the thermodynamic limit. This fact is easily understood as follows. In the imaginary-time formalism of the finite-temperature field theories, the temporal (thermal) extent is compactified to the inverse temperature. Then, in addition to the pair excitations of a quark and an antiquark as depicted in Fig. 2 (a), isolated quark excitations as shown in Fig. 2 (a) also satisfy the quark current conservation and thus the quark number can fluctuate thermally. In the thermodynamic limit, however, the quark number conservation is recovered because only the states with zero quark number dominate the partition function at vanishing chemical potential. This point of view has already been discussed in Ref. azc98  and it seems reasonable from the principle of the statistical mechanics.

Usually it does not matter in calculating physical observables whether such fractional excitations of dynamical quarks are present. The situation may be totally different for the manifestation of the center symmetry. In the language of the lattice gauge theory, the center symmetry is defined as the invariance under the center transformation, that is, the transformation from the temporal link variables $`U_4(x_4=N_\tau )`$ on a time slice altogether into $`zU_4(x_4=N_\tau )`$, where $`z`$ is an element of the center of the gauge group (Z(3) for the SU(3) gauge theory). For the quark thermal contributions the Z(3) factor $`z`$ is multiplied every time dynamical quarks wind around the thermal torus. Thus thermal excitations gathered typically in Fig. 2 are center symmetric because $`zz^{}=z^3=1`$ for $`zZ(3)`$. Thermal excitations shown in Fig. 2 are, on the other hand, typical examples which break the center symmetry. Then, what happens if we exclude any thermal excitation like those in Fig. 2 by using the canonical description where the quark number is fixed? Once we employ the canonical ensemble with zero quark number, for example, only such excitations as shown in Fig. 2 (a) (b) are allowed so that the canonical partition function would regain the center symmetry. Certainly the symmetry itself is explicitly broken in the Lagrangian and nevertheless the physical states remain center symmetric. Owing to the recovered symmetry, we can presume that the Polyakov loop will provide a criterion for the deconfinement transition in the same way as in the case of pure gauge theories. Essentially, this is the idea proposed originally by DeTar and McLerran in Ref. det82 .

The canonical description of QCD at finite density is available by the Legendre transformation from the quark chemical potential $`\mu `$ to the quark number $`n_q`$ eng99 . (In this paper we simply write the quark number to mean the quark number minus the antiquark number.) For later convenience, we explain the transformation procedures in some mathematical details. To make our calculation tangible, we adopt the Wilson fermion on the lattice from now on. The final results, of course, can be understood not only in a specific formalism but in a general way also.

The quark number density is the fourth component of the conserved current given by

$$j_4(x)=\kappa \{\overline{\psi }(x)(1\gamma _4)U_4(x)\psi (x+\widehat{4})\overline{\psi }(x+\widehat{4})(1+\gamma _4)U_4^{}(x)\psi (x)\},$$
(1)
where $`\kappa `$ is the hopping parameter and the Wilson parameter is chosen as $`r=1`$ for the sake of simplicity. By imposing the constraint onto the configurations we can write the canonical partition function with respect to the quark number $`n_q`$ as

$`Z_{\text{CE}}(n_q)`$ $`={\displaystyle 𝒟U𝒟\overline{\psi }𝒟\psi \mathrm{e}^{S_\text{G}[U]S_\text{F}[U,\overline{\psi },\psi ]}\delta \left(\frac{1}{N_\tau }\underset{x}{}j_4(x)n_q\right)}`$
$`={\displaystyle _0^{2\pi }}{\displaystyle \frac{\mathrm{d}\varphi }{2\pi }}\mathrm{e}^{\mathrm{i}n_q\varphi }{\displaystyle 𝒟U𝒟\overline{\psi }𝒟\psi \mathrm{e}^{S_\text{G}[U]S_\text{F}[U,\overline{\psi },\psi ]+\mathrm{i}\frac{\varphi a}{\beta }{\scriptscriptstyle j_4(x)}}},`$ (2)
where $`S_\text{G}[U]`$ is the gluon action and $`S_\text{F}[U,\overline{\psi },\psi ]`$ is the action of the Wilson fermion with $`r=1`$. The lattice spacing $`a`$ and the number of lattices along the temporal direction $`N_\tau `$ give the inverse temperature, i.e., $`\beta =N_\tau a`$. The quark parts of the action are put together as

$`S_\text{F}[U,\overline{\psi },\psi ]+\mathrm{i}{\displaystyle \frac{\varphi a}{\beta }}{\displaystyle \underset{x}{}}j_4(x)`$
$`=`$ $`{\displaystyle \underset{x}{}}\{\overline{\psi }(x)\psi (x)\kappa {\displaystyle \underset{j}{}}[\overline{\psi }(x)(1\gamma _j)U_j(x)\psi (x+\widehat{j})+\overline{\psi }(x+\widehat{j})(1+\gamma _j)U_j^{}(x)\psi (x)]`$
$`\kappa [(1+\mathrm{i}{\displaystyle \frac{\varphi a}{\beta }})\overline{\psi }(x)(1\gamma _4)U_4(x)\psi (x+\widehat{4})+(1\mathrm{i}{\displaystyle \frac{\varphi a}{\beta }})\overline{\psi }(x+\widehat{4})(1+\gamma _4)U_4^{}(x)\psi (x)]\}.`$ (3)
Here $`\mathrm{i}\varphi `$ can be regarded as an imaginary chemical potential and the integration over $`\varphi `$ corresponds to the Legendre transformation from the chemical potential to the quark number eng99 . $`\widehat{j}`$ runs from $`\widehat{1}`$ to $`\widehat{3}`$ in the spatial directions. In the same way as the usual prescription to treat the chemical potential in the lattice gauge theories has83  we exponentiate the imaginary chemical potential as follows;

$$1+\mathrm{i}\frac{\varphi a}{\beta }\mathrm{exp}\left(\mathrm{i}\frac{\varphi a}{\beta }\right),1\mathrm{i}\frac{\varphi a}{\beta }\mathrm{exp}\left(\mathrm{i}\frac{\varphi a}{\beta }\right),$$
(4)
which makes only higher order corrections of $`\mathrm{O}(a^2)`$. The fundamental reason to demand the above prescription is because of the gauge invariance, or the center symmetry in the present case. With these alterations the partition function given by Eq. (2) becomes exactly center symmetric for $`n_q=0,\pm 3,\pm 6,\mathrm{}(0mod3)`$ as it should be. In order to see it apparently, transforming the quark fields by

$$\psi ^{}(x)=\mathrm{e}^{\mathrm{i}\frac{\varphi a}{\beta }x_4}\psi (x),\overline{\psi }^{}(x)=\mathrm{e}^{\mathrm{i}\frac{\varphi a}{\beta }x_4}\overline{\psi }(x),$$
(5)
we rewrite the quark action in terms of the transformed fields as eng99

$`\stackrel{~}{S}_\text{F}[U,\overline{\psi }^{},\psi ^{}]`$
$`+\kappa {\displaystyle \underset{\stackrel{}{x}}{}}\left[\mathrm{e}^{\mathrm{i}\varphi }\overline{\psi }^{}(\stackrel{}{x},N_\tau )(1\gamma _4)U_4(\stackrel{}{x},N_\tau )\psi ^{}(\stackrel{}{x},1)+\mathrm{e}^{\mathrm{i}\varphi }\overline{\psi }^{}(\stackrel{}{x},1)(1+\gamma _4)U_4^{}(\stackrel{}{x},N_\tau )\psi ^{}(\stackrel{}{x},N_\tau )\right].`$ (6)
$`\stackrel{~}{S}_\text{F}[U,\overline{\psi },\psi ]`$ is the rest of the action from which the terms involving $`U_4(\stackrel{}{x},N_\tau )`$ are subtracted. It is obvious that the center transformation $`U_4(\stackrel{}{x},N_\tau )\mathrm{e}^{\mathrm{i2}\pi k/3}U_4(\stackrel{}{x},N_\tau )`$ is equivalent to the shift of $`\varphi `$ by $`2\pi k/3`$. We can immediately recognize from Eq. (2) that the canonical partition function is certainly invariant under $`\varphi \varphi +2\pi k/3`$ as long as the quark number satisfies $`n_q=0mod3`$.

The integration of the partition function in terms of $`\varphi `$ exactly corresponds to the decomposition of the Dirac determinant into each part with fixed quark winding number. This is seen transparently if the Dirac determinant obtained from the quark action Eq. (6), that is denoted by $`detD[U,\varphi ]`$ from now on, is calculated by means of the hopping parameter expansion. The expanded terms are represented by closed quark paths similarly as depicted in Figs. 2 and 2. Every time the quark path winds around the thermal torus, it picks up the factor $`\mathrm{e}^{\mathrm{i}\varphi }`$ so that the Dirac determinant can be decomposed according to the quark number $`n_q`$;

$$detD[U,\varphi ]=\underset{n_q=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{e}^{\mathrm{i}n_q\varphi }detD^{(n_q)}[U],$$
(7)
where $`detD^{(n_q)}[U]`$ is the part of the Dirac determinant with the quark path winding $`n_q`$ times around the thermal torus. Thus the integration over $`\varphi `$ in Eq. (2) singles out some specific combinations of the operators coming from the Dirac determinant under the condition that the quark winding number is fixed. This result is not particular to the case of the Wilson fermion formalism on the lattice. In general, the canonical partition function is given by

$$Z_{\text{CE}}(n_q)=_0^{2\pi }\frac{\mathrm{d}\varphi }{2\pi }\mathrm{e}^{\mathrm{i}n_q\varphi }𝒟UdetD[U,\varphi ]\mathrm{e}^{S_\text{G}[U]}$$
(8)
with the Dirac determinant $`detD[U,\varphi ]`$ under the twisted boundary conditions,

$$\psi (\stackrel{}{x},x_4=\beta )=\mathrm{e}^{\mathrm{i}\varphi }\psi (\stackrel{}{x},x_4=0),\overline{\psi }(\stackrel{}{x},x_4=\beta )=\mathrm{e}^{\mathrm{i}\varphi }\overline{\psi }(\stackrel{}{x},x_4=0).$$
(9)
It leads to the canonical partition function in the same sense as in Ref. ole92 , that is, the partition function with the Dirac determinant replaced by its specific part with some fixed quark excitations. It is interesting that the constraint onto the configurations, imposed by Dirac’s delta function, results in the restriction onto the combination of the operators after all. This is because of the non-local character of the quark (fermion) fields.

### II.2 DeTar–McLerran’s order parameter

Here we will explicitly construct DeTar–McLerran’s order parameter for the deconfinement transition. They introduced the quark triality $`t_q`$ defined as $`n_qmod3`$ and considered the canonical partition function with respect to $`t_q`$. It is available from the canonical partition function with respect to $`n_q`$ as

$`\stackrel{~}{Z}_{\text{CE}}(t_q)`$ $`={\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}Z_{\text{CE}}(n_q=3m+t_q)`$
$`={\displaystyle \frac{1}{3}}{\displaystyle \underset{\varphi =0,\pm 2\pi /3}{}}\mathrm{e}^{\mathrm{i}t_q\varphi }{\displaystyle 𝒟UdetD[U,\varphi ]\mathrm{e}^{S_\text{G}[U]}}.`$ (10)
Then DeTar–McLerran’s order parameter is given by

$$\stackrel{~}{𝒫}_𝒩=\frac{\stackrel{~}{Z}_{\text{CE}}(t_q=𝒩)}{_{t_q}\stackrel{~}{Z}_{\text{CE}}(t_q)}$$
(11)
for arbitrary $`𝒩0`$. When the quark number is fixed instead of the quark triality, we can define also in the same way as

$$𝒫_𝒩=\frac{Z_{\text{CE}}(n_q=𝒩)}{_{n_q}Z_{\text{CE}}(n_q)},$$
(12)
which is expected to behave similarly to $`\stackrel{~}{𝒫}_𝒩`$. In the confined phase where all the physical excitations are color-neutral, any thermal excitation of isolated quarks would cost infinitely large energy. Therefore it follows that $`\stackrel{~}{Z}_{\text{CE}}(t_q=𝒩)=Z_{\text{CE}}(n_q=𝒩)=0`$ for $`𝒩0mod3`$. After the dynamical breaking of the center symmetry occurs, fractional quark excitations are allowed at finite energies so that the canonical partition function can take a non-vanishing value. Hence, $`\stackrel{~}{𝒫}_𝒩`$ and $`𝒫_𝒩`$ are anticipated to serve as order parameters for the deconfinement transition.

We can manipulate the order parameter in other forms comparable with the conventional definition, namely the Polyakov loop. As discussed in Refs. ole92 ; fab95 , the expectation value of the Polyakov loop itself can be expected to be a proper order parameter in the canonical ensemble. Because the center symmetry becomes manifest in the canonical partition function with the quark number $`n_q=0mod3`$ (or the quark triality $`t_q=0`$), the deconfinement transition could be characterized by the dynamical breaking of the center symmetry. Thus another candidate for an order parameter for the deconfinement transition is given by $`\stackrel{~}{\mathrm{tr}L(\stackrel{}{x})}_{\text{CE}}`$ with $`t_q=0`$, that is,

$$\stackrel{~}{\mathrm{tr}L(\stackrel{}{x})}_{\text{CE}}=\frac{\frac{1}{3}_\varphi 𝒟UdetD[U,\varphi ]\mathrm{e}^{S_\text{G}[U]}\mathrm{tr}L(\stackrel{}{x})}{\stackrel{~}{Z}_{\text{CE}}(t_q=0)}.$$
(13)
In the same way as in the case of $`\stackrel{~}{𝒫}_𝒩`$ a similar variation is also possible for $`\stackrel{~}{\mathrm{tr}L(\stackrel{}{x})}_{\text{CE}}`$, i.e., we can consider $`\mathrm{tr}L(\stackrel{}{x})_{\text{CE}}`$ with $`n_q=0`$ for example as

$$\mathrm{tr}L(\stackrel{}{x})_{\text{CE}}=\frac{_0^{2\pi }\frac{\mathrm{d}\varphi }{2\pi }𝒟UdetD[U,\varphi ]\mathrm{e}^{S_\text{G}[U]}\mathrm{tr}L(\stackrel{}{x})}{Z_{\text{CE}}(n_q=0)},$$
(14)
where $`L(\stackrel{}{x})`$ stands for the Polyakov loop defined by

$$L(\stackrel{}{x})=𝒯\mathrm{exp}\left[\mathrm{i}g_0^\beta dx_4A_4(x_4)\right].$$
(15)
Here $`𝒯`$ represents the time-ordering and $`A_4(x)`$ is the temporal component of gauge fields.

We can explicitly evaluate the integration with respect to $`\varphi `$ once the actual form of $`detD[U,\varphi ]`$ is given. Now we will exploit the hopping parameter expansion. Following the calculation of Ref. gre84  the Dirac determinant amended by $`\varphi `$ can be readily expanded in the lowest order of the hopping parameter expansion as follows;

$$detD[U,\varphi ]=\mathrm{exp}\left[H\underset{\stackrel{}{x}}{}\left\{\mathrm{e}^{\mathrm{i}\varphi }\mathrm{tr}L(\stackrel{}{x})+\mathrm{e}^{\mathrm{i}\varphi }\mathrm{tr}L^{}(\stackrel{}{x})\right\}\right],$$
(16)
where $`H=2(2\kappa )^{N_\tau }`$ and $`L(\stackrel{}{x})`$ here is the Polyakov loop defined on the lattice by

$$L(\stackrel{}{x})=\underset{x_4=a}{\overset{N_\tau a}{}}U_4(\stackrel{}{x},x_4),$$
(17)
which can be graphically represented like in Fig. 2 (a). With the notations $`M=\mathrm{tr}L(\stackrel{}{x})=|M|\mathrm{e}^{\mathrm{i}\theta _M}`$, the integration in terms of $`\varphi `$ in Eq. (2) results in mil87

$$_0^{2\pi }\frac{\mathrm{d}\varphi }{2\pi }\mathrm{e}^{\mathrm{i}n_q\varphi }\mathrm{e}^{H(\mathrm{e}^{\mathrm{i}\varphi }M+\mathrm{e}^{\mathrm{i}\varphi }M^{})}=\mathrm{e}^{\mathrm{i}n_q\theta _M}I_{n_q}(2H|M|),$$
(18)
where $`I_\nu (x)`$ is the modified Bessel function of the first kind. By differentiating the integrand with respect to $`\varphi `$, we can have the expression for the expectation value of the quark number as

$$n_q=HM_{n_q1}HM^{}_{n_q+1}.$$
(19)
This relation can be confirmed analytically also by the formula $`I_{\nu 1}(x)I_{\nu +1}(x)=(2\nu /x)I_\nu (x)`$. Also it can be expressed as $`n_q=2\mathrm{i}H\mathrm{Im}\mathrm{tr}L(\stackrel{}{x})\mathrm{e}^{\mathrm{i}\varphi }_{n_q}`$. This form of the relation has been found already in Ref. mil87 . It would be worth mentioning that a different but similar relation has been discussed also in a numerical way in Ref. for00 , which suggests that the above relation could persist beyound the leading order of the hopping parameter expansion.

In the lowest order of the hopping parameter expansion, we can write down DeTar–McLerran’s order parameter as

$$𝒫_𝒩=\frac{𝒟U\mathrm{e}^{S_\text{G}[U]+\mathrm{i}𝒩\theta _M}I_𝒩(2H|M|)}{𝒟U\mathrm{e}^{S_\text{G}[U]+2H|M|\mathrm{cos}\theta _M}}.$$
(20)
As for the expectation value of the Polyakov loop in the canonical ensemble, we can immediately write down as

$$\mathrm{tr}L_{\text{CE}}=\frac{1}{V}\frac{𝒟U\mathrm{e}^{S_\text{G}[U]+\mathrm{i}\theta _M}I_0(2H|M|)|M|}{𝒟U\mathrm{e}^{S_\text{G}[U]}I_0(2H|M|)}.$$
(21)
In the next section, we address the problem of taking the thermodynamic limit for these would-be order parameters (20) and (21). The point will turn out that the sum of the Polyakov loop $`M=\mathrm{tr}L(\stackrel{}{x})`$ is an extensive quantity and as large as proportional to the volume.

## III THERMODYNAMIC LIMIT

### III.1 Failure as order parameters

We consider the thermodynamic limit here, that is, we make the system volume go to infinity. Then we will find that both $`𝒫_𝒩`$ and $`\mathrm{tr}L_{\text{CE}}`$ fail serving as order parameters for the deconfinement transition. Roughly speaking, the failure as order parameters stems from fractional quark excitations permitted even in the canonical ensemble under the infinite volume limit. As schematically depicted in Fig. 3, a pair of a quark and an antiquark can be excited (drawn by the thin curves) because such an excitation is mesonic and the net quark number is zero. When the distance between source quarks (drawn by the thick curves) gets large, it would be regarded as the system composed of two mesonic excitations as shown in the left figure of Fig. 3. In the limit of the infinite distance between source quarks, two mesonic excitations should become individual (clustering decomposing property as shown in the right figure of Fig. 3), which means that the Polyakov loop associated with a single quark excitation can be effectively screened by the fractional excitation of dynamical quarks. Although any fractional excitation of dynamical quarks is prohibited at first by the definition of the canonical ensemble, it would be allowed eventually in the infinitely large volume limit due to the clustering decomposing property. As a result, $`𝒫_{𝒩0}`$ amounts to non-zero $`𝒫_{𝒩=0}`$ and $`\mathrm{tr}L_{\text{CE}}`$ is reduced to the Polyakov loop expectation value in the conventional grand canonical ensemble under the thermodynamic limit.

Let us confirm the above intuitive understanding in the actual expressions given by Eqs. (20) and (21). Without any spontaneous breaking of the center symmetry, $`𝒫_𝒩`$ is expected to vanish for $`𝒩0mod3`$. This expected property derives from the integration with respect to $`\theta _M`$ in Eqs. (20) and (21). Owing to the projection by the integration over $`\varphi `$, the partition function depends on $`\theta _M`$ only through the factor $`\mathrm{e}^{\mathrm{i}𝒩\theta _M}`$. Therefore, if $`𝒩`$ takes a non-zero value, the integration with respect to $`\theta _M`$ makes zero as a result of the whole average over the phase factor. Unless external perturbation is introduced, $`𝒫_{𝒩0}=0`$ should be maintained even in the thermodynamic limit. In other words, any spontaneous breaking of the symmetry cannot be described without external perturbation. Thus even after taking the thermodynamic limit, we have vanishing $`𝒫_{𝒩0}`$ in the absence of external perturbation regardless of whether the system lies in the confined or deconfined phase. This is what Meyer-Ortmanns found in her analysis.

In the presence of small perturbation denoted by $`h`$, the phase factor is modified to $`\mathrm{e}^{\mathrm{i}𝒩\theta _M+2h|M|\mathrm{cos}\theta _M}`$. In the thermodynamic limit, $`|M|`$ may become as large as proportional to the volume and then only the stationary point $`\theta _M=0`$ dominates the integration with respect to the Polyakov loop. Consequently $`𝒫_𝒩`$ loses any $`𝒩`$ dependence to lead to $`𝒫_𝒩0`$ for arbitrary $`𝒩`$. In other words, the regained center symmetry is always broken and the partition function is reduced into that in the grand canonical ensemble in which the center symmetry is explicitly broken whenever small perturbation is introduced in order to realize the spontaneous symmetry breaking. Thus the situation is sharply contrast to Meyer-Ortmanns’ analysis. Once external perturbation is applied, the canonical description resolves itself into the grand canonical description. If we deal carefully with the thermodynamic limit and try to retain the canonical description, even the spontaneous breaking of the center symmetry becomes unavailable. In any case, the state of affairs is somewhat of antinomy.

In order to gain a deeper insight, we will deal with a toy model in which the Polyakov loop dynamics is emulated by the spin system in the next subsection.

### III.2 A toy model

In this subsection we investigate the thermodynamic limit of the canonical ensemble by means of the Ising model as a toy tool. The argument of the center symmetry tells us that the Polyakov loop dynamics can be effectively described by the Ising model with the Z(2) symmetry for the SU(2) gauge theory and by the Potts model with the Z(3) symmetry for the SU(3) gauge theory. In fact, such an effective description has been established also in the numerical study sve97 .

For simplicity, we will focus here on the Ising model with the Z(2) symmetry, though the discussion is easily extended to the case of the Potts model. The counterpart for the order parameter (20) can be immediately expressed in terms of the spin variables as

$$𝒫_1^{\text{Ising}}=\frac{_s\mathrm{e}^{J{\scriptscriptstyle s_ns_m}}\frac{1}{2}\left(\mathrm{e}^{H{\scriptscriptstyle s_n}}\mathrm{e}^{H{\scriptscriptstyle s_n}}\right)}{_s\mathrm{e}^{J{\scriptscriptstyle s_ns_m}+H{\scriptscriptstyle s_n}}}$$
(22)
with the exchange interaction $`J`$ between the nearest neighbor sites under a magnetic field $`H`$. Also the counterpart for the order parameter (21) can be written as

$$s_{\text{CE}}^{\text{Ising}}=\frac{1}{V}\frac{_s\mathrm{e}^{J{\scriptscriptstyle s_ns_m}}\frac{1}{2}\left(\mathrm{e}^{H{\scriptscriptstyle s_n}}+\mathrm{e}^{H{\scriptscriptstyle s_n}}\right)s_n}{_s\mathrm{e}^{J{\scriptscriptstyle s_ns_m}}\frac{1}{2}\left(\mathrm{e}^{H{\scriptscriptstyle s_n}}+\mathrm{e}^{H{\scriptscriptstyle s_n}}\right)}.$$
(23)
For the purpose of investigating what happens in the thermodynamic limit, it will be sufficient to consider the simplest case in which the exchange interaction is absent ($`J=0`$). Then the order parameters given by Eqs. (22) and (23) can be exactly evaluated, namely, they are trivially zero. The naive expectation is that these order parameters remain zero while the magnitude of the exchange interaction is small, regardless of the presence of an external magnetic field $`H0`$. They will come to take a finite value when the strength of the interaction grows large enough to bring about the dynamical symmetry breaking. Actually, however, this naive expectation is not realized at all. We must follow the standard procedure to describe the spontaneous symmetry breaking; we add an infinitesimal external field denoted by $`h`$ and then take the thermodynamic limit ($`V\mathrm{}`$) and finally turn off the external field ($`h0`$). The order parameters under an infinitesimal external field $`h`$ is calculated for $`J=0`$ as

$`𝒫_1^{\text{Ising}}`$ $`={\displaystyle \frac{_s\frac{1}{2}\left(\mathrm{e}^{(H+h){\scriptscriptstyle s_n}}\mathrm{e}^{(Hh){\scriptscriptstyle s_n}}\right)}{_s\mathrm{e}^{(H+h){\scriptscriptstyle s_n}}}}`$
$`={\displaystyle \frac{\left(2\mathrm{cosh}(H+h)\right)^V(2\mathrm{cosh}(Hh))^V}{2\left(2\mathrm{cosh}(H+h)\right)^V}}`$
$`{\displaystyle \frac{1}{2}}(V\mathrm{},Hh>0).`$ (24)
We can draw the similar result for another candidate, that is,

$`s_{\text{CE}}^{\text{Ising}}`$ $`={\displaystyle \frac{1}{V}}{\displaystyle \frac{_s\frac{1}{2}\left(\mathrm{e}^{(H+h){\scriptscriptstyle s_n}}+\mathrm{e}^{(Hh){\scriptscriptstyle s_n}}\right)s_n}{_s\frac{1}{2}\left(\mathrm{e}^{(H+h){\scriptscriptstyle s_n}}+\mathrm{e}^{(Hh){\scriptscriptstyle s_n}}\right)}}`$
$`={\displaystyle \frac{\left(2\mathrm{cosh}(H+h)\right)^{V1}2\mathrm{sinh}(H+h)\left(2\mathrm{cosh}(Hh)\right)^{V1}2\mathrm{sinh}(Hh)}{\left(2\mathrm{cosh}(H+h)\right)^V+\left(2\mathrm{cosh}(Hh)\right)^V}}`$
$`\mathrm{tanh}(H+h)(V\mathrm{},Hh>0).`$ (25)
The above results of Eqs. (24) and (25) clearly signify that the would-be order parameters given by Eqs. (22) and (23) are to be disqualified under the thermodynamic limit because they always remain finite. In other words, any expectation value calculated in the canonical ensemble in which the center symmetry is seemingly recovered would be reduced into that calculated in the grand canonical ensemble in which the symmetry is explicitly broken. The point is that the contribution among the projective superposition with respect to $`\varphi `$ survives only when it is parallel to the direction of the external magnetic field. As a result, the projection into the canonical ensemble is diminished only to result in the same description as in the grand canonical ensemble. We can intuitively understand the situation in the following way: In the thermodynamic limit, the projective superposition in Eq. (23) for example becomes

$$\mathrm{e}^{H{\scriptscriptstyle s_n}}+\mathrm{e}^{H{\scriptscriptstyle s_n}}\mathrm{e}^{H|{\scriptscriptstyle s_n}|},$$
(26)
which follows that the spin configurations for which $`\mathrm{e}^{H{\scriptscriptstyle s_n}}`$ dominates should be definitely separated from the spin configuration for which $`\mathrm{e}^{H{\scriptscriptstyle s_n}}`$ does, that is, the system becomes non-ergodic. Which term survives in the thermodynamic limit depends on which direction an external field is introduced in.

The same mechanism makes the canonical ensemble in QCD be unstable under the thermodynamic limit. For the configuration of the Polyakov loop with a macroscopic order, $`M=\mathrm{tr}L(\stackrel{}{x})=|M|\mathrm{e}^{\mathrm{i}\theta _M}\mathrm{O}(V)\mathrm{}`$, the integration with respect to $`\varphi `$ in Eq. (18) becomes

$`{\displaystyle _0^{2\pi }}{\displaystyle \frac{\mathrm{d}\varphi }{2\pi }}\mathrm{e}^{\mathrm{i}n_q\varphi }\mathrm{e}^{H(\mathrm{e}^{\mathrm{i}\varphi }M+\mathrm{e}^{\mathrm{i}\varphi }M^{})}={\displaystyle _0^{2\pi }}{\displaystyle \frac{\mathrm{d}\varphi }{2\pi }}\mathrm{e}^{\mathrm{i}n_q\varphi }\mathrm{e}^{2H|M|\mathrm{cos}(\varphi +\theta _M)}`$
$``$ $`{\displaystyle \frac{\sqrt{2\pi }\mathrm{e}^{\mathrm{i}n_q\varphi _0}\mathrm{e}^{2H|M|\mathrm{cos}(\varphi _0+\theta _M)}}{2\pi \left|2H|M|\mathrm{cos}(\varphi _0+\theta _M)\right|^{\frac{1}{2}}}}|_{\varphi _0=\theta _M}(|M|\mathrm{O}(V)\mathrm{}).`$ (27)
In the last line of the above equations only the leading order of the saddle-point approximation ($`2H\mathrm{cos}(\varphi _0+\theta _M)`$ is maximized at $`\varphi _0=\theta _M`$) is left, which becomes exact when $`|M|`$ goes to infinity. The meaning of the saddle-point approximation in Eq. (27) is absolutely the same as what Eq. (26) means. Then infinitesimal perturbation like the $`h`$ field in Eqs. (24) and (25) forces only specific $`\theta _M`$ to be favored in the thermodynamic limit. Thus neither DeTar–McLerran’s order parameter nor the Polyakov loop expectation value in the canonical ensemble would serve as a proper order parameter for the deconfinement transition.

## IV IDEA TO OVERCOME THE PROBLEM

### IV.1 Prospect as order parameters

So far we have clarified why DeTar–McLerran’s order parameter does not work as expected under the thermodynamic limit. Nevertheless, we have prospects of making use of DeTar–McLerran’s order parameter to identify the deconfinement transition. The idea is quite simple. The spin-system counterpart for DeTar–McLerran’s order parameter given in Eq. (22) can be equivalently expressed in the following form:

$$𝒫_1^{\text{Ising}}=\frac{\frac{1}{2}\mathrm{e}^{H{\scriptscriptstyle s_n}}\mathrm{e}^{H{\scriptscriptstyle s_n}}_J}{\mathrm{e}^{H{\scriptscriptstyle s_n}}_J},$$
(28)
where $`\mathrm{}_J`$ denotes the ensemble average with the action $`Js_ns_m`$. Let us imagine that we compute the order parameter of Eq. (28) by the method of the Monte-Carlo simulation. First we generate the ensemble of configurations with the probability specified by the action $`Js_ns_m`$, which has both the symmetric phase with random alignments and the broken phase with a spontaneous magnetization. Next we calculate the ensemble average with those generated configurations. In the symmetric phase, the numerator of Eq. (28) should vanish, while it take a finite value in the broken phase. Hence it is apparent that the would-be order parameter given in Eq. (28) will serve as a proper order parameter. If we prefer the spin expectation value to DeTar–McLerran’s formulation, it would be enough to rewrite Eq. (23) into

$$s_{\text{CE}}^{\text{Ising}}=\frac{1}{V}\frac{\left(\mathrm{e}^{H{\scriptscriptstyle s_n}}+\mathrm{e}^{H{\scriptscriptstyle s_n}}\right)s_n_J}{\mathrm{e}^{H{\scriptscriptstyle s_n}}+\mathrm{e}^{H{\scriptscriptstyle s_n}}_J}.$$
(29)
What is important here is that the probability for configurations is specified by the symmetric part of the action. If we adopt the whole action, for example $`Js_ns_m+\mathrm{ln}\mathrm{cosh}(Hs_n)`$ for Eq. (23), to generate the configurations, the ensemble of configurations becomes non-ergodic in the thermodynamic limit and the would-be order parameters should inevitably fail as demonstrated in Eqs. (24) and (25). In other words, the ensemble of configurations given by $`Js_ns_m+\mathrm{ln}\mathrm{cosh}(Hs_n)`$ is unstable against configurations with $`s_n\mathrm{O}(V)`$ because of the long-ranged nature of the interaction term $`\mathrm{ln}\mathrm{cosh}(Hs_n)`$. As a result the spin configurations with $`s_n>0`$ is completely decoupled from those with $`s_n<0`$ and the canonical description is reduced to the grand canonical one.

The above idea is quite simple in itself. In the rest of this subsection, we will discuss the physical meaning of the above prescription to generate the ensemble of configurations with the probability specified by the symmetric part of the action.

Looking back at the argument around Eq. (27), let us reconsider the validity for the saddle-point approximation which is expected to be exact in the thermodynamic limit. Once $`|M|`$ is really as large as of order $`\mathrm{O}(V)`$, the argument given in Eq. (27) has no suspicious controversy at all. However, going back to the second line of Eq. (2), we realize that the coefficient in front of $`\varphi `$ is at most $`n_q`$ which is now of order $`\mathrm{O}(1)`$. If we considered the thermodynamic limit where we take $`n_q\mathrm{O}(V)\mathrm{}`$ with the number density $`n_q/V`$ fixed, the aforementioned argument of the saddle-point approximation would become exact and the calculation of the canonical ensemble is absolutely identical as that in the grand canonical one, as closely discussed in Ref. hag85 . Then the expectation value of the phase of the Polyakov loop is finite and almost proportional to the quark number density for00 , as seen in Eq. (19).

As we have already emphasized, the canonical ensemble for the purpose of characterizing the deconfinement transition is different. What we should be confronted with is the ensemble where $`n_q`$ is fixed at so small number of order $`\mathrm{O}(1)`$ that the saddle-point approximation would be inapplicable. So as to realize this situation an additional condition is needed; $`|M|=|\mathrm{tr}L(\stackrel{}{x})|`$ must be kept of order $`\mathrm{O}(1)`$. This condition can be achieved in the form of the ensemble average in Eqs. (28) and (29) as numerically demonstrated in the next subsection.

We shall summarize our findings and assertions here again.

* Even when the quark number is fixed in the canonical formulation, it is the quark number density that is actually fixed under the thermodynamic limit because the system becomes unstable against configurations with a macroscopic order. Thus any fluctuation of the quark number of order $`\mathrm{O}(1)`$ is allowed even in the canonical description, though it is unintended. As a result the notion of confinement becomes obscured.
* The idea is simple: If we completely get rid of the fluctuation of the quark number, which is negligible in the thermodynamic limit but responsible for screening the Polyakov loop, the canonical ensemble with respect to the quark number is expected to recover its meaning. Then DeTar–McLerran’s order parameter and also the Polyakov loop should serve as order parameters for the deconfinement transition.
* The fluctuation of the quark number is induced by the infinitely non-local interaction arising from the projective superposition. Such fluctuation can be excluded by taking the ensemble average with configurations generated by the symmetric part of the action. The asymmetric part of the action is regarded as included in the operator part to be averaged over.

### IV.2 Numerical tests

By using the Monte-Carlo simulation, we numerically calculate Eq. (29) which can be comparable with an ordinary magnetization in the presence of an external magnetic field. The lattice size is chosen as $`100\times 100`$ and we adopt the standard Metropolis algorithm. The strengths of the exchange interaction and the magnetic field are $`J=0.5/T`$ and $`H=0.1/T`$ respectively, where $`T`$ is the temperature. The thermalization of the spin system is well achieved after 10000 times sweeps, as can be confirmed by the magnetization distributions shown in Fig. 4 for $`T=0.5`$ and $`T=2.0`$. It should be noted that the magnetization is typically of order $`\mathrm{O}(V)`$ in the ordered phase ($`T=0.5`$) and of order $`\mathrm{O}(1)`$ in the disordered phase ($`T=2.0`$).

The resultant $`s_{\text{CE}}^{\text{Ising}}`$ is presented in Fig. 5. Filled circles stand for the spontaneous magnetization in the case of $`H=0`$ and indicate that the spin system goes through the second-order phase transition at $`T_\text{c}=2.27J=1.13`$ (theoretical value). Blank circles represent the magnetization in the presence of the magnetic field $`H0`$ and show a smooth crossover as the temperature raises. We plot $`s_{\text{CE}}^{\text{Ising}}`$ by crosses with error-bars. Although we take 10000 samples by each 10 intervals, the statistical error is still large because the term $`\mathrm{e}^{H{\scriptscriptstyle s_n}}+\mathrm{e}^{H{\scriptscriptstyle s_n}}`$ in both the numerator and the denominator of Eq. (29) can become so large that the operator to be taken the ensemble average of must furiously fluctuate to lead to large statistical errors. In a sense we can say that this problem of slow convergence arises from the same origin as that of the sign problem for simulations at finite density, in particular, the situation is alike the Glasgow method in lattice QCD at finite density. This parallelism to the sign problem will be more apparent for the order parameter given by Eq. (28) or Eq. (20). As discussed in Sec. III the canonical partition function can vanish due to the average over the phase factor coming from the Dirac determinant.

The results shown in Fig. 5 seem to suggest that the would-be order parameter $`s_{\text{CE}}^{\text{Ising}}`$ should work as expected. Unfortunately, at this stage, we cannot address any stronger statement due to large statistical errors. In order to improve the accuracy, some essential ingenuity such as a reweighting method fod01  are needed rather than simple accumulation of many samples. That is beyond the scope of the present paper.

### IV.3 Construction of the order parameter in QCD

It is easy to write down the QCD counterparts for the order parameters given by Eqs. (28) and (29) in the spin system. Here we present the expressions only for DeTar–McLerran’s order parameter defined by Eq. (11) and the Polyakov loop expectation value given by Eq. (13). The formulation can be straightforwardly extended to general cases.

The quark contribution from the Dirac determinant can be expressed as the effective action in the following form;

$$detD[U]=\mathrm{e}^{S_{\text{eff}}[U,U_4(x_4=N_\tau )]},$$
(30)
where the link variable $`U_4(x_4=N_\tau )`$ receiving a modification under the center transformation is explicitly written for convenience. We can decompose the effective action according to the triality, that is,

$$S_{\text{eff}}^{(k)}[U]=\frac{1}{3}\underset{n=0,\pm 1}{}\mathrm{e}^{\mathrm{i2}\pi kn/3}S_{\text{eff}}[U,z_0^nU_4(x_4=N_\tau )](k=0,\pm 1)$$
(31)
with $`z_0=\mathrm{e}^{\mathrm{i2}\pi /3}\mathrm{Z}(3)`$. Among the decomposed parts, all the contributions with non-trivial triality ($`k=\pm 1`$) are bound to be included in the operator to be averaged over. Therefore the would-be order parameter given by Eq. (11) can be written for $`𝒩=1`$ as

$$\stackrel{~}{𝒫}_1=\frac{\frac{1}{3}\mathrm{e}^{S_{\text{eff}}^{(1)}S_{\text{eff}}^{(1)}}+z_0^1\mathrm{e}^{z_0S_{\text{eff}}^{(1)}z_0^1S_{\text{eff}}^{(1)}}+z_0\mathrm{e}^{z_0^1S_{\text{eff}}^{(1)}z_0S_{\text{eff}}^{(1)}}_{\text{sym}}}{\mathrm{e}^{S_{\text{eff}}^{(1)}S_{\text{eff}}^{(1)}}_{\text{sym}}}.$$
(32)
The ensemble average is to be taken with the probability weight specified by the symmetric part of the action, i.e., $`S_\text{G}[U]S_{\text{eff}}^{(0)}[U]`$. In the same way, the expression for Eq. (13) can be written as

$$\stackrel{~}{\mathrm{tr}L(\stackrel{}{x})}_{\text{CE}}=\frac{\left(\mathrm{e}^{S_{\text{eff}}^{(1)}S_{\text{eff}}^{(1)}}+\mathrm{e}^{z_0S_{\text{eff}}^{(1)}z_0^1S_{\text{eff}}^{(1)}}+\mathrm{e}^{z_0^1S_{\text{eff}}^{(1)}z_0S_{\text{eff}}^{(1)}}\right)\mathrm{tr}L(\stackrel{}{x})_{\text{sym}}}{\mathrm{e}^{S_{\text{eff}}^{(1)}S_{\text{eff}}^{(1)}}+\mathrm{e}^{z_0S_{\text{eff}}^{(1)}z_0^1S_{\text{eff}}^{(1)}}+\mathrm{e}^{z_0^1S_{\text{eff}}^{(1)}z_0S_{\text{eff}}^{(1)}}_{\text{sym}}}.$$
(33)
## V CONCLUDING REMARKS

In order to identify the deconfinement phase transition we employed and elaborated the idea of the canonical ensemble with respect to the quark number or triality. We clarify why the canonical ensemble would be reduced into the grand canonical ensemble eventually in the thermodynamic limit. In order to overcome the problem of the thermodynamic limit, we propose a prescription to compute the ensemble average. The idea is tested by means of an effective model in terms of spin variables in the presence of the magnetic field. The results seem prosperous but turn out to suffer from the severe sign problem. Although our definition of the order parameter for the deconfinement transition can be applied to the Monte-Carlo simulation on the lattice, it is inevitably necessary to get over the sign problem inherent in the lattice QCD at finite density. We find it interesting that two serious questions – one is the criterion for the deconfinement transition at finite temperature, and the other is the sign problem in the lattice QCD at finite density – are closely related to each others.

###### Acknowledgements.

The author, who is supported by Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists, thanks Professor J. Polonyi for stimulating discussions.


format: text

text:
# Determining the Saddle Point in Micromagnetic Models of Magnetization Switching

## Abstract

A numerical model of single-domain nanoscale iron magnets fabricated using scanning-microscope-assisted chemical vapor deposition is simulated using finite-temperature micromagnetics. A Projective-dynamics method is used to determine the magnetization at the saddle point as a function of temperature. This magnetization is found to decrease linearly as the temperature is raised.

An interesting problem in nonequilibrium statistical mechanics, with numerous applications in condensed-matter physics and materials science, is how a system approaches the global minimum of the free energy. One simple, but common, situation is that of a local free-energy minimum which is separated from the global minimum by a free-energy barrier. In other words, any path between the metastable local minimum and the equilibrium global minimum involves an initial increase in the free energy. The most probable path between the minima involves the minimum free energy increase since spontaneous increases in free energy are improbable. The maximum free energy along that most probable path corresponds to the saddle point. Often, the behavior of the nonequilibrium dynamics are dominated by properties near the minima and the saddle point. While the minima can be determined by examining histograms of the state of the system, the saddle point is much harder to determine. In this paper we present a technologically important example of finding the saddle point using a projective-dynamics technique.

The results described here are for a numerical model of single-domain nanoscale iron magnets that have been fabricated using chemical vapor deposition directed by a scanning-tunneling microscope . These magnetic pillars have cross-sectional dimensions on the order of $`10`$ nm and extend on the order of $`100`$ nm perpendicular to the substrate. The magnetic pillars are modeled by a one-dimensional array of magnetization density vectors $`\stackrel{}{M}(\stackrel{}{r})`$ with fixed length $`M_s`$. Each vector precesses around a local field $`\stackrel{}{H}(\stackrel{}{r})`$ according to the Landau-Lifshitz-Gilbert equation

$$\frac{d\stackrel{}{M}(\stackrel{}{r})}{dt}=\frac{\gamma _0}{1+\alpha ^2}\stackrel{}{M}(\stackrel{}{r})\times \left[\stackrel{}{H}(\stackrel{}{r})\frac{\alpha }{M_s}\stackrel{}{M}(\stackrel{}{r})\times \stackrel{}{H}(\stackrel{}{r})\right],$$
(1)
where the electron gyromagnetic ratio $`\gamma _0=1.76\times 10^7`$ Hz/Oe. The phenomenological damping parameter $`\alpha =0.1`$ was chosen to give underdamped dynamics. Properties of bulk iron were assumed with the saturation magnetization $`M_s=1700`$ emu/cm3 and the exchange length $`\mathrm{}_x=2.6`$ nm. The model pillar has a square cross section with area $`4\mathrm{}_x^2,`$ and is $`34\mathrm{}_x`$ long. Details of the numerical approach, including the inclusion of thermal fluctuations in the local field, are given in Ref. .

These nanopillars have a strong shape anisotropy, and the magnetization is most favorably oriented along the long axis of the magnet, which is taken to be the $`z`$ direction. The nanomagnets are prepared in a metastable state via equilibration in an externally-applied magnetic field $`+H_0\widehat{z}`$, which is subsequently varied rapidly in magnitude (without changing its orientation) to $`H_0\widehat{z}.`$ If the field $`H_0`$ is less than the coercive field, here $`H_c1500`$ Oe , a free-energy barrier separates the positively-oriented magnetization from the equilibrium negative orientation. Previous simulations  suggest that the saddle point in these pillars corresponds to the nucleation of a region of reversed magnetization at one of the ends of the pillar. The nucleated region grows until the entire magnetization is reversed.

A useful technique for determining the magnetization corresponding to the free-energy maximum of the barrier, also called the saddle point, is the projective-dynamics technique . The essence of this technique is projecting the original dynamics described in a high-dimensional phase space onto a probabilistic dynamic in a phase space of much lower dimension. A specific example is helpful. The dynamics of an Ising model can be projected onto the number of overturned spins, a measure of the global magnetization. The dynamics of the one-dimensional model are then described in terms of the probability that the number of overturned spins increases or decreases. Thinking in terms of the number of spins in the stable orientation, these probabilities are referred to as the probability of growing ($`P_{grow}`$) and shrinking ($`P_{shrink}`$), respectively .

Projective dynamics for the model nanopillars is more complicated than that for the Ising model for two reasons. First, the magnetization is a continuous, three-dimensional variable. This is handled by projecting the $`z`$-component of the global magnetization into uniformly sized bins . Then $`P_{grow}`$, the probability of the region of stable magnetization “growing,” is the probability that during an integration step the magnetization moves into a bin corresponding to a smaller magnetization. Similarly, $`P_{shrink}`$ is the probability that the magnetization moves into a bin corresponding to a larger magnetization. The second complication involves the persistent, subcritical regions of reduced magnetization that develop at each end of the pillar. To accomodate this, during the projective-dynamics analysis the pillar is divided into its top and bottom parts with the normalized, “global” magnetization $`M_z`$ of each half considered separately.

The measured $`P_{grow}`$ and $`P_{shrink}`$ at $`H_0`$$`=`$$`1000`$ Oe are shown in Fig. 1 for temperatures of $`20`$ K, $`50`$ K, and $`100`$ K. For each temperature the results represent averaging over more than $`10^9`$ integration steps and on the order of $`10^4`$ switches. Consider the probabilistic dynamics at $`M_z`$$`=`$$`0.85`$ for $`T`$$`=`$$`100`$ K. Here $`P_{shrink}`$$`>`$$`P_{grow}`$ and on average the magnetization moves to the right. For $`M_z`$$`=`$$`0.98`$ the situation is reversed, and on average the magnetization moves to the left. The crossing point where $`P_{shrink}`$$`=`$$`P_{grow}`$ corresponds to a locally stable fixed point since, on average, here the magnetization moves towards the crossing. In this case the crossing is the metastable local free-energy minimum. For the crossing near $`M_z`$$`=`$$`0.74`$ the magnetization moves away from the point on average, and it corresponds to an unstable local maximum in the free energy. This is the saddle point. Measuring the value of $`M_z`$ for which $`P_{shrink}`$$`=`$$`P_{grow}`$ is a convenient method for finding the saddle point.

From Fig. 1 it is obvious that $`M_z`$ at the saddle point is different for magnetization reversal at different temperatures. The value of $`M_z`$ at the saddle point is presented vs temperature in Fig. 2. The point where $`P_{shrink}`$$`=`$$`P_{grow}`$ is estimated using lines fit via a least-squares method to each probability in the region near the crossing. The error bars are taken large enough to include nearly all $`M_z`$ for which the difference in the probabilities is less than the fluctuations in their estimates.

There is a clear linear trend in $`M_z`$ at the saddle point with respect to the temperature. The solid line in Fig. 2 is an unweighted least-squares fit to the data. This fit indicates that $`M_z(T0)`$$`=`$$`0.860`$ at the saddle point. The dependence on temperature is quite strong, since $`M_z(100`$ K$`)`$$`=`$$`0.72`$.

In summary, micromagnetic simulations at finite temperature have been used to investigate magnetization reversal in single-domain nanoscale magnets. The saddle point in the free energy was determined using the projective-dynamics technique, and the magnetization at the saddle point was found to depend linearly on the temperature for all the temperatures investigated, $`T`$$``$$`100`$ K.

This work was supported in part by the Laboratory Directed Research and Development Program of Oak Ridge National Laboratory, managed by UT-Battelle, LLC for the U.S. Department of Energy under Contract No. DE-AC05-00OR22725, and by the Ames Laboratory, which is operated for the U.S. Department of Energy by Iowa State University under Contract No. W-7405-82. This work was also supported by NSF grant No. DMR-0120310 and by FSU/CSIT. Computer resources were provided by the FSU Deptartment of Physics.


text:
# Breathers in a Discrete Nonlinear Schrödinger type Model : Exact Stability Results.

## Abstract

Exact stability analysis of 1-site breathers in an NLS type model (ref. , see below) indicates destabilisation through a growth rate becoming positive as a stability border is crossed, while above a critical spatial decay rate ($`\lambda _c`$) the breather is always unstable. Similar analysis on 2-site breathers with phase difference $`\pi `$ (see below) indicates destabilisation through Krein collision of a pair of growth rates. All eigenmodes can be exactly calculated and conform to numerical estimates.
PACS number(s): 05.45.-a, 63.20.Pw, 63.20.Ry

The question of dynamical stability of discrete breathers has been addressed in the literature and a number of useful results are now known. The stability of a monochromatic breather depends on the eigenvalues of a symplectic matrix, the so called Floquet matrix. Aubry  introduced a novel extension to the concept of the Floquet matrix , thereby greatly simplifying the task of computing the Floquet multipliers and of looking for the onset of instability. In particular, he established that, one-site breathers sufficiently close to the ‘anti-integrable’ limit are dynamically stable. However, few exact results relating to the multipliers and the associated eigenmodes are known (, ; see references therein for a background to the problem) notably because the breather solutions around which the linearisation is to be performed are themselves obtained numerically, even though with high precision. In this context it is useful to have models that admit of exact breather solutions as also exact stability results. In  we introduced a number of piecewise smooth (PWS) models including a nonlinear Schrödinger type one and worked out exact one-site breather solutions where one of the sites is excited to a value above a threshold while all other sites remain below the threshold ( we denote these as respectively high and low). In the present paper we work out exact stability results for these breathers as also for 2-site breathers where two sites with an arbitrary intervening gap are high while all other sites are low. While the destabilisation of 1-site breathers always occur through the growth rate of a mode becoming positive, the 2-site breathers get distabilised through a Krein collision of the growth rates (see below).
More precisely, the PWS Schrödinger type model introduced in  is

$`i{\displaystyle \frac{d\psi _n}{dt}}+V(\psi _{n+1}+\psi _{n1})=\gamma \psi _nf(\psi _n),(n=0,\pm 1,\mathrm{})`$ (1)
$`f(z)=(1a/|z|)\mathrm{\Theta }(|z|a),`$ (2)
where $`a(>0),V(>0),\gamma `$ are appropriate constants and $`\mathrm{\Theta }`$ denotes the Heaviside step function. We scale the nearest neighbour coupling strength and the threshold parameter to $`V=1`$, $`a=1`$. Fig.1 compares $`f(z)`$ of (2) with the more usual form $`f(z)=|z|^2`$ of the discrete non-linear Schrödinger equation. The one-site breather solution centred at $`n=0`$ is then

$`\overline{\psi }_n`$ $`=`$ $`\overline{\varphi }_ne^{i\omega t},`$ (3)
$`\overline{\varphi }_n`$ $`=`$ $`{\displaystyle \frac{\gamma \lambda ^{|n|}}{\gamma +\lambda 1/\lambda }}.`$ (4)
The spatial decay rate $`\lambda (|\lambda |<1)`$ and the breather frequency $`\omega `$ satisfy

$$|\omega |>2,\omega =(\lambda +1/\lambda ),\omega \gamma <0,|\gamma |>1+\frac{1}{|\lambda |},$$
(5)
(note: Equations (10), (11) in  expressing the above relations contained errors due to oversight, which we regret).We distinguish between breathers of type A $`(\lambda >0,\omega <2)`$ and type B $`(\lambda <0,\omega >2)`$), referring to the former in presenting our results below (results for type B are obtained analogously). The overbar in (3),(4) is used to distinguish the breather solution from neighbouring perturbed solutions (see below).
We consider these perturbations in the ‘rotating frame’ in the form

$$\psi _n=(\overline{\varphi }_n+u_n(t))e^{i\omega t}.$$
(6)
Spliting $`u_n`$ into real and imaginary parts, $`u_n=x_n+iy_n`$, one observes that the time variations of $`x_n,y_n`$ are obtained from the Hamiltonian

$`H={\displaystyle \underset{n}{}}\{(x_{n+1}x_n+y_{n+1}y_n)+{\displaystyle \frac{\omega }{2}}(x_{n}^{}{}_{}{}^{2}+y_{n}^{}{}_{}{}^{2})\}+`$ (7)
$`+{\displaystyle \frac{\gamma }{2}}\{x_{0}^{}{}_{}{}^{2}+y_{0}^{}{}_{}{}^{2}(1{\displaystyle \frac{1}{b}})\},`$ (8)
where $`b=\frac{\gamma }{\gamma +\lambda \frac{1}{\lambda }}`$. Eigenvalues of the matrix $`L`$ determining the evolution of the $`x_n`$’s and $`y_n`$’s will be termed growth rates (we denote a typical growth rate by the symbol $`p`$) while a more convenient approach in the present context is to eleminate the $`y_n`$’s or $`x_n`$’s, arriving at

$$\ddot{X}=AX,\ddot{Y}=A^TY(\mu 0),$$
(9)
where we define $`X=(\mathrm{}\mathrm{}..x_2,x_1,x_0,x_1,x_2,\mathrm{}\mathrm{}.)^T`$, $`Y=(\mathrm{}\mathrm{}\mathrm{}y_2,y_1,y_0,y_1,y_2,\mathrm{}\mathrm{}..)^T`$, and $`A`$ is a banded matrix with elements $`a_{mn}(\mathrm{}<m,n<\mathrm{})`$ differing from zero only when $`|mn|2`$ :

$`a_{0,0}=2(\gamma \lambda 2\lambda ^2),`$
$`a_{0,1}=a_{0,1}=3\lambda +1/\lambda ,`$
$`a_{1,0}=a_{1,0}=2\lambda +2/\lambda \gamma ,`$
$`a_{n,n}=(\lambda ^2+1/\lambda ^2+4),(n0),`$
$`a_{n,n+1}=a_{n+1,n}=2(\lambda +1/\lambda ),(n0,1),`$
$`a_{n,n+2}=a_{n+2,n}=1`$.
The eigenvalues $`(\mu )`$ of $`A`$ are related to growth rates $`(p)`$ as $`\mu =p^2`$. The spectrum of eigenvalues $`\mu `$ and the corresponding eigenmodes can all be worked out exactly for any given parameter values $`\lambda ,\gamma `$ satisfying (5). There exists a band of eigenvalues $`(\omega 2)^2<\mu <(\omega +2)^2`$ with extended eigenmodes, there being a symmetric and an antisymmetric eigenmode associated with each eigenvalue. Additionally, there exists a single localised symmetric eigenmode for $`\gamma >\gamma _1(\lambda )`$ where

$`\gamma _1(\lambda )={\displaystyle \frac{(\lambda \frac{1}{\lambda })\sqrt{((\lambda +\frac{1}{\lambda })^22(\lambda +\frac{1}{\lambda }))}}{(\lambda \frac{1}{\lambda })+\sqrt{((\lambda +\frac{1}{\lambda })^22(\lambda +\frac{1}{\lambda })}}}.`$ (10)
For any given $`\lambda ,\gamma _1`$ is thus the threshold for the localised mode: for $`\gamma <\gamma _1`$ there is no localised mode and the band edges are empty, while for $`\gamma =\gamma _1`$ an eigenmode appears at the band-edge, and dissociates from the band as the localised mode for $`\gamma >\gamma _1`$. One notes in this context that there always exists a trivial eigenmode $`\mu =0`$ with $`x_n=0,y_n\overline{\varphi }_n`$, which is to be distinguished from the eigenvectors of $`A,A^T`$ with $`\mu =0`$, the latter corresponding to a spurious mode in the present context. While a typical localised eigenvector of $`A`$ is of the form

$`x_n=\alpha \chi _1^{|n|}+\beta \chi _2^{|n|},`$ (11)
where the spatial decay rates $`\chi _1,\chi _2`$ (each of modulus $`<1`$) as also the ratio $`\alpha /\beta `$ can be obtained exactly, the eigenvector corresponding to $`\mu =0`$ looks like

$`x_n=(\alpha +\beta |n|)\lambda ^{|n|}.`$ (12)
The growth rates for the localised mode are imaginary $`(p=\pm i\sqrt{|}\mu |)`$ so long as $`\gamma _1(\lambda )<\gamma <\overline{\gamma }(\lambda )`$ where the transition value $`\overline{\gamma }`$ can also be obtained exactly in the model:

$`\overline{\gamma }(\lambda )={\displaystyle \frac{(1+4\lambda ^2+\lambda ^4)(1\lambda ^2)}{2\lambda ^3}}.`$ (13)
As $`\gamma `$ crosses $`\overline{\gamma }`$ from below, the isolated eigenvalue $`\mu `$ crosses $`\mu =0`$ and becomes positive, giving rise to a positive growth rate $`(p_+=\sqrt{\mu };`$ this is associated with a negetive growth rate $`p_{}=\sqrt{\mu }`$, making up a pair) and implying destabilisation of the breather.

Fig.2 depicts the disposition of eigenvalues on either side of the destabilisation, while

Fig.3 indicates how the disposition of growth rates in the complex plane changes across the destabilisation border.

Fig.4 depicts the $`\lambda \gamma `$ parameter space for the model together with the curves $`C_1(\gamma =\gamma _0(\lambda )1+\frac{1}{\lambda }),C_2(\gamma =\overline{\gamma }(\lambda )),C_3(\gamma =\gamma _1(\lambda ))`$. The breather exists only for points lying above $`C_1`$ and is stable only for points lying below $`C_2`$, $`i.e.,`$ for $`(\lambda ,\gamma )`$ lying in between $`C_1`$ and $`C_2`$ the model admits of a stable breather solution. The curve $`C_3`$ gives the threshold for the localised mode. Interestingly, the model predicts the existence of a critical value $`(\lambda _c)`$ of the spatial decay rate corresponding to the intersection of $`C_1`$ and $`C_2`$ $`(\lambda _c0.6948)`$ so that the breathers with $`\lambda >\lambda _c`$ are unstable for all values of the strength of non-linearity $`\gamma (>\gamma _0(\lambda ))`$. All these features resemble those of ‘internal’ or ‘breathing’ modes of DNLS breathers whose relevance has been emphasized in the literature (,). However, the antisymmetric ‘pinning’ mode , known to exist for DNLS breathers, is marginal in the present PWS model $`(\lambda =1,\mu =0)`$ and possesses a linear spatial growth.

We show in Fig.$`5a`$ the temporal evolution according to Eqs. (1), (2) of a profile initially (t=0) coinciding with (3), with $`(\lambda ,\gamma )`$ in the region between $`C_1`$ and $`C_2`$ while Fig.5b shows the similar evolution for a breather with $`\lambda >\lambda _c`$, above $`C_1`$. One notes that the breather in Fig. 5a performs stable oscillations while that in Fig. 5b breaks up within a short time, confirming our results. In this context see also Figs. 4,5 of .
The eigenmodes can be obtained exactly in the model. For instance, Fig.6 shows the localised eigenmode for $`\lambda =0.5,\gamma =6`$ calculated from the theory as compared with the eigenmode obtained numerically by diagonalising the matrix $`A`$.

As we mentioned in , the PWS model (1), (2) is a veritable little ‘laboratory’ yielding exact breather solutions of a wide variety. For instance, one can construct two-site breathers with various possible lengths of intervening gap $`N`$. One other feature of crucial relevance is the relative phase between the two sites. This relates to the fact that (4) admits of an arbitrary phase factor $`e^{i\delta }`$ that becomes relevant in 2-site or multisite breather solutions. Thus, considering a breather with two sites high ($`i.e.,`$ with $`|\psi _n|>1)`$ and all the other sites low $`(|\psi _n|<1)`$ with the high sites having opposite phases ($`\delta =\pi `$) one obtains a monochromatic breather solution of the form (3) with the breather profile $`\overline{\varphi }_n`$ now given by

$`\overline{\varphi }_n={\displaystyle \frac{b}{1\lambda ^N}}(\lambda ^n\lambda ^{Nn})(0nN)`$ (14)
$`=b\lambda ^n(n<0)`$ (15)
$`=b\lambda ^{nN}(n>N),`$ (16)
where $`b=\frac{\gamma (1\lambda ^N)}{\gamma (1\lambda ^N)(\lambda ^1\lambda )}`$ and $`\lambda ,\gamma `$ are to satisfy $`\gamma >\frac{1+\lambda ^1}{1\lambda ^N}.`$ Here the two high sites have been chosen to be at $`n=0`$ and $`n=N`$ for any given separation $`N`$. It is now easy to construct the stability matrix $`A`$ of Eq.(9) and to obtain the spectrum of eigenvalues together with the exact eigenmodes. There is once again a band extending from $`(\omega 2)^2`$ to $`(\omega +2)^2`$ (recall that we are considering only type A breathers, with $`\omega <2`$) containing extended eigenmodes, together with a number of isolated eigenvalues corresponding to localised modes. The latter appear above a certain threshold value of $`\gamma `$ (for any given $`\lambda `$), and are now typically four in number of which one is the trivial eigenvalue $`\mu =0`$ while the others lie in the range $`(\omega +2)^2<\mu <0`$. Ordering these as $`\mu _1<\mu _2<\mu _3<\mu _4=0,`$ the eigenmodes corresponding to $`\mu _1`$ and $`\mu _4`$ are antisymmetric (with respect to the sites $`n=0,N`$) while the other two are symmetric. Fig.7 shows the four eigenmodes for $`N=8,\gamma =4.9443,\lambda =0.5`$. For a given $`N,\lambda `$, as $`\gamma `$ is made to increase from the threshold value, a boundary $`\gamma _H(\lambda )`$ is reached when $`\mu _2`$ and $`\mu _3`$ become coincident and then become complex. In terms of the growth rates $`p`$ $`(p=\pm \sqrt{\mu })`$ there occurs a Krein collision of a quartet of growth rates on the imaginary axis and the subsequent spliting of the growth rates out of the complex plane, as depicted schematically in Fig.8.

In a model with sufficiently smooth $`f(z)`$ (ref. Eq.2) this leads to the Hamiltonian Hopf bifurcation (, see references therein; the Krein collision has been discussed in  in terms of the Floquet multipliers related to the breather solution - the latter are related to our eigenvalues ($`\mu `$) and the corresponding growth rates $`(p)`$ in a simple manner; there exists a large body of literature relating to Krein collision in symplectic mappings - see, $`e.g.`$, ). In the present context, the collision signifies a loss of stability. For instance, with $`N=8,\lambda =0.5,`$ the model predicts a loss of stability of the two-site breather through a Krein collision at $`\gamma _H4.9444`$.

Fig.9 depicts the time evolution of an initial profile chosen in accordance with (16) with $`\gamma `$ on either side of the collision: one notes that the breather indeed remains stable for $`\gamma <\gamma _H`$ and breaks up quickly for $`\gamma >\gamma _H`$.

As mentioned above, in Fig.8 we depict the sequence of changes in the disposition of the growth rates $`p`$ as $`\gamma `$ is made to increase through $`\gamma _H`$ for fixed $`\lambda `$. One observes that at some value $`\stackrel{~}{\gamma }(>\gamma _H)`$ of $`\gamma `$, the root $`\mu `$ becomes positive (this is similar to the crossing of a negative eigenvalue through the value zero at $`\gamma =\overline{\gamma }(\lambda ))`$ with $`\mu _2,\mu _3`$ remaining complex; subsequently at some value $`\gamma =\gamma _H^{}`$ the latter collide at a point on the positive real axis and move off in opposite directions - Krein collision ‘in the reverse’. In the present context of stability of the breathers, however, this reverse collision does not have a direct relevance since the breather, having become unstable at $`\gamma =\gamma _H`$ continues to remain unstable even at $`\gamma =\stackrel{~}{\gamma }`$ and $`\gamma =\gamma _H^{}`$. For large $`N`$, one finds that $`\gamma _H,\stackrel{~}{\gamma }`$ and $`\gamma _H^{}`$ are all contained in a small interval around $`\overline{\gamma }`$ as expected, since the two high sites are then almost unaffected by each other, each essentially behaving like a one-site breather (3), (4).
Interestingly, a 2-site breather with the two high sites in phase with each other can also be constructed in the model, but turns out to be always unstable.
The details of the calculations together with numerical and other results on breather stability will be presented elsewhere.