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.

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 flows5–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 used9,10 for low-beta reduced MHD11 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 MHD13 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 equilibria17 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.

In this paper, the low-beta reduced MHD model11 is investigated. We consider a cylindrical plasma with minor radius a and length 2πR0. The low-beta reduced MHD equations are

Ut=[U,φ]+[ψ,J]εJζ,
(1)
ψt=[ψ,φ]εφζ,
(2)

where all physical quantities are normalized by using the length a, the magnetic field in the axial direction B0, the Alfvén velocity vA:=B0/μ0ρ0 with μ0 and ρ0 being vacuum permeability and typical mass density, respectively, and the Alfvén time τA:=a/vA. The cylindrical coordinates (r,θ,z) as well as the toroidal angle ζ:=z/R0 are used. The inverse aspect ratio is defined as ε:=a/R0. The fluid velocity is v=ẑ×φ, and the magnetic field is B=ẑ+ψ×ẑ, where the unit vector in the z direction is denoted by ẑ. The vorticity is U:=φ and the current density is J:=ψ, where the Laplacian in the rθ plane is denoted by . The Poisson bracket for two functions f and g is defined as [f,g]:=ẑ·f×g.

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

H[u]:=Dd3x12{|(1U)|2+|ψ|2},
(3)

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

Cm:=Dd3xψ(x,t),
(4)
Cv:=Dd3xU(x,t).
(5)

Furthermore,

C1:=Dd3xf(ψh),
(6)
C2:=Dd3xUg(ψh)
(7)

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 ψh(x,t):=ψ(x,t)+ε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.

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:

Ut=[U,φ̃]+[ψ,J̃]εJ̃ζ,
(8)
ψt=[ψ,φ̃]εφ̃ζ,
(9)

where the advection fields are defined as

φ̃(x,t):=Dd3xK1j(x,x)fj(x,t),
(10)
J̃(x,t):=Dd3xK2j(x,x)fj(x,t).
(11)

Here, j =1 or 2, f1 is the right-hand side of Eq. (1), f2 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

Kii(x,x)=αiig(x,x),
(12)

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

g(x,x):=δ3(xx).
(13)

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

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 φ̃ and J̃ 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 φ̃ and J̃, 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.

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/(1r2/2) with q0=1.75 and shown in Fig. 1. It has q =2 surface at r =1/2.

FIG. 1.

Safety factor q profile.

FIG. 1.

Safety factor q profile.

Close modal

This equilibrium is spectrally stable against ideal MHD modes, although it is unstable to a resonant (m,n)=(2,1) tearing mode. Note that we take the perturbed quantities proportional to ei(mθ+nζ). 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

φ̃(r,θ,ζ)=Aφr(1r)e(rr0L)2sin(2θζ),
(14)
J̃(r,θ,ζ)=AJr(1r)e(rr0L)2cos(2θζ),
(15)

where Aφ=AJ=103,r0=0.5, and L =0.1 are all fixed in time. The Fourier-decomposed radial profiles with this procedure are shown in Fig. 2.

FIG. 2.

Radial profiles of advection fields φ̃ and J̃ used for generating dynamically accessible perturbations, i.e., perturbations on the same Casimir leaf as the equilibrium under consideration. (a) φ̃ (real part is zero). (b) J̃ (imaginary part is zero).

FIG. 2.

Radial profiles of advection fields φ̃ and J̃ used for generating dynamically accessible perturbations, i.e., perturbations on the same Casimir leaf as the equilibrium under consideration. (a) φ̃ (real part is zero). (b) J̃ (imaginary part is zero).

Close modal

The initial condition for the time evolution is the cylindrically symmetric equilibrium. Thus, the (m,n)=(2,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)=(2,1) components and the (m,n)=(2,1) modes of the advection fields generate the (m,n)=(4,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.

FIG. 3.

Time evolution of the total energy and the magnetic helicity, where the arbitrary chosen advection fields (14) and (15) are used for perturbing the cylindrically symmetric equilibrium. (a) The total energy increases in time, where Ek and Em are kinetic and magnetic energies, respectively. (b) Magnetic helicity Cm does not change in time.

FIG. 3.

Time evolution of the total energy and the magnetic helicity, where the arbitrary chosen advection fields (14) and (15) are used for perturbing the cylindrically symmetric equilibrium. (a) The total energy increases in time, where Ek and Em are kinetic and magnetic energies, respectively. (b) Magnetic helicity Cm does not change in time.

Close modal

In the reminder of Sec. IV, we perform SA to minimize the total energy of the system by using the advection fields [Eqs. (10) and (11)]. The initial condition is chosen as the state at t =10 in Fig. 3. The radial profiles are shown in Fig. 4.

FIG. 4.

Radial profiles of the perturbed state that is used as an initial condition for SA. The (m,n)=(2,1) components are dominant, but they still have small amplitudes. (a) Radial profile of ℑUmn. The real part is zero. (b) Radial profile of ℑφmn. The real part is zero. (c) Radial profile of ℜψmn. The imaginary part is zero. (d) Radial profile of ℜJmn. The imaginary part is zero.

FIG. 4.

Radial profiles of the perturbed state that is used as an initial condition for SA. The (m,n)=(2,1) components are dominant, but they still have small amplitudes. (a) Radial profile of ℑUmn. The real part is zero. (b) Radial profile of ℑφmn. The real part is zero. (c) Radial profile of ℜψmn. The imaginary part is zero. (d) Radial profile of ℜJmn. The imaginary part is zero.

Close modal

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.

FIG. 5.

Time evolution plots of the total energy and magnetic helicity during the SA evolution. (a) The total energy monotonically decreases in time. (b) Magnetic helicity Cm does not change in time.

FIG. 5.

Time evolution plots of the total energy and magnetic helicity during the SA evolution. (a) The total energy monotonically decreases in time. (b) Magnetic helicity Cm does not change in time.

Close modal

Figure 6 shows the time evolution of radial profiles of the (m,n)=(2,1) components of (a) U, (b) φ, (c) ψ, and (d) J, respectively. Note that U2,1,φ2,1,ψ2,1, and J2,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×104. On the other hand, the magnetic components start to decrease in the early stage of SA. Although J2,1 still remains finite at t=3×104, the magnetic components are negligibly small at t=5×104. Thus, the rates of change in the velocity and magnetic components are significantly different.

FIG. 6.

Radial profiles of the (m,n)=(2,1) components are plotted at several times during the SA evolution. The real parts of U2,1 and φ2,1 and the imaginary parts of ψ2,1 and J2,1 remain zero. The magnetic parts ψ2,1 and J2,1 disappear relatively quickly, while the kinetic parts U2,1 and φ2,1 are changing in the almost stationary state of the energy. (a) Radial profile of ℑU−2,1, (b) radial profile of ℑφ−2,1, (c) radial profile of ℜψ-2,1 and (d) radial profile of ℜJ-2,1.

FIG. 6.

Radial profiles of the (m,n)=(2,1) components are plotted at several times during the SA evolution. The real parts of U2,1 and φ2,1 and the imaginary parts of ψ2,1 and J2,1 remain zero. The magnetic parts ψ2,1 and J2,1 disappear relatively quickly, while the kinetic parts U2,1 and φ2,1 are changing in the almost stationary state of the energy. (a) Radial profile of ℑU−2,1, (b) radial profile of ℑφ−2,1, (c) radial profile of ℜψ-2,1 and (d) radial profile of ℜJ-2,1.

Close modal

This may be explained as follows. It is first pointed out that the (m,n)=(2,1) components of U and ψ, the dominant Fourier modes of the perturbation, are O(106). Then, φ is O(108), while J is O(103). This is because φ 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)=(2,1) component of f1, or the right-hand sides of Eq. (1), is generated mainly by the [ψ00,J2,1] term. On the other hand, the (m,n)=(2,1) component of f2, or the right-hand sides of Eq. (2), is generated by the [ψ00,φ2,1] term. Thus, the (m,n)=(2,1) component of f1 is much larger than that of f2.

The artificial advection fields φ̃ and J̃ are given by Eqs. (10) and (11), respectively. Since (Kij) is taken to be diagonal in this study, the (m,n)=(2,1) component of φ̃ is much larger than that of J̃. Then, in the SA Eqs. (8) and (9), the [ψ00,φ̃2,1] term in Eq. (9) is much larger than the [ψ00,J̃2,1] term in Eq. (8). Therefore, ψ2,1 starts to change much faster than U2,1. Note that the [U,φ̃] term is negligible especially in the early stage of the SA since U00 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 α11=0 or α22=0.

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

FIG. 7.

Time evolution of (a) total energy, (b) kinetic energy, and (c) magnetic energy. If φ̃ is dropped by setting α11=0, the dominant magnetic energy does not change and thus no relaxation occurs. (a) Total energy, (b) kinetic energy, and (c) magnetic energy.

FIG. 7.

Time evolution of (a) total energy, (b) kinetic energy, and (c) magnetic energy. If φ̃ is dropped by setting α11=0, the dominant magnetic energy does not change and thus no relaxation occurs. (a) Total energy, (b) kinetic energy, and (c) magnetic energy.

Close modal

In the case of α11=0 and α22=100,φ̃ 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 φ̃. Although the kinetic energy changes slightly due to the finite J̃, 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 α11=100 and α22=0,J̃ becomes zero. Then, the magnetic energy changes in time due to finite φ̃. The kinetic energy remains almost unchanged as seen in Fig. 7(b), although it can change slowly due to the nonlinear term [U,φ̃] 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 α11=α22=100.

Summarizing, the magnetic energy can never decrease if φ̃ is completely dropped; thus, the system does not relax to a minimum energy state. In the case of vanishing J̃, 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

φ̃(x,t):=α11(t)Dd3xg(x,x)f1(x,t),
(16)
J̃(x,t):=α22(t)Dd3xg(x,x)f2(x,t),
(17)

where αjj(t) are chosen according to

αjj(t)=min{Fmaxmaxr,m,n| fmnj(r,t)|,αmax}.
(18)

Here, maxr,m,n| fmnj(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 φ̃ and J̃. 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 αjj. In the final phase of the relaxation, maxr,m,n| fmnj(r,t)| becomes vanishingly small and the upper limit on αmax avoids the divergence caused by this.

Figure 8 shows the time evolution of (a) total energy, (b) α11 and α22, (c) kinetic energy and (d) magnetic energy for Fmax=103 and 102. The time evolution of energy without the time dependence of αjj, fixed α11=α22=100, is also plotted in (a), (c), and (d). The upper limit for αjj was set to αmax=107. The convergence to the stationary state is clearly seen to be accelerated; for example, when Fmax=102, the required time for reaching the stationary state is considerably faster than the fixed αjj case.

FIG. 8.

Time evolution of (a) total energy, (b) α11 and α22, (c) kinetic energy, and (d) magnetic energy for Fmax=103 or 102. The upper limit was set to αmax=107. The convergence to the stationary state is accelerated. (a) Total energy, (b) α11 and α22, (c) kinetic energy, and (d) magnetic energy.

FIG. 8.

Time evolution of (a) total energy, (b) α11 and α22, (c) kinetic energy, and (d) magnetic energy for Fmax=103 or 102. The upper limit was set to αmax=107. The convergence to the stationary state is accelerated. (a) Total energy, (b) α11 and α22, (c) kinetic energy, and (d) magnetic energy.

Close modal

Especially, the time relaxation evolution of the kinetic energy starts significantly earlier than the α11=α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 αjj is varied in time.

As seen in Fig. 8(b), α22 is significantly larger than α11. This is because f2 is much smaller than f1. Accordingly, the rapid decrease in the kinetic energy is obtained.

The choice of Fmax and α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 αmax also has similar tendency, although the effect on the divergence is indirect since α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 αjj for each m and n as αjj,mn. However, this is not so straightforward since maxr| fmnj(r,t)| becomes smaller for larger m and n and thus α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 Umn or ψmn becomes large. Then, the simulation is likely to diverge. In this case, Fmax and αmax should be also chosen to have an appropriate order of magnitude for each Umn and ψ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 f1 and f2 and can make the advection fields φ̃ and J̃ 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(107)) than the perturbed kinetic energy (O(1012)). In order to examine another case where the perturbed kinetic energy is larger than the perturbed magnetic energy, we chose a smaller φ̃ in Eq. (14) and a larger J̃ 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(107), which is larger than the perturbed magnetic energy of O(109). They are still in a linear regime.

FIG. 9.

Radial profiles of the perturbed state that is used as an initial condition for SA. The (m,n)=(2,1) components are dominant, but they still have small amplitudes. (a) Radial profile of ℑUmn. The real part is zero. (b) Radial profile of ℑφmn. The real part is zero. (c) Radial profile of ℜψmn. The imaginary part is zero. (d) Radial profile of ℜJmn. The imaginary part is zero.

FIG. 9.

Radial profiles of the perturbed state that is used as an initial condition for SA. The (m,n)=(2,1) components are dominant, but they still have small amplitudes. (a) Radial profile of ℑUmn. The real part is zero. (b) Radial profile of ℑφmn. The real part is zero. (c) Radial profile of ℜψmn. The imaginary part is zero. (d) Radial profile of ℜJmn. The imaginary part is zero.

Close modal

Starting from the initial state shown in Fig. 9, SA was performed. For the acceleration, Fmax=103 and α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 108, 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.

FIG. 10.

The total energy decreased in time by SA. The perturbation amplitudes become smaller in time as shown in Fig. 11.

FIG. 10.

The total energy decreased in time by SA. The perturbation amplitudes become smaller in time as shown in Fig. 11.

Close modal

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(109) 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.

FIG. 11.

Radial profiles of the (m,n)=(2,1) components are plotted at several times during SA evolution. The perturbation amplitudes decreased in time. (a) Radial profile of ℑU−2,1, (b) radial profile of ℑφ−2,1, (c) radial profile of ℜψ−2,1, and (d) radial profile of ℜJ−2,1.

FIG. 11.

Radial profiles of the (m,n)=(2,1) components are plotted at several times during SA evolution. The perturbation amplitudes decreased in time. (a) Radial profile of ℑU−2,1, (b) radial profile of ℑφ−2,1, (c) radial profile of ℜψ−2,1, and (d) radial profile of ℜJ−2,1.

Close modal

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

vθ(r)=vθmax(α+1)α+1ααr(1r)α,
(19)

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

FIG. 12.

Poloidal rotation velocity vθ profile.

FIG. 12.

Poloidal rotation velocity vθ profile.

Close modal

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.

FIG. 13.

The total energy decreased in time by the dynamically accessible perturbation.

FIG. 13.

The total energy decreased in time by the dynamically accessible perturbation.

Close modal

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)=(2,1) components are dominant, although the (m,n)=(4,2) and (m,n)=(6,3) components are still visible. The amplitudes are small and in a linear regime.

FIG. 14.

Radial profiles of the perturbed state that is used as an initial condition for SA. The (m,n)=(2,1) components are dominant, but they still have small amplitudes. (a) Radial profile of ℜUmn, (b) radial profile of ℑUmn, (c) radial profile of ℜφmn, (d) radial profile of ℑφmn. (e) Radial profile of ℜψmn. The imaginary part is zero. (f) Radial profile of ℜJmn. The imaginary part is zero.

FIG. 14.

Radial profiles of the perturbed state that is used as an initial condition for SA. The (m,n)=(2,1) components are dominant, but they still have small amplitudes. (a) Radial profile of ℜUmn, (b) radial profile of ℑUmn, (c) radial profile of ℜφmn, (d) radial profile of ℑφmn. (e) Radial profile of ℜψmn. The imaginary part is zero. (f) Radial profile of ℜJmn. The imaginary part is zero.

Close modal

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=103 and αmax=107 were used. After t100, 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 t100. 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.

FIG. 15.

The total energy decreased in time by SA, although the perturbation amplitudes grew in time as shown in Fig. 16.

FIG. 15.

The total energy decreased in time by SA, although the perturbation amplitudes grew in time as shown in Fig. 16.

Close modal

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 t100, 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.

FIG. 16.

Radial profiles of the (m,n)=(2,1) components are plotted at several times during the SA evolution. The perturbation amplitudes grew in time. (a) Radial profile of ℜU−2,1, (b) radial profile of ℑU−2,1, (c) radial profile of ℜφ−2,1, (d) radial profile of ℑφ−2,1, (e) radial profile of ℜψ−2,1, (f) radial profile of ℑψ−2,1, (g) radial profile of ℜJ−2,1, and (h) radial profile of ℑJ−2,1.

FIG. 16.

Radial profiles of the (m,n)=(2,1) components are plotted at several times during the SA evolution. The perturbation amplitudes grew in time. (a) Radial profile of ℜU−2,1, (b) radial profile of ℑU−2,1, (c) radial profile of ℜφ−2,1, (d) radial profile of ℑφ−2,1, (e) radial profile of ℜψ−2,1, (f) radial profile of ℑψ−2,1, (g) radial profile of ℜJ−2,1, and (h) radial profile of ℑJ−2,1.

Close modal

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(1015). 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)=(2,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(1010) for the case of Sec. IV A 5 and O(109) for the case of Sec. IV A 6, of which relative magnitudes to the total energy of the original equilibrium are O(109) and O(108), 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)=(2,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.

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.

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.

The authors have no conflicts to disclose.

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).

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

1.
G. R.
Flierl
and
P. J.
Morrison
, “
Hamiltonian-Dirac simulated annealing: Application to the calculation of vortex states
,”
Physica D
240
,
212
(
2011
).
2.
P. J.
Morrison
, “
Hamiltonian description of the ideal fluid
,”
Rev. Mod. Phys.
70
,
467
(
1998
).
3.
M. D.
Kruskal
and
C. R.
Oberman
, “
On the stability of plasma in static equilibrium
,”
Phys. Fluids
1
,
275
280
(
1958
).
4.
V. I.
Arnol'd
, “
Variational principle for three-dimensional steady-state flows of an ideal fluid
,”
Prikl. Math. Mech.
29
,
846
851
(
1965
)
V. I.
Arnol'd
[
J. Appl. Math. Mech
29
,
1002
1008
(
1965
)].
5.
G. K.
Vallis
,
G. F.
Carnevale
, and
W. R.
Young
, “
Extremal energy properties and construction of stable solutions of the Euler equations
,”
J. Fluid Mech.
207
,
133
152
(
1989
).
6.
G. F.
Carnevale
and
G. K.
Vallis
, “
Pseudo-advective relaxation to stable states of inviscid two-dimensional fluids
,”
J. Fluid Mech.
213
,
549
571
(
1990
).
7.
T. G.
Shepherd
, “
A general method for finding extremal states of Hamiltonian dynamic systems, with applications to perfect fluids
,”
J. Fluid Mech.
213
,
573
(
1990
).
8.
P. J.
Morrison
and
J. M.
Greene
, “
Non-canonical Hamiltonian density formulation of hydrodynamics and ideal magnetohydrodynamics
,”
Phys. Rev. Lett.
45
,
790
(
1980
).
9.
Y.
Chikasue
and
M.
Furukawa
, “
Simulated annealing applied to two-dimensional low-beta reduced magnetohydrodynamics
,”
Phys. Plasmas
22
,
022511
(
2015
).
10.
Y.
Chikasue
and
M.
Furukawa
, “
Adjustment of vorticity fields with specified values of Casimir invariants as initial condition for simulated annealing of an incompressible, ideal neutral fluid and its MHD in two dimensions
,”
J. Fluid Mech.
774
,
443
(
2015
).
11.
H. R.
Strauss
, “
Nonlinear, three-dimensional magnetohydrodynamics of noncircular tokamaks
,”
Phys. Fluids
19
,
134
(
1976
).
12.
M.
Furukawa
and
P. J.
Morrison
, “
Simulated annealing for three-dimensional low-beta reduced MHD equilibria in cylindrical geometry
,”
Plasma Phys. Controlled Fusion
59
,
054001
(
2017
).
13.
H. R.
Strauss
, “
Dynamics of high β tokamaks
,”
Phys. Fluids
20
,
1354
1360
(
1977
).
14.
M.
Furukawa
,
T.
Watanabe
,
P. J.
Morrison
, and
K.
Ichiguchi
, “
Calculation of large-aspect-ratio tokamak and toroidally-averaged stellarator equilibria of high-beta reduced magnetohydrodynamics via simulated annealing
,”
Phys. Plasmas
25
,
082506
(
2018
).
15.
P. J.
Morrison
, “
Structure and structure-preserving algorithms for plasma physics
,”
Phys. Plasmas
24
,
055502
(
2017
).
16.
P. J.
Morrison
, “
A paradigm for joined Hamiltonian and dissipative systems
,”
Physica D
18
,
410
(
1986
).
17.
C.
Bressan
,
M.
Kraus
,
P. J.
Morrison
, and
O.
Maj
, “
Relaxation to magnetohydrodynamics equilibria via collision brackets
,”
J. Phys.: Conf. Ser.
1125
,
012002
(
2018
).
18.
V. D.
Shafranov
, “
On magnetohydrodynamical equilibrium configurations
,”
Sov. Phys. JETP
6
,
545
(
1958
).
19.
R.
Lüst
and
A.
Schlüter
,
Z. Naturforsch.
12A
,
850
(
1957
).
20.
H.
Grad
and
H.
Rubin
, in
Proceedings of the Second United Nations International Conference on the Peaceful Uses of Atomic Energy
(
United Nations
,
Geneva
,
1958
), Vol.
31
.
21.
J. A.
Nelder
and
R.
Mead
, “
A simplex method for function minimization
,”
Comput. J.
7
,
308
313
(
1965
).
22.
J. H.
Holland
, “
Adaptation in natural and artificial systems: An introductory analysis with applications to biology
,”
Control and Artificial Intelligence
(
The University of Michigan Press
,
Ann Arbor
,
1975
).
23.
P. J.
Morrison
and
D.
Pfirsch
, “
Free energy expressions for Vlasov equilibria
,”
Phys. Rev. A
40
,
3898
3910
(
1989
).
24.
P. J.
Morrison
and
D.
Pfirsch
, “
The free energy of Maxwell-Vlasov equilibria
,”
Phys. Fluids B
2
,
1105
1113
(
1990
).
25.
T.
Andreussi
,
P. J.
Morrison
, and
F.
Pegoraro
, “
Hamiltonian magnetohydrodynamics: Lagrangian, Eulerian, and dynamically accessible stability—Examples with translation symmetry
,”
Phys. Plasmas
23
,
102112
(
2016
).
26.
P. J.
Morrison
and
R. D.
Hazeltine
, “
Hamiltonian formulation of reduced magnetohydrodynamics
,”
Phys. Fluids
27
,
886
897
(
1984
).
27.
J. E.
Marsden
and
P. J.
Morrison
, “
Noncanonical Hamiltonian field theory and reduced MHD
,”
Contemp. Math.
28
,
133
150
(
1984
).