Simulated annealing (SA) is a kind of relaxation method for finding equilibria of Hamiltonian systems. A set of evolution equations is solved with SA, which is derived from the original Hamiltonian system so that the energy of the system changes monotonically while preserving Casimir invariants inherent to noncanonical Hamiltonian systems. The energy extremum reached by SA is an equilibrium. Since SA searches for an energy extremum, it can also be used for stability analysis when initiated from a state where a perturbation is added to an equilibrium. The procedure of the stability analysis is explained, and some examples are shown. Because the time evolution is computationally time consuming, efficient relaxation is necessary for SA to be practically useful. An acceleration method is developed by introducing time dependence in the symmetric kernel used in the double bracket, which is part of the SA formulation described here. An explicit formulation for low-beta reduced magnetohydrodynamics (MHD) in cylindrical geometry is presented. Since SA for low-beta reduced MHD has two advection fields that relax, it is important to balance the orders of magnitude of these advection fields.

## I. INTRODUCTION

Simulated annealing (SA)^{1} is a kind of relaxation method for obtaining equilibria, stationary states, of Hamiltonian systems. In noncanonical Hamiltonian systems, such as non-dissipative fluid dynamics (see Ref. 2 for review), there exists Casimir invariants that foliate the infinite-dimensional phase space. The surface in the phase space on which each Casimir invariant takes a same value is called a Casimir leaf. A stationary state or an equilibrium of the system is given by an energy extremum on the Casimir leaf.^{3,4}

With SA, a set of artificial evolution equations is solved, which is derived from the original equations of the Hamiltonian system. The set of artificial evolution equations is constructed so that the energy of the system changes monotonically, while preserving the Casimir invariants. The idea of the artificial dynamics was proposed originally for obtaining stationary states of two-dimensional fluid flows^{5–7} and generalized in Ref. 1. In Ref. 1, the formulation of the artificial dynamics was considerably extended and made more workable by providing a means of including additional constraints and by introducing a smoothing kernel. There a variety of stationary states were obtained for two-dimensional vortex dynamics and associated two-layer problems. The terminology “simulated annealing” was introduced in this paper.

The noncanonical Hamiltonian structure of ideal magnetohydrodynamics (MHD) was first given in Ref. 8. For testing the idea of SA, the MHD Hamiltonian structure was first used^{9,10} for low-beta reduced MHD^{11} in a two-dimensional rectangular domain with doubly periodic boundary conditions. Since low-beta reduced MHD has two fields, i.e., vorticity and magnetic flux function, the relaxation path to a stationary state becomes rather complex. However, the application was successful. Later, it was extended to cylindrical geometry, and an equilibrium with magnetic islands was obtained.^{12} Furthermore, SA was applied to high-beta reduced MHD^{13} in toroidal geometry and successfully used to obtain tokamak as well as toroidally averaged stellarator equilibria.^{14}

The above-mentioned results are based on the double bracket formulation of SA.^{15} However, another type of artificial dynamics method for calculating stationary states is metriplectic dynamics.^{15,16} This relaxation method was successfully applied to two-dimensional vortex dynamics and further to axisymmetric MHD equilibria^{17} described by the Grad–Shafranov equation.^{18–20}

From a more general point of view, minimization of a functional like the energy functional for SA in this paper may be done by various methods such as a classical steepest-descent method and its extensions, the Nelder–Mead method,^{21} genetic algorithms,^{22} and many others developed in the context of machine learning. Although they can be used for minimizing the energy functional of the MHD systems, it may not be straightforward to preserve the Casimir invariants since these methods are not based on the Hamiltonian structure of the governing equations. Our method is a natural choice for the energy minimization on a Casimir leaf.

It is pointed out that because SA searches for an energy extremum or an energy minimum in the MHD problems, it can be used not only for equilibrium calculations but also for stability analysis. Suppose an equilibrium is given by solving the Grad–Shafranov equation in toroidal geometry or, e.g., just as a trivial solution in cylindrical geometry. Then, suppose that a perturbation is given to the equilibrium followed by the SA procedure. If the equilibrium is an energy minimum, SA will recover the original equilibrium. If not, SA will lead to another state that is an energy minimum or at least will leave the original equilibrium in a sense of linear stability. Note that it is important to stay on the same Casimir leaf as the original equilibrium when giving the perturbation. Such perturbations were termed dynamically accessible in Refs. 23–25. If a perturbation is not dynamically accessible, SA will not recover the original equilibrium because the perturbed state is on another Casimir leaf. The SA procedure will find another equilibrium that is on the Casimir leaf where the perturbed state exists. In the present paper, the procedure for the stability analysis, especially how the perturbation is given without leaving the Casimir leaf, is developed.

On the basis of the previous studies, one would expect that SA can calculate a wide class of ideal MHD equilibria, including magnetic islands and/or magnetic chaotic regions. Also, SA can be used for analysis of not only linear but also nonlinear stability. However, it is not so straightforward because SA solves an initial-value problem, and the numerical simulation of an initial-value problem generally takes time. Even in the linear stability analysis, it takes time to show that an equilibrium is stable since we have to observe that the given perturbation disappears. Therefore, an acceleration method is necessary for SA to be practically useful. In the present paper, an acceleration method is developed for the double bracket formulation of SA. The method is explicitly described for the low-beta reduced MHD model in cylindrical geometry and applied to the stability analysis.

The paper is organized as follows. In Sec. II, the double bracket formulation of SA for low-beta reduced MHD is introduced. Then, the procedure for the stability analysis is explained in Sec. III. Section IV presents numerical results. Linear stability of a stable equilibrium is analyzed in Sec. IV A, while Sec. IV B is for an unstable equilibrium. An equilibrium, the stability of which is to be subsequently analyzed, is given in Sec. IV A 1, and a dynamically accessible perturbation is given to the equilibrium in Sec. IV A 2. Here, it is explained how the dynamically accessible perturbation stays on the original Casimir leaf. Then, a preliminary result, one before an acceleration method is introduced, is presented, and the problem of slow convergence is explicitly shown in Sec. IV A 3. The slow convergence comes from an imbalance of relaxation speeds between the relaxing fields. In order to confirm this, simulation results where one field is eliminated are presented in Sec. IV A 4. Next, in Sec. IV A 5, the acceleration method is developed that balances the relaxation speeds, and numerical results are shown. Section IV A 6 shows numerical results with another initial condition where the kinetic energy perturbation is larger. As explained earlier, Sec. IV B shows that a given perturbation grows while the energy decreases by SA for an unstable equilibrium. Stability analyses via SA in this paper are examined in Sec. IV C. Section V concludes the paper.

## II. MODEL

### A. Low-beta reduced MHD and normalization

In this paper, the low-beta reduced MHD model^{11} is investigated. We consider a cylindrical plasma with minor radius *a* and length $2\pi R0$. The low-beta reduced MHD equations are

where all physical quantities are normalized by using the length *a*, the magnetic field in the axial direction *B*_{0}, the Alfvén velocity $vA:=B0/\mu 0\rho 0$ with *μ*_{0} and *ρ*_{0} being vacuum permeability and typical mass density, respectively, and the Alfvén time $\tau A:=a/vA$. The cylindrical coordinates $(r,\theta ,z)$ as well as the toroidal angle $\zeta :=z/R0$ are used. The inverse aspect ratio is defined as $\epsilon :=a/R0$. The fluid velocity is $v=z\u0302\xd7\u2207\phi $, and the magnetic field is $B=z\u0302+\u2207\psi \xd7z\u0302$, where the unit vector in the *z* direction is denoted by $z\u0302$. The vorticity is $U:=\u25b3\u22a5\phi $ and the current density is $J:=\u25b3\u22a5\psi $, where the Laplacian in the *r*–*θ* plane is denoted by $\u25b3\u22a5$. The Poisson bracket for two functions *f* and *g* is defined as $[f,g]:=z\u0302\xb7\u2207f\xd7\u2207g$.

The Hamiltonian structure for low-beta reduced MHD, which follows from that of Ref. 8, was first given in Refs. 26 and 27. The Hamiltonian is the summation of kinetic and magnetic energies as

where $D$ is the domain of the cylindrical plasma. The Casimir invariants are

Furthermore,

are also Casimir invariants when only single helicity components are included in the dynamics. Here, *f* and *g* are arbitrary functions of a helical flux $\psi h(x,t):=\psi (x,t)+\epsilon r2/(2qs)$ for a specified safety factor of a rational number $qs=m/n$ with *m* and *n* being a poloidal and a toroidal mode number, respectively.

### B. Double bracket formulation of simulated annealing

We adopt simulated annealing, formulated by using a double bracket.^{1,12,15} For low-beta reduced MHD, the set of artificial evolution equations are conveniently given in the reduced following form:

where the advection fields are defined as

Here, *j *=* *1 or 2, *f* ^{1} is the right-hand side of Eq. (1), *f* ^{2} is the right-hand side of Eq. (2), and $(Kij)$ is a kernel with a definite sign.

Note that the equations have the same form as those of the original low-beta reduced MHD; however, the advection fields are replaced by the artificial ones. The advection fields must be chosen so that the energy of the system changes monotonically; however, there are still a variety of choices through the kernel.

If $(Kij)$ is positive definite, the energy of the system decreases monotonically. On the other hand, Casimir invariants are preserved. In this study, we choose $(Kij)$ to be diagonal with its diagonal components given by

for *i *=* *1 and 2, where $g(x,x\u2032)$ is a Green's function defined through

Here, $\u25b3$ is the Laplacian and $\delta 3(x)$ is the Dirac's delta function in three dimensions, respectively.

## III. STABILITY ANALYSIS PROCEDURE

In this section, the procedure for using SA to assess stability is explained. It consists of three steps.

The first step is to choose an equilibrium. Although SA has been used for obtaining equilibria; here, the equilibrium can be given by any method. For example, any cylindrically symmetric state is a trivial equilibrium in a cylindrical plasma. Although unbeknownst beforehand, the chosen equilibrium can be either stable or unstable. Similarly, any solution to the Grad–Shafranov equation is a candidate equilibrium in axisymmetric toroidal plasma, one for which the SA stability method could be applied. Of course, we could also choose an equilibrium that is obtained by SA.

The second step is to perturb the equilibrium. The key point is to make the perturbation dynamically accessible, i.e., have it stay on the same Casimir leaf as the equilibrium whose stability is being analyzed. This is accomplished by time evolution of *U* and *ψ* under Eqs. (8) and (9) using arbitrary advection fields $\phi \u0303$ and $J\u0303$ that are different from Eqs. (10) and (11). This is possible because the Casimir invariants are preserved for any time evolution by the equations of the form (8) and (9) since they can be generated by the noncanonical Poisson bracket operator, the co-kernel of which projects onto the symplectic leaves. If we choose these advection fields in the form (10) and (11) with the kernel $(Kij)$ of a definite sign, then the energy changes monotonically and the Casimir invariants are preserved. For the arbitrary choices of $\phi \u0303$ and $J\u0303$, the Casimir invariants are still preserved, although the energy need not change monotonically.

The final step is to perform SA. If the chosen equilibrium is an energy minimum state, which by Dirichlet's theorem is a stable equilibrium (cf. Ref. 2), SA starting from the perturbed state in the second step will recover the original minimum energy state. This is at least true when the perturbation is in a linear regime.

Note that SA can identify a linearly unstable case without a big computational cost since a given perturbation grows in time while the energy is decreased by SA. On the other hand, it is time consuming for a stable case since we have to observe that the given perturbation really disappears by SA. Even if the perturbation becomes smaller by the SA time evolution, the system cannot be said to be stable if the perturbation remains with a small amplitude.

Let us emphasize the importance of the perturbation being dynamically accessible, where the time evolution in the final step occurs on a Casimir leaf determined by the initial condition used for the SA procedure. If the perturbed state were on a different Casimir leaf, the time evolution of SA in the final step would not recover the original equilibrium.

## IV. NUMERICAL RESULTS

### A. A case of stable equilibrium

#### 1. A case study equilibrium for stability analysis

In order to clearly explain the efficacy of our method, in the present paper, we choose an equilibrium with cylindrical symmetry. Any cylindrically symmetric state is an equilibrium of Eqs. (1) and (2). We assume that the equilibrium has no plasma rotation. The safety factor profile is assumed to be $q(r)=q0/(1\u2212r2/2)$ with $q0=1.75$ and shown in Fig. 1. It has *q *=* *2 surface at *r *=* *1/2.

This equilibrium is spectrally stable against ideal MHD modes, although it is unstable to a resonant $(m,n)=(\u22122,1)$ tearing mode. Note that we take the perturbed quantities proportional to $ei\u200a(m\theta +n\zeta )$. Since SA does not change magnetic field line topology, we expect that SA will recover this equilibrium as an energy minimum state.

#### 2. Dynamically accessible perturbation of the equilibrium

The equilibrium in Subsection IV A 1 is perturbed by using arbitrary advection fields, as explained in Sec. III. We choose the advection fields to be

where $A\phi =AJ=10\u22123,\u2009r0=0.5$, and *L *=* *0.1 are all fixed in time. The Fourier-decomposed radial profiles with this procedure are shown in Fig. 2.

The initial condition for the time evolution is the cylindrically symmetric equilibrium. Thus, the $(m,n)=(\u22122,1)$ components of *U* and *ψ* are generated in the beginning phase just after the time evolution is started. Then, the nonlinear terms between the generated $(m,n)=(\u22122,1)$ components and the $(m,n)=(\u22122,1)$ modes of the advection fields generate the $(m,n)=(\u22124,2)$ and (0, 0) components. When the cylindrically symmetric equilibrium is perturbed in this way, the result stays on the same Casimir leaf.

The time evolution of the total energy, the sum of the kinetic energy $Ek$ and the magnetic energy $Em$, and the magnetic helicity $Cm$, which is one of the Casimir invariants, are shown in Fig. 3. When the arbitrary chosen advection fields in Eqs. (14) and (15) are used for perturbing the cylindrically symmetric equilibrium, we observe that the total energy increases in time, while the magnetic helicity remains constant.

#### 3. SA slow convergence of the velocity

Figure 5 shows the time evolution of total energy and the magnetic helicity during the SA evolution. Observe that the total energy monotonically decreases until the system appears to reach a steady state when the energy remains nearly constant for a long time. Note that the horizontal axis is a log scale. Therefore, the state seems to be an equilibrium that has the minimum energy. Also, Fig. 5 shows that the magnetic helicity does not change during the SA evolution.

Figure 6 shows the time evolution of radial profiles of the $(m,n)=(\u22122,1)$ components of (a) $\u2111\u2009U$, (b) $\u2111\u2009\phi $, (c) $\u211c\u2009\psi $, and (d) $\u211c\u2009J$, respectively. Note that $\u211c\u2009U\u22122,1,\u2009\u211c\u2009\phi \u22122,1,\u2009\u2111\u2009\psi \u22122,1$, and $\u2111\u2009J\u22122,1$ remain almost zero during SA. Even though the total energy starts to decrease just after the SA is started, the velocity components do not change visibly in the early stage. The change in the velocity components become visible even after the total energy reaches the stationary state, with these components changing even at $t=5\xd7104$. On the other hand, the magnetic components start to decrease in the early stage of SA. Although $\u211c\u2009J\u22122,1$ still remains finite at $t=3\xd7104$, the magnetic components are negligibly small at $t=5\xd7104$. Thus, the rates of change in the velocity and magnetic components are significantly different.

This may be explained as follows. It is first pointed out that the $(m,n)=(\u22122,1)$ components of *U* and *ψ*, the dominant Fourier modes of the perturbation, are $O(10\u22126)$. Then, $\phi $ is $O(10\u22128)$, while *J* is $O(10\u22123)$. This is because $\phi $ is obtained by integrating *U* especially in the radial direction, while *J* is obtained by differentiating *ψ*. Let us consider the very early stage of SA. The $(m,n)=(\u22122,1)$ component of *f* ^{1}, or the right-hand sides of Eq. (1), is generated mainly by the $[\psi 00,J\u22122,1]$ term. On the other hand, the $(m,n)=(\u22122,1)$ component of *f* ^{2}, or the right-hand sides of Eq. (2), is generated by the $[\psi 00,\phi \u22122,1]$ term. Thus, the $(m,n)=(\u22122,1)$ component of *f* ^{1} is much larger than that of *f* ^{2}.

The artificial advection fields $\phi \u0303$ and $J\u0303$ are given by Eqs. (10) and (11), respectively. Since $(Kij)$ is taken to be diagonal in this study, the $(m,n)=(\u22122,1)$ component of $\phi \u0303$ is much larger than that of $J\u0303$. Then, in the SA Eqs. (8) and (9), the $[\psi 00,\phi \u0303\u22122,1]$ term in Eq. (9) is much larger than the $[\psi 00,J\u0303\u22122,1]$ term in Eq. (8). Therefore, $\psi \u22122,1$ starts to change much faster than $U\u22122,1$. Note that the $[U,\phi \u0303]$ term is negligible especially in the early stage of the SA since *U*_{00} is zero initially.

The SA equations for reduced MHD have two advection fields. Therefore, depending on the relative magnitudes, the relaxation path can change. In the present settings, the relaxation of the magnetic component is much faster than the velocity component. This may be understood schematically by supposing the energy level corresponds to the height of a mountain with topography. Then, the perturbed state corresponds to a point somewhere on a hillside of the mountain and the equilibrium is the place of the lowest height. There can be a variety of paths to the equilibrium. In the present case, the system goes down to the valley quickly and then moves slowly along a valley to the place of the lowest height. Thus, the question is whether the path can be controlled without losing the characteristics of SA: monotonic change in the energy while preserving the Casimir invariants. This would be realized by forcing the two advection fields to have same order of magnitudes, which is what we will propose in Sec. IV A 5.

#### 4. SA relaxation with a forbidden direction

Before devising the advection fields that will speed up SA relaxation, let us first examine what happens if one of the advection fields is fully dropped. Since the system loses one of its fields for the relaxation, we expect the relaxation to be incomplete. Dropping one of the fields can be accomplished by setting either $\alpha 11=0$ or $\alpha 22=0$.

Figure 7 shows the time evolution of (a) total energy, (b) kinetic energy, and (c) magnetic energy, respectively. The $\alpha 11=\alpha 22=100$, the case of Fig. 5, is plotted here for comparison. Note that Fig. 7 is plotted onward from $t=10\u22122$ for showing clearly that the initial energy is the same for all cases.

In the case of $\alpha 11=0$ and $\alpha 22=100,\u2009\phi \u0303$ becomes zero. Then, the magnetic energy does not change in time as shown in Fig. 7(c). This is clear from Eq. (9); *ψ* can only evolve with finite $\phi \u0303$. Although the kinetic energy changes slightly due to the finite $J\u0303$, as seen in Fig. 7(b), the total energy remains almost unchanged since the magnetic energy is dominant.

On the other hand, in the case of $\alpha 11=100$ and $\alpha 22=0,\u2009J\u0303$ becomes zero. Then, the magnetic energy changes in time due to finite $\phi \u0303$. The kinetic energy remains almost unchanged as seen in Fig. 7(b), although it can change slowly due to the nonlinear term $[U,\phi \u0303]$ in Eq. (8). Since the relaxation of the dominant magnetic energy occurs in this case, the time evolution of the total energy almost overlaps the case of $\alpha 11=\alpha 22=100$.

Summarizing, the magnetic energy can never decrease if $\phi \u0303$ is completely dropped; thus, the system does not relax to a minimum energy state. In the case of vanishing $J\u0303$, the system can relax; however, it must be quite slow. Therefore, both advection fields should remain finite and of comparable magnitude for obtaining efficient relaxation. This is the task addressed in Sec. IV A 5.

#### 5. SA accelerated relaxation

Let us now consider a method for maintaining the advection fields to be of the same order of magnitude. In this study, we do this by introducing time dependence in the quantity $(Kij)$ of Eqs. (10) and (11), while retaining positive definiteness. Explicitly, we choose

where $\alpha jj(t)$ are chosen according to

Here, $maxr,m,n|\u2009fmnj(r,t)|$ is the maximum absolute value of the Fourier components of *f ^{j}* in

*r*and $Fmax$ is the desired magnitude of the advection field. The parameter $Fmax$ is taken to be the same for

*j*=

*1 and 2 for balancing the order of magnitude of $\phi \u0303$ and $J\u0303$. Therefore, the first term in the curly braces controls the magnitude of the advection fields, forcing them to have essentially the same order of magnitude. The second term in the curly braces is introduced to avoid divergence of the*

*α*. In the final phase of the relaxation, $maxr,m,n|\u2009fmnj(r,t)|$ becomes vanishingly small and the upper limit on $\alpha max$ avoids the divergence caused by this.

_{jj}Figure 8 shows the time evolution of (a) total energy, (b) *α*_{11} and *α*_{22}, (c) kinetic energy and (d) magnetic energy for $Fmax=10\u22123$ and $10\u22122$. The time evolution of energy without the time dependence of *α _{jj}*, fixed $\alpha 11=\alpha 22=100$, is also plotted in (a), (c), and (d). The upper limit for

*α*was set to $\alpha max=107$. The convergence to the stationary state is clearly seen to be accelerated; for example, when $Fmax=10\u22122$, the required time for reaching the stationary state is considerably faster than the fixed

_{jj}*α*case.

_{jj}Especially, the time relaxation evolution of the kinetic energy starts significantly earlier than the $\alpha 11=\alpha 22=100$ case, as shown in Fig. 8(c), although the order of magnitude of the kinetic energy is small. In the case of fixed *α _{jj}*, the kinetic energy starts to decrease even after reaching the almost stationary state of the total energy. On the other hand, the kinetic energy starts to decrease in the early stage of SA when

*α*is varied in time.

_{jj}As seen in Fig. 8(b), *α*_{22} is significantly larger than *α*_{11}. This is because *f* ^{2} is much smaller than *f* ^{1}. Accordingly, the rapid decrease in the kinetic energy is obtained.

The choice of $Fmax$ and $\alpha max$ can largely affect the convergence efficiency. If $Fmax$ is too large, the simulation can diverge easily without an appropriate step-size control since the right-hand sides of the evolution equations (8) and (9) become too large. A large $\alpha max$ also has similar tendency, although the effect on the divergence is indirect since $\alpha max$ is an upper limit for the right-hand sides of the evolution equation (8) and (9).

Related to the convergence efficiency, there can be a better choice for *α _{jj}* than that of Eq. (18). For example, one can expect defining

*α*for each

_{jj}*m*and

*n*as $\alpha jj,mn$. However, this is not so straightforward since $maxr\u2009|\u2009fmnj(r,t)|$ becomes smaller for larger

*m*and

*n*and thus $\alpha jj,mn$ can become huge or reach the upper limit. This means that the relative magnitude of the right-hand side associated with the evolved quantity

*U*or

_{mn}*ψ*becomes large. Then, the simulation is likely to diverge. In this case, $Fmax$ and $\alpha max$ should be also chosen to have an appropriate order of magnitude for each

_{mn}*U*and

_{mn}*ψ*.

_{mn}Another possible choice for accelerating convergence is to take the symmetric kernel $(Kij)$ in Eqs. (10) and (11) with nonzero off diagonal elements. This mixes *f* ^{1} and *f* ^{2} and can make the advection fields $\phi \u0303$ and $J\u0303$ for SA to be of comparable magnitudes. We have not tried this choice yet. Time dependence of the kernel may be necessary even with this choice to accelerate the convergence near the energy minimum state.

#### 6. SA with another initial perturbation with larger kinetic energy

The SA results shown in the previous subsection IV A 5 were performed for the initial perturbation that has a larger perturbed magnetic energy ($O(10\u22127)$) than the perturbed kinetic energy ($O(10\u221212)$). In order to examine another case where the perturbed kinetic energy is larger than the perturbed magnetic energy, we chose a smaller $\phi \u0303$ in Eq. (14) and a larger $J\u0303$ in Eq. (15) for dynamically accessible perturbations. The total energy of the system increased by the advection. The radial profiles of the initial condition for SA are plotted in Fig. 9. The velocity part has about 200 times larger amplitudes and the magnetic part has about 10 times smaller amplitudes compared with Fig. 4. The perturbed kinetic energy is $O(10\u22127)$, which is larger than the perturbed magnetic energy of $O(10\u22129)$. They are still in a linear regime.

Starting from the initial state shown in Fig. 9, SA was performed. For the acceleration, $Fmax=10\u22123$ and $\alpha max=107$ were used. The total energy decreased by SA as shown in Fig. 10. The simulation was stopped since the Fourier modes of the right-hand sides of Eqs. (8) and (9) for SA as well as the original low-beta reduced MHD Eqs. (1) and (2) became smaller than a threshold $10\u22128$, although the energy looks still to be decreasing. Note that the horizontal axis is a log scale; thus, the decrease in energy is not so rapid.

The perturbation amplitudes decreased in time as shown in Fig. 11. Since the initial magnetic perturbation was small in this case, the amplitudes of the magnetic part quickly became negligibly small. The velocity part is still changing, decreasing the kinetic energy, although it is already small. In fact, not shown in the presented figures, the magnetic energy became almost the same value as the cylindrically symmetric state, while the kinetic energy of $O(10\u22129)$ still remained when the simulation was stopped. However, we observe that the amplitudes of the velocity part were disappearing. The decrease in total energy with the disappearance of perturbation again indicate that the equilibrium is stable.

### B. A case of unstable equilibrium

In this subsection, SA is performed for linear stability analysis of an unstable equilibrium. The safety factor profile is the same as Fig. 1. The poloidal rotation profile is given by

where *α* is a positive parameter. For the following numerical results, $v\theta max=0.01$ and *α* = 3 were used. The profile is plotted in Fig. 12.

This equilibrium was perturbed by the same advection fields of Eqs. (14) and (15) with the amplitudes shown in Fig. 2. The total energy of the system decreased by the dynamically accessible perturbation as shown in Fig. 13. This already means that there exists a state with a lower energy than the cylindrically symmetric equilibrium.

From a perturbed state with *t *=* *10 in Fig. 13, we performed SA. The radial profiles of the initial condition for SA are shown in Fig. 14. The $(m,n)=(\u22122,1)$ components are dominant, although the $(m,n)=(\u22124,2)$ and $(m,n)=(\u22126,3)$ components are still visible. The amplitudes are small and in a linear regime.

The time evolution of the total energy by SA is plotted in Fig. 15. The total energy decreased by SA. For the kernel of the double bracket, $Fmax=10\u22123$ and $\alpha max=107$ were used. After $t\u2243100$, the radial profiles of perturbed quantities started to oscillate in *r*, and the conservation of magnetic helicity was apparently violated. Therefore, the time evolution is plotted only for $t\u2264100$. Linear eigenmode analysis without dissipation for this equilibrium gave us multiple unstable modes, although the numerical technique used was very simple. The eigenmodes look singular. The physical situation seems to be numerically difficult, and oscillation in *r* occurred.

Although the total energy decreased in time by SA, the perturbation amplitudes grew in time as shown in Fig. 16. As noted above, the radial profiles started to oscillate after $t\u2243100$, and thus, another equilibrium without cylindrical symmetry was not obtained. However, the growth of the perturbation under decreasing energy shows that the original equilibrium is at least linearly unstable.

### C. Stability of equilibrium via SA

Let us now consider how SA can be used to find a stable equilibrium. Finding a stable equilibrium is rather time consuming since we have to observe that the given perturbation actually disappears by SA. We tried two kinds of dynamically accessible perturbations: one with a larger perturbed magnetic energy, and the other with a larger perturbed kinetic energy, although both are still in a linear regime. In particular, we consider the equilibrium state obtained in Sec. IV A 5, which has almost zero kinetic energy, it being $O(10\u221215)$. Also the magnetic energy reaches the stationary value and does not seem to further decrease. The magnetic components with finite *m* and/or *n* disappear as time progresses as shown in Fig. 6, although only the $(m,n)=(\u22122,1)$ components are shown. We obtained similar results in Sec. IV A 6 for the other initial perturbation. Therefore, SA seems to recover the original cylindrical symmetric equilibrium. This is of course reasonable since the original equilibrium is linearly stable against ideal MHD modes. In fact, the total energy is almost the same as that of the original equilibrium; it is larger than the original one by $O(10\u221210)$ for the case of Sec. IV A 5 and $O(10\u22129)$ for the case of Sec. IV A 6, of which relative magnitudes to the total energy of the original equilibrium are $O(10\u22129)$ and $O(10\u22128)$, respectively.

Thus, here, we have shown two cases where SA recovers the original equilibrium if it is ideal MHD stable. Another case was performed with $r0=0.8$ in Eqs. (14) and (15), with similar amplitudes as in Fig. 2. The perturbation in this case exists mainly in the region outside the *q *=* *2 surface, which is resonant with the family of $(m,n)=(\u22122,1)$ harmonics. For this case we obtained qualitatively similar results.

For an unstable equilibrium case, SA finds growth of the dynamically accessible perturbation in time under decreasing total energy. We expect the SA procedure, in principle, will converge to another equilibrium if there is one on the same Casimir leaf and the perturbation lies within its basin of attraction. In the case tried in this paper, however, the physical situation seems to be tough for numerical study, and only linear stability was examined.

## V. CONCLUSIONS

In this paper, simulated annealing (SA) was applied to low-beta reduced magnetohydrodynamics (MHD) in cylindrical geometry. Specifically, SA was used to verify the stability of a cylindrically symmetric equilibrium, which is known to be linearly stable or unstable against ideal MHD modes. To this end, a dynamically accessible perturbation was added to the cylindrically symmetric equilibrium by introducing advection fields for perturbing the equilibrium without leaving the Casimir leaf on which the original symmetric equilibrium exists. The dynamically accessible perturbation increased (decreased) the total energy of the system if the original equilibrium is linearly stable (unstable). Then, SA was performed to monotonically decrease the total energy. It was found that the SA reasonably recovers the original, cylindrically symmetric equilibrium when it is stable. On the other hand, the perturbation grew when the original equilibrium is unstable. Since low-beta reduced MHD is a two-field model, it has two advection fields for relaxation. It was found that for practical implementation of SA, it is important to balance the order of magnitude of the two advection fields for SA relaxation to occur on a short timescale. To achieve such efficient relaxation, time dependence was introduced in the symmetric kernel of SA, in such a way as to balance the orders of magnitude of the two advection fields, thereby leading to accelerated convergence to the minimum energy state. The acceleration is particularly necessary for a stable equilibrium, since we have to observe that the perturbation actually disappears.

In essence, the SA method provides a general prescription for building from a Hamiltonian system an alternative non-Hamiltonian system with asymptotic stability to an equilibrium state, with the stability provided by an energy function serving as a Lyapunov function. Convergence to the equilibrium requires that the Lyapunov function have a local extremum, which is enough to show stability for the original Hamiltonian system. Systems that are spectrally stable need not have Lyapunov functions (because of negative energy modes, see e.g., Ref. 2) so the SA method here provides a method that singles out systems that are likely to be nonlinearly stable. Given this general framework, there are many avenues for further explorations in the context of fluid and kinetic plasma models.

## ACKNOWLEDGMENTS

M.F. was supported by JSPS KAKENHI (Grant No. JP21K03507) while P.J.M. was supported by the U.S. Department of Energy (Contract No. DE-FG05–80ET-53088) and the Humboldt Foundation. The authors both would like to acknowledge the JIFT program for the support of MF to visit IFS in the Spring of 2019 when a portion of this work was carried out.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Masaru Furukawa:** Conceptualization (equal); Data curation (lead); Formal analysis (lead); Funding acquisition (equal); Investigation (lead); Methodology (lead); Project administration (equal); Resources (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). **Philip J. Morrison:** Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.