CC BY
65
Instantons describing tunneling between magnetic states at finite temperature
S.M. Vlasov1'2, P. F. Bessarab1'2, V.M. Uzdin2'3, H. Jonsson1'4 1Science Institute and Faculty of Physical Sciences, Univ. of Iceland 107 Reykjavik, Iceland 2ITMO University, Kronverkskiy, 49, St. Petersburg, 197101, Russia 3St. Petersburg State University, St. Petersburg, 198504, Russia 4Center for Nonlinear Studies, Los Alamos, NM 87545, USA [email protected], [email protected], [email protected]
DOI 10.17586/2220-8054-2017-8-6-746-759
A method is presented for finding instantons in magnetic systems - optimal paths corresponding to tunneling from one magnetic state to another at a finite temperature. The method involves analytical continuation of the energy to allow for complex values of the angle variables. First, a set of discretization points are placed equally spaced on a chosen energy contour. Then, an estimate of the corresponding temperature is obtained using Landau-Lifshitz dynamics in imaginary time along the contour. Finally, the distribution of the discretization points as well as the energy are systematically refined by converging on the nearest stationary point of the Euclidean action, thereby obtaining a discrete representation of the closest instanton at the given temperature. The method is illustrated with an application to a system consisting of a single spin subject to uniaxial anisotropy and transverse external magnetic field. First-order and second-order crossovers from over-the-barrier mechanism to tunneling are found depending on the applied field, and the difference in the dependence of the instanton temperature on the energy illustrated for the two cases. By comparing the Boltzmann factors for over-the-barrier and tunneling transitions, the crossover temperature between the two mechanisms is estimated for both first- and second-order crossover. Keywords: magnetic transitions, tunneling, instanton, path optimization. Received: 15 November 2017 Revised: 24 November 2017
1. Introduction
Transitions between (meta)stable states of atomic systems are typically thermally activated, i.e. they occur as a result of thermal fluctuations due to coupling to a heat bath. The states are separated and can be identified because of an energy barrier, an energy ridge on the multidimensional energy surface. At high enough temperature, the transitions rely on large energy fluctuations in the right degrees of freedom that enable the system to overcome the energy barrier. This is referred to as over-the-barrier transition mechanism. For high barriers, these are rare events compared with the timescale of vibrations. The rate of such rare events can be estimated from statistical rate theories, in particular, transition state theory and Kramers/Langer theory [1-3]. These theoretical tools make it possible to predict and analyze the typically observed Arrhenius dependence of the transition rate on temperature, where the activation energy is essentially the height of the energy barrier. This methodology was initially developed in the context of atomic rearrangements such as chemical reactions, diffusion events and conformational changes of molecules. More recently, statistical rate theories have been extended to describe transitions from one magnetic state to another [4-8]. More generally, atomic coordinates and orientation of magnetic moments should be treated on an equal footing.
At low enough temperature a different transition mechanism will eventually become dominant. There, the system advances from one state to another by quantum mechanical tunneling through the energy barrier. The crossover from over-the-barrier mechanism to tunneling is evident from a rapid drop in the activation energy as the temperature is lowered. For atomic rearrangements and temperature close to room temperature, this kind of behavior is mostly observed in transitions involving light atoms, such as hydrogen. Transition state theory for atomic rearrangements has been extended by use of Feynman path integrals to take into account the possibility of tunneling [9-13]. Full free energy calculations for determining transition rates have been carried out as well as sampling based on dynamical approaches [14-17].
Most of the time, transition state theory is used within the harmonic approximation where the calculations are simplified by approximating the energy surface with a quadratic expansion in the vicinity of the initial state minimum and at a first-order saddle point on the energy ridge separating the initial state from the final state(s). For over-the-barrier transitions, this leads to the widely used harmonic transition state theory expression for the rate constant [18]. The main challenge in such calculations is finding the relevant first-order saddle point(s) and
evaluating the energy there as well as its second derivatives. If both the initial and final states of the transition are known, the minimum energy path on the energy surface connecting the two local energy minima can be found and thereby the first-order saddle point - the highest energy point on that path. The nudged elastic band (NEB) method is commonly used for this task [19-21]. A more challenging task is to find saddle points on the energy ridge given only the initial state minimum. For this, the minimum mode following method can, for example, be used [22,23].
An analogous harmonic approximation to quantum transition state theory has been developed [10-13,24]. It relies on finding special transition paths corresponding to thermally activated tunneling. Such paths are first-order saddle points on an extended quantum mechanical energy surface underlying the statistical properties of Feynman path integrals. These first-order saddle points are often referred to as 'instantons'. Evaluation of the rate of thermally activated tunneling within instanton theory requires much fewer computations than full quantum transition state theory calculations. Efficient implementations of instanton theory calculations for atomic rearrangements have recently been presented and applied to various chemical reactions [25-27]. Instantons can be found using saddle point search methods, such as the minimum mode following method, or by using a line integral extension of the NEB method.
Much less has been done on quantum mechanical rate theory for magnetic transitions. Most of this work has relied on a mapping of the magnetic system onto a particle system with an effective Hamiltonian, however, this can only be done for certain simple spin systems [28,29]. One such system consists of a single spin with uniaxial anisotropy and a transverse magnetic field [30]. In the particle mapping method, a spin wave function using the Sz eigenstates is constructed and then transformed to an eigenvalue equation i= Ea Schrodinger equation with an effective potential and possibly coordinate dependent mass. The energy spectrum of the spin system coincides with the first 2j+1 levels of the corresponding particle system, where j is the quantum number giving the length of the magnetic momentum vector. Thereby, one can use techniques developed for particle systems to study spin systems. This approach, however, is not universal, since a general way to construct the corresponding particle Hamiltonian is not known. In fact, there are no known strategies for systems with multiple spins, as well as for single spin systems with a Hamiltonian containing higher-order terms (cubic and higher) as is often the case in molecular magnets, for example [31,32]. As far as we know, there is only one study based on spin-particle mapping for a two-spin system [33]. Furthermore, the action of the corresponding particle system does not contain a topological term (14) that affects quantum interference in magnetic systems [34-36]. An efficient method for finding instantons for magnetic systems without relying on spin-particle mapping is, therefore, needed.
Interesting results have, nevertheless, been obtained from studies using the spin-particle mapping technique. It has, in particular, been shown that crossover from over-the-barrier to tunneling mechanism can be of two different types, a smooth second-order crossover where the activation energy changes rather gradually and a more abrupt first-order crossover, where the activation energy drops suddenly as temperature is lowered [30]. A general expression has been presented recently for the second-order crossover temperature applicable to any system for which the energy can be evaluated, including self-consistent field calculations [37,38]. The calculations can be carried out based on the second derivatives of the energy at the first-order saddle point. Such an estimate of the crossover temperature can be useful when assessing the lifetime of magnetic states since it shows whether tunneling needs to be taken into account, while classical rate theory would overestimate the lifetime. An estimate of a first-order crossover temperature is more difficult, however, since it cannot be obtained from the vicinity of the first-order saddle point. For this purpose, it is necessary to find the magnetic instantons for the temperature of interest, but we are not aware of such calculations in the literature.
In this article, we present a method for finding instantons representing tunneling in magnetic systems. The formulation is based on the Euclidean action obtained from the coherent state path integral formulation of the partition function. A numerical technique is presented for finding instantons as stationary points of the action. This approach is applied to a system consisting of a single spin with uniaxial anisotropy in the presence of a transverse magnetic field, and a comparison made with results obtained from spin-particle mapping. The approach can in principle be applied to any magnetic system, including multiple spins, as long as the energy of the system can be described by an analytical function.
2. Theoretical background
A brief review of the theoretical background for coherent spin states and the formulation of path integrals and Euclidean action for spin systems is given below for completeness. The reader is referred to the work by Fradkin for more detailed discussion [39]. An alternative derivation has been given by Kochetov [40].
2.1. Coherent spin states
Coherent states are well known in the context of the quantum harmonic oscillator, where they have been used to relate quantum mechanical dynamics as closely as possible to classical dynamics. They are in many cases a useful addition to other descriptions of a quantum mechanical system, but in the case of spins, they are the only known way to develop a path integral representation. The reason is that the Hilbert space of a spin does not have continuously parameterized eigenstates corresponding to classical variables such as coordinates and momenta. For a system with spin j, (for example, the total magnetic moment of an atom) there are 2j + 1 eigenstates of the Sz operator: Szj) = j), where ^ runs from -j to j. They can be used to construct coherent states, \n), that correspond to particular orientation of the spin specified by continuous spherical polar coordinates, 0 and ^ as:
in = M = E
(2j )!
/ \ i / „, , sin -(j + m)!(j - m) ! v 2j
6
cos -2
j-m
|m, j).
The coherent states satisfy the following closure relation:
2j + 1 ,
1
4n
dn |n)(n|.
Matrix elements of the spin operators in terms of coherent states are:
0
(1) (2)
(3)
(4)
{Sz) = -j + 2j sin2 ^ = -j cos 0,
(Sx) = -j sin 0 cos {Sy) = - j sin 0 sin
2.2. Path integral representation of the partition function
The partition function for a spin system can be written as:
Z = Tr(e-^) = J dn(n|e-^\n). As for Feynman path integrals for particles, the interval p is divided up into N parts:
{n\e-^ii\n) = lim {n\ (1 - eJS ) \n),
N —^^o V /
where e = P/N. Then, the resolution of the identity (2) is inserted between the segments of the interval and the expression for the partition function takes the form:
(5)
(6)
N
Z = Jim / n )(ûfc i (l - eî) |nfc-i),
k=i
(7)
where d^(nk )
2j + 1 4n
dM(n2 - 1) is the invariant integration measure and n0 = nN, i.e. periodic boundary
conditions are fulfilled.
2.3. Euclidean action
The matrix elements in (7) can be approximated up to quadratic terms as:
(nk| (1 - eiï) |nk-i) = (nk|nk-i) 1 - e
(flk |HH |nk-i) (nk |nk-i)
(nk |ii |nk-i) (nk |nk-i)
- (nk |nk-i) exp -e
(8)
The first term in (8) does not equal zero since coherent states are not orthogonal. It leads to an important feature of the spin path integral, a topological phase equivalent to Berry phase [41]. This term can be rewritten as:
(nk|nk-i) = 1 - e (nk|(|nk)-|nk-i)) - exp (-e ^^) - |fik-i))
(9)
m=-j
which in the limit e ^ 0 gives the path integral expression for the partition function as:
Z = J D[n(r)]e-SE, (10)
where the Euclidean action:
¡3
i \ d SE (H) = dr
(n| — |n> + <n|H|n>
(11)
0
is an integral over a closed path Q = (6(r), )) = w(t) on the energy surface E(w). The integration variable, t, can be interpreted as imaginary time.
The first term in the integrand can be rewritten using the expression for the coherent states, (1), as:
H = (cos 2) ^ 2) ' (—j + j (12)
£> = jr^icos i) 2J (tan i)' e-iP^
1 dT> p^V P!(2j -p)! V y V V
|- j + p,j>, (13)
and
p csc 0 — j tan i ) 0 — ip<
d 0 • 0 (n| — |n> = — j tan — 0 + j tan — 0 — ¿j (1 — cos 0)< = —ij (1 — cos 0)<. (14)
Hence, the partition function (10) affords the phase:
exp | f dr(n|7^|n> | = exp | —ijf dr(1 — cos 0)<
= exp | —ijj(1 — cos 0)d< I , (15)
where r is a closed path described by the vector n = [9(t), ^(t)}. This term has a geometric interpretation. It is proportional to the area bounded by the curve r.
The final expression for the Euclidean action of a system with a spin j can be written as:
3/2
S[H] = J dr —ij(1 — cos0)< + E(w) .
(16)
-3/2
2.4. Instantons
Instantons are stationary paths for (16), i.e. they correspond to SS = 0. Differentiation gives the equations of motion:
0 i dE
j sin 0 d<
__i_d_E_ K '
j sin 0 d0 '
where 0 and < satisfy the boundary condition H(0) = H(3). These are the Landau-Lifshitz equations for spin dynamics in imaginary time.
These equations of motion have two types of solutions. The first type corresponds simply to a stationary point on the energy surface, H = w0, with energy E(w0). This can, in particular, be the first-order saddle point on the energy surface wt = The corresponding action is:
S jump = 3E (wt) = j3E t. (18)
This solution characterizes the classical over-the-barrier transition mechanism which dominates at high enough temperature.
The second type of solutions are extended closed paths, so-called instantons. They correspond to tunneling between two potential wells. In the limit of zero temperature T ^ 0, i.e. p ^ to, the instanton corresponds to tunneling from the ground state, but at finite temperature, it corresponds to thermally assisted tunneling. An instanton is a trajectory in imaginary time that conserves energy (analogous to Landau-Lifshitz dynamics in real time) as can be shown by taking the time derivative of E and inserting into (17). This gives:
dE _ dE • dE • i dEdE i dE dE _ ()
dT = "d0 + = j sin o - j sin 0 = 0. (19)
In contrast to particle systems, magnetic systems have no kinetic energy. As a result, in order to find paths that satisfy the boundary conditions and solve (17), it is necessary to analytically continue the function E(w) to allow for complex values of the angles, w.
3. Method for finding magnetic instantons
This section describes the method we have developed for finding instantons for magnetic transitions. Since the instanton corresponds to a stationary point of the action, it corresponds to a periodic orbit described by Landau-Lifshitz equations in imaginary time. As the dynamics conserve energy, the instantons lie on equienergy contours. The first step is, therefore, to identify the contour for a given value of the energy. In order to map out the dependence on temperature, the calculations need to be carried out for a grid of energy values in the interval from the initial state minimum to the first-order saddle point. Once the energy contour has been identified and represented by a set of discretization points, the corresponding temperature is estimated by evaluating the period of the Landau-Lifshitz dynamics along the closed contour. The distribution of the discretization points along the energy contour as well as a refinement of the value of the energy is then carried out using iterative zeroing of the derivatives of the action. The resulting distribution of discretization points then represents an instanton for the given temperature.
3.1. Energy contours
The method for finding and representing energy contours with an even distribution of discretization points is essentially the one previously presented by Einarsdottir et al. [26]. This procedure is a modification of the NEB method [19,20]. The contour is represented with a set of discretization points, k = {1, N}. Each point corresponds to a configuration of the system, i.e. a value of all variables, "k. For a system with a single magnetic moment, "k = (6k, ). More generally, for systems with more than one magnetic moment, "k is a vector with values of 6 and ^ for each one of the magnetic moments. For simplicity, the presentation here will be for a system with a single magnetic moment, but it can be generalized easily. In order to find non-trivial instantons, the angles need to be able to have complex values so a configuration can be written as "k = (6R, 6k, ). The energy of the
system is some known function of the orientation, E("), analytically continued to the complex plane.
In order to find the energy contour corresponding to some particular value of the energy, E', an objective function, SiSO, is defined as the sum of quadratic deviations of the energy:
1N
Sis>i, ...,WN ) = 1 E |E ("k ) - E'|2. (20)
k=0
Given some initial distribution of the discretization points, an iterative algorithm is used to converge on the energy contour corresponding to energy E'. The iterative algorithm is based on zeroing the derivatives of SiSO with respect to the variables.
The negative gradient of the objective function with respect to the variables "k :
gk = -VkSiso = -(E("k) - E')VkE("k), (21)
gives the direction of steepest descent for discretization point k. It represents a force acting on the discretization point that can be used to move it towards the energy contour. As in the NEB method, the displacements of the discretization point parallel and perpendicular to the path are treated separately. A tangent to the path, 4, described by the current location of discretization points is evaluated as in the NEB method [42]. It is used to distinguish between displacements that represent changes in the shape of the path and displacements that are related to the distribution of the discretization points along the path. The component of the gradient representing a change in the shape of the path is given by the force component perpendicular to the path:
g~t = gk - (gk • 4)4. (22)
An even distribution of the discretization points along the path is obtained by including a spring force between adjacent images in the direction parallel to the tangent:
9kP = (kSp(|wfc+1 - uk| - |wfc - wfc_i|)) 4, (23)
where ksp is a spring force chosen in such a way that g^f is roughly of the same magnitude as g^. The total force acting on each discretization point is then:
g°kpt = gi + glp. (24)
A minimization procedure based on the velocity Verlet algorithm is then used to move all the discretization points simultaneously in the direction given by the total force, analogous to what is often done in NEB calculations of minimum energy paths [20].
Figure 1 shows an example calculation where an energy contour is found using the algorithm described above. First, the discretization points are arranged arbitrarily in a circle. Then, they are moved iteratively in the direction of g°pt until an even distribution is obtained on the desired energy contour. This is an arbitrary test case where the energy surface is generated by adding up a few Gaussians for illustration purposes only.
Fig. 1. Illustration of the method for finding and representing an energy contour of a model energy surface. The target contour is shown by the bold red line. Starting with a discrete representation of a circle in configuration space (left), the iterative algorithm moves the discretization points towards the specified energy contour (middle) and eventually converges on it (right).
3.2. Temperature estimate
Given an even distribution of discretization points along an energy contour, the next task is to estimate what temperature an instanton on this energy contour would correspond to. For this purpose, an approximate Landau-Lifshitz dynamics calculation is carried out to estimate the imaginary time it takes to go from one discretization point to another. From the total time it takes to go along the full path, the period, t, can be estimated, and then the temperature as T = 1/t .
Let the closed path described by the set of N discretization points be denoted as Qis° = ms°, ... The time it takes to go from one discretization point to another can be estimated as:
At* = M, (25)
|di|
where the distance between the two points is |d*| = 1^*+° - w*s°| and the velocity can be estimated from Landau-Lifshitz equations as:
( i dE\
j sin 6* d^i -i dE
(26)
\ j sin 6i ddi )
Recall that 0 and $ are complex variables. The total time period corresponding to the closed path along the contour is:
N
EAt,
(27)
i=i
v
This approximate calculation gives an estimate of the instanton temperature associated with the given energy contour. The even distribution of the discretization points along the contour and the temperature estimate are subsequently used as initial conditions for a calculation where an accurate distribution of the discretization points corresponding to an instanton is found, as well as a revised value of the energy for which the obtained temperature value corresponds more closely. This refinement is described in the following section.
3.3. Refinement of discretization points and energy
Given the even distribution of discretization points along an energy contour and the value of the temperature obtained as described above, a refinement calculation is carried out where the temperature is kept fixed but the distribution of the discretization points and the value of the energy are refined so as to converge on a stationary point of the action. This is done by zeroing the derivatives of the action with respect to all variables describing the path, Q = {"i,"2... "N}.
The scaled action, S = PS(Q), is approximated using the set of N discretization points:
s(Q) = " [k(T) ( 1 - cos6k + cos6k+i) (^k+i - ^) + ) + E(6k+i, '^k+i)
k=i
2N
(28)
where k(T) = —ijTkB. The task is to find a stationary point of the function S for a given value of T, that is, to solve the system of differential equations:
dS 1 dE
— = 1 k(T) sin0fc(^fc+i — ¿fc_i) + — =0,
d C 1 d E
— = -k(T)(cos0fc+1 — cos0fc_i) + — =0, k e [1,N]. (29)
d^fc 2 d0fc
Since the action is complex, it is convenient here to define a real objective function, F(Q), by summing up the magnitude squared of the derivatives of S with respect to all the variables. The minimum value of F is zero and corresponds to a stationary point of the action. Standard minimization techniques can be applied to find the minimum of F since it is a real function.
We first define the vector function f (H ) as:
i dS(Q)\
f (Q)
(30)
\ d^N '/
the objective function can be written as the squared norm of f (fi):
F(Q) = 2f(Q) • f(Q).
(31)
The task is to find the N discretization points that give F(Q) = 0 for the given temperature value (obtained as described in the previous section). Since the temperature estimate does not correspond exactly to the energy contour initially chosen, the discretization points will converge on a slightly different energy contour and their distribution will not be uniform.
The BFGS algorithm is used for the iterative minimization of F. It is important to start with the discretization points located on a contour of roughly the right energy to avoid converging on the trivial solution where the discretization points are all at the same location (the classical limit).
The BFGS algorithm requires the first derivatives of F with respect to all the variables. These, in turn, involve the second derivatives of S:
N
dF _ f„ ^
dflR Z^ ^ f)0,.
k = 1
MR
d 2s
00k ddk80j
N
dF _ w dS
doi 9 no,.
dF
dtfR
k = l N
k=l k
N
dF _
d^l ¿I9 d0kJ
+ 9
dS „ d2S
d2S
d0k d0k d0j dS „ d2S
dS „ d2S
d0k d0k 30j d2S
d0k 30k30j dS „ d2S
+ 9- ~ d0k d^j d0k d0k d^j
+ 9
+
O^k d^kdOj dS _ d2S
+ 9
dS ^ d2S
O^k d^kdOj dS „ d2S
, d2S d0k O^j
dS^ d2S
del9d0kd4
+9
O^k d^k d0j
fi-dF
d^k d^kd^j
dS „ d2S
j
+9
d^k d^k d0j' öS? „ d2S
d^k d^k d^j'
d2S
d4k d4k d4j d4k d4k d4j ' Vj e [1,N ].
(32)
4. Model
The methodology described above is applied to a uniaxial spin in transverse magnetic field. The anisotropy axis is taken to be the z-axis and the applied field is pointing along the x-axis. The Hamiltonian of the system is:
H = DS2Z + BSx,
where D is the anisotropy constant and B the strength of the field. The corresponding energy surface is:
B2
E(0, 4) _ Dj2 cos2 0 + Bj sin 0 cos 4 + Dj2 + —.
(33)
(34)
The energy is shifted here by the constant DjZ + BZ/4D to set the minimum energy of the system to zero.
This system has been studied by Chudnovsky using the particle mapping technique [30]. For small field, B < 0.5Dj, the crossover from over-the-barrier mechanism to tunneling is first-order, while it is second-order for larger field. The second-order crossover temperature can be obtained from the second derivatives of the energy at the first-order saddle point on the energy surface [37,38,43]. The expression is [43]:
T
(2)
B(2Dj - B) 2nkB
(35)
The temperature for the first-order crossover can, however, only be found from calculations of the instantons, as illustrated below.
For convenience, the energy function is rather written in terms of dimensionless quantities:
E(0,4) _ - cos2 0 - 2hx sin 0 cos 4 + 1 + h2x,
(36)
where E = DjZE is the scaled energy and hx = B/2Dj is the scaled field strength.
The energy surface is shown in Fig. 2 for two values of the scaled field, hx = 0.50 and 0.05. The two minima corresponding to the ground states of the system as well as the minimum energy path connecting them are located at $ = 0. The larger the field, the further away the minima are from the anisotropy axis. The analytical continuation to complex angles is:
E(d, <£) = - cosZ(dR + id1) - 2hx sin(dR + id1 )cos($R + i4>1) where dR and are the real parts and d1 and are the imaginary parts of the angular variables.
(37)
5. Results
Figure 3 shows calculated instantons for an applied field of hx for all these instantons and the energy is real:
2R
0.50 and 0.05. It turns out that 4R _ 01 _ 0
E_ E _ - cos2(0R) - 2hx sin(0R) cosh(41),
(38)
so the instantons can be visualized in the (dR, ) plane. For higher temperature values, the replicas are nearly equally distributed along the equipotential contours. But, for the lowest temperature, the distribution is distinctly uneven, with larger density of discretization points near the extreme values of dR.
e, rad e, rad
Fig. 2. The energy surface of a system with a single spin subject to anisotropy along the z-axis and an applied magnetic field in x-direction of magnitude hx = B/2Dj = 0.50 (left) and hx = 0.05 (right), see eqn. (34). The local minima corresponding to initial and final states are marked with a white x. For the weaker field, they are located closer to the anisotropy axis. The minimum energy path for the transition is a straight line at ^ = 0 between the minima. The maximum along the minimum energy path, a first-order saddle point on the energy surface, is marked with a red x.
The temperature corresponding to each of the instantons is given in the legend of Fig. 3. The classical first-order saddle point on the energy surface appears as a maximum in this plane. An instanton corresponding to a slightly lower energy contour has all discretization points at nearly the same location, as the path just barely opens up. The corresponding temperature is denoted as T(2) (1.592 K for the higher field, 0.694 K for the lower field).
Instantons on lower energy contours correspond to lower temperatures in the high field case. There is a monotonic decrease in the temperature as the energy is decreased. This behavior can be seen for the high field case in Fig. 4. The crossover from over-the-barrier mechanism to tunneling then occurs at a temperature of T(2) and this is referred to as second-order crossover. The relevant instantons for the onset of tunneling appear at an energy just below E* and correspond to an energy contour lying in the vicinity of the classical first-order saddle point on the real energy surface. There is a continuous transition from over-the-barrier to tunneling mechanism. The value of T(2) can be obtained from the second derivatives of the energy at the saddle point [37,38]. In the limit of zero energy, the temperature approaches zero and the instanton lies on the energy contour that goes through the initial and final state minima.
For the low field case, hx = 0.05, the situation is different. From Fig. 3, it can be seen that the shapes of the instantons are quite different from the high field case. The change in the instanton temperature as the energy is lowered from the saddle point energy is non-monotonic. The temperature first increases as the energy is decreased, reaching a maximum at a value Tm (see Fig. 4) and then decreases, but eventually reaches zero at the energy of the initial state minimum. Two of the instantons shown in Fig. 3 for the low field case have a temperature that is higher than T(2). Even though the instantons are more spread out, and the corresponding periodic orbits longer, the imaginary time period is shorter. Apparently, this is related to the long vertical segments that can be seen in Fig. 3 for the instantons in low field. There, the partial derivative of the energy with respect to 6R is large and the acceleration in given by the Landau-Lifshitz equations has a correspondingly large magnitude.
Figure 5 compares the exponents of the Boltzmann factors for the two mechanisms: over-the-barrier, E* /kB T, and tunneling, Sins/kB T. The smaller this exponent is, the larger the transition rate. For the large field case, hx = 0.50, a smaller exponent is obtained for tunneling as the temperature drops below T(2), consistent with a second-order crossover from over-the-barrier to tunneling mechanism. After a slight increase in the value of Sins/ksT as the temperature is lowered below T(2), it remains relatively independent of temperature.
Fig. 3. Analytically continued energy surface in (0R, ^) plane for field strength of hx = B/2Dj = 0.50 (left) and 0.05 (right). For the system studied here, the instantons are at = 0 and 30 = 0 so it is sufficient to show them in the (3^, K0) plane. Instantons for various energy contours are shown (red, orange, brown and blue lines with dots, upper triangles, and filled circles indicating the location of the discretization points) and the corresponding values of the temperature are given in the legend. The blue instantons are near the first-order saddle point on the real energy surface, corresponding to over-the-barrier mechanism and temperature close to the second-order crossover temperature, T(2). For the stronger field, hx = 0.50, the lower energy and more extended instantons correspond to longer period for the periodic orbit and lower temperature. The crossover from over-the-barrier mechanism to tunneling is second-order. For the weaker field, hx = 0.05, the period decreases and the temperature increases as the energy is lowered from the value at the first-order saddle point until a maximum is reached (notice the high temperature associated with the brown curve with circles). The crossover is first-order and occurs at a temperature T(1) that is higher than T(2).
For the small field case, hx = 0.05, the instantons on energy contours just below the saddle point energy turn out to have a larger Boltzmann exponent than over-the-barrier transitions so they do not lead to crossover to tunneling. The instanton temperature continues to increase as the energy of the contour is lowered until a maximum value is obtained at Tm. Further lowering of the energy then gives a second set of instantons for each temperature value between T(2) and Tm. The Boltzmann exponent for these instantons is smaller and at a certain temperature, T(i), an instanton has the same value of the exponent as the over-the-barrier mechanism, indicating a crossover to tunneling (see Fig. 5). The stationary point of the action that gives the fastest transition rate then changes abruptly from the classical saddle point to an instanton that has significantly lower energy.
Since the first-order crossover corresponds to an instanton that is far removed from the first-order saddle point on the energy surface, it is not possible to estimate T(1) from the energy and its derivatives at the saddle point, as can be done for second-order crossover. Identification of first-order crossover to tunneling requires the evaluation of instantons over the relevant energy range from E = 0 to E^
The difference between first and second-order crossover can be seen from calculations carried out at and slightly below the critical field, hx = B/2Dj = 0.25, as shown in Fig. 6. The second-order crossover is characterized by a monoton decrease of the instanton temperature as the energy decreases. On the other hand, the first-order crossover is characterized by non-monotone behavior. Over a certain range in temperature, instantons on two different energy contours correspond to the same value of the temperature. The instanton lying on the lower energy contour corresponds to a smaller exponent in the Boltzmann factor and therefore higher tunneling rate.
The critical value of hx = 0.25 between first- and second-order crossover had previously been determined for this system by Chudnovsky using the particle mapping technique [30]. The calculations presented here of finite temperature instantons are clearly in good agreement with this previous estimate.
140-
120-
t
; ioo-
ta.
80- l/f(2)
60- E t
0.00 0.05 0.10 0.15
0 0.2 0.4 0.6 0.8 1.0
Fig. 4. Inverse scaled temperature, l/T = 2Dj/TkB, found from the period of Landau-Lifshitz orbits in imaginary time vs. scaled energy, E = Dj2E, of the contour where the instanton lies for field strength of hx = B/2Dj = 0.50 (left) and 0.05 (right). W is the scaled energy of the first-order saddle point on the energy surface corresponding to the activation energy of the classical over-the-barrier mechanism. T(2) is the scaled temperature for an instanton at an energy contour just below . For strong field, hx = 0.50, this is the highest instanton temperature and a second-order crossover to tunneling occurs there. For weak field, hx = 0.05, the crossover is first-order and the instanton involved is not located near the saddle point. There, the instanton temperature increases as the energy is lowered from W and reaches a maximum of Tm for an energy that is intermediate between that of the initial state and the first-order saddle point. The crossover temperature turns out to be higher than T(2) in this case (see Fig. 5).
\ \ 1 \ \ \ \ \
\ \
\
\ \ N
f(2)
0.004 0.006 0.008 0.010 0.012 0.014 0.016 0.018
\ \ \ s N X / \ V t, n. \ / fm
\ \
i i i i 'N
0.0070 0.0075 0.0080 0.0085 0.0090 f
Fig. 5. Exponents of the Boltzmann factor, Sins/kBT for instantons (blue curve) and E^/kBT for over-the-barrier mechanism (dashed orange curve), as a function of scaled temperature, T = TkB/2Dj, for field strength of hx = B/2Dj = 0.50 (left) and 0.05 (right). For strong field, the highest temperature instanton is obtained for energy contour just below the saddle point energy and the crossover temperature, T(2), can be deduced from the second derivatives of the energy at the saddle point [37,38]. For weak field values, the instanton temperature increases as the energy is lowered from the saddle point energy and the Boltzmann exponent for those instantons is larger than for over-the-barrier transitions. After reaching a maximum at Tm, the instanton temperature decreases as the contour energy is lowered and at a temperature of T(i) the exponent for the instanton equals that of the over-the-barrier mechanism indicating a first-order crossover to tunneling.
«n.
<D
1.0
1.5-
1.4-
1.3-
1.2 -
1.1-
0.0
0.2
0.4
0.6
0.8
1.0
P
Fig. 6. Illustration of the change from first-order to second-order crossover to tunneling as the applied field is increased from hx = B/2Dj = 0.20 (orange line (b)) to 0.25 (blue line (a)). Scaled inverse temperature © = Tc/T is shown as a function of scaled energy P = E/Ef. The instanton temperature is a monotone function of the energy when the the crossover to tunneling is second-order. When the crossover is first-order, the temperature corresponding to the instantons first increases as the energy of the contour is lowered below E t, but then decreases for even lower energy.
6. Discussion
A method is presented for finding instantons corresponding to tunneling between magnetic states at finite temperature. The method is based on analytical continuation of the energy of the system to complex angular variables and identification of an energy contour, followed by evaluation of the corresponding temperature and finally refinement of the distribution of discretization points to converge on a stationary point of the action. This method takes advantage of the fact that the Landau-Lifshitz equations of motion for spin systems conserve energy. Since the instanton corresponds to a stationary point of the action, it also corresponds to a periodic orbit along an energy contour. The refinement is based on local minimization of the magnitude of the gradients of the action with respect to the degrees of freedom in the system. Since a good initial guess for this iterative refinement procedure can be obtained from an even distribution of discretization points along the energy contour, such a local minimization most likely will converge on the appropriate instanton, rather than some other stationary point of the action.
This method is quite different from the ones used to find instantons for atomic rearrangements [25,26] where the equations of motion conserve the total energy, i.e. kinetic plus potential energy. In some sense, it is easier to find instantons for magnetic systems since their shape can so clearly be identified on the energy surface. The additional complication, however, lies in the analytical continuation of the energy expression to allow for complex values of the angular variables. While this is relatively straightforward when the energy of the system is given by an analytical expression, it is not clear how to approach this problem when self-consistent calculations (such as the non-collinear Alexander-Andersson model [44]) are used to evaluate the energy.
By finding the instantons as a function of temperature, the crossover temperature for both first-order and second-order crossover can be estimated. Previously, a general expression for the second-order crossover temperature based on second derivatives of the energy at the first-order saddle point on the energy surface had been presented [37,38], however, calculation of the first-order crossover temperature is more challenging because it does not relate directly to the first-order saddle point. There is not a continuous transition from the over-the-barrier to tunneling mechanism in that case. By finding the instantons, however, the first-order case can be treated.
The method has been illustrated here with an application to a system with a single spin. It is, however, easily generalized to systems with an arbitrary number of spins. The method for finding the energy contour, an extension
of the NEB method which is routinely used for multidimensional systems, the Landau-Lifshitz dynamics and the BFGS minimization of the magnitude of the gradient of the action can all be carried out in a straightforward way for many degrees of freedom. We anticipate that tunneling of, for example, magnetic skyrmions and other localized non-collinear states can be studied with the technique presented here.
Acknowledgements
We thank Igor Lobanov for helpful discussions. This work was supported by the Icelandic Research Fund, the Academy of Finland (grant 278260) and the Government of the Russian Federation (grant 074U01).
Appendix
Then, as in the NEB algorithm, using some iterative procedure of minimization, like steepest descent, conjugate gradients method and etc., all the points will be displaced simultaneously at each iteration. The example of searching an isocontour f (x, y) = 4.5 for the function
1 3
f (x, y) = -(kxx2 + ky y2) + £ bi e-a((x-xi)2+(y-yi)2) (39)
2 i=i
using equations (20) and (24) is shown in Fig. 1. The parameters from (39) are shown in Table. 1.
Table 1. Parameters used for the test example.
xi = 1.5 yi =0.9 b1 = 10
x2 = 0.8 y2 = -1 b2 = 5
x3 = -1.2 y3 = 1 b3 = 6
References
[1] Wigner E. The transition state method. Trans. Faraday Soc., 1938, 34, P. 29-41.
[2] Kramers H. A. Brownian motion in a field of force and the diffusion model of chemical reactions. Physica, 1940, 7(4), P. 284-304.
[3] Langer J. S. Statistical theory of the decay of metastable states. Ann. Phys., 1969, 54, P. 258-275.
[4] Braun H.-B. Kramers's rate theory, broken symmetries and magnetization reversal. J. Appl. Physics, 1994, 76, P. 6310-6315.
[5] Coffey W. T., Garanin D. A. and McCarthy D.J. Crossover formulas in the Kramers theory of thermally activated escape rates-application to spin systems. Adv. Chem. Phys., 2001, 117, P. 483-765.
[6] Fiedler G., Fidler J., Lee J., Schrefl T., Stamps R. L., Braun H. B.,Suess D. Direct calculation of the attempt frequency of magnetic structures using the finite element method. J. Appl. Phys., 2012, 111(9), P. 093917(7 pp.)
[7] Bessarab P. F., Uzdin V. M. and Jonsson H. Harmonic Transition State Theory of Thermal Spin Transitions. Phys. Rev. B, 2012, 85, P. 184409(4 pp.)
[8] Bessarab P.F., Uzdin V.M. and Jonsson H. Potential Energy Surfaces and Rates of Spin Transitions. Z. Phys. Chem., 2013, 227, P. 15431557.
[9] Feynman R. P., Hibbs A. R. Quantum Mechanics and Path Integrals. McGraw-Hill, NY, 1965.
[10] Miller W.H. Semiclassical limit of quantum mechanical transition state theory for nonseparable systems. J. Chem. Phys., 1975, 62, P. 1899-1906.
[11] Coleman S. Fate of the false vacuum: Semiclassical theory Phys. Rev. D, 1977, 15, P. 2929-2936.
[12] Callan C., Coleman S. Fate of the false vacuum. II. First quantum corrections. Phys. Rev D, 1977, 16, P. 1762-1768.
[13] Benderskii V.A., Makarov D. E., Wight C. A. Chemical dynamics at low temperatures, John Wiley & Sons, New York, 1994.
[14] Mills G., Schenter G. K., Makarov D. E. and Jonsson H. Generalized path integral based quantum transition state theory. Chem. Phys. Letters, 1997, 278, P. 91-96.
[15] Mills G., Schenter G. K., Makarov D. E., and Jonsson H. RAW quantum transition state theory. in Classical and Quantum Dynamics in Condensed Phase Simulations, edited by Berne B. J., Ciccotti G., and Coker D. F., World Scientific, 1998, P. 405-421.
[16] Craig I.R., Manolopoulos D.E. J. Chem. Phys., 2005, 123, P. 034102.
[17] Habershon S., Manolopoulos D.E., Markland T.E. and Miller III T.F. Annu. Rev. Phys. Chem., 2013, 64, P. 387-413.
[18] Vineyard G. H. Frequency factors and isotope effects in solid state rate processes. J. Phys. Chem. Solids, 1957, 3(1), P. 121-127.
[19] Mills G., Jonsson H., and Schenter G. K. Reversible work based transition state theory: application to H2 dissociative adsorption. Surf. Sci., 1995, 324, P. 305-337.
[20] Jonsson H., Mills G., and Jacobsen K. W., "Nudged elastic band method for finding minimum energy paths of transitions," in Classical and Quantum Dynamics in Condensed Phase Simulations, edited by Berne B. J., Ciccotti G., and Coker D. F., World Scientific, 1998, P. 385-404.
[21] Bessarab P.F., Uzdin V.M. and Jonsson H. Method for finding mechanism and activation energy of magnetic transitions, applied to skyrmion and antivortex annihilation. Comp. Phys. Commun., 2015, 196, P. 335-347.
[22] Henkelman G., Jonsson H. A Dimer Method for Finding Saddle Points on High Dimensional Potential Surfaces Using Only First Derivatives. J. Chem. Phys., 1999, 111, P. 7010-7022.
[23] Gutiérrez M. P., Argáez C. and Jonsson H. Improved Minimum Mode Following Method for Finding First Order Saddle Points J. Chem. Theor. Comput., 2017, 13, P. 125-134.
[24] Richardson J. O. Derivation of instanton rate theory from first principles. J. Chem. Phys., 2016, 144, P. 114106(5 pp.).
[25] Andersson S., Nyman G., Arnaldsson A., Manthe U. and Jonsson H. Comparison of Quantum Dynamics and Quantum Transition State Theory Estimates of the H + CH4 Reaction Rate. J. Phys. Chem. A, 2009, 113, P. 4468-4478.
[26] Einarsdottir D. M., Arnaldsson A., Oskarsson F. and Jonsson H. Path Optimization with Application to Tunneling. Lecture Notes in Computer Science, 2012, 7134, P. 45-55.
[27] Rommel J. B., Liu Y., Werner H.-J. and Kastner J. Role of tunneling in the enzyme glutamate mutase. J. Phys. Chem. B, 2012, 116, P. 13682-13689.
[28] Zaslavskii O. Spin tunneling and the effective potential method. Phys. Lett. A, 1990, 145(8), P. 471-475.
[29] Ulyanov V. V., Zaslavskii O. B. New methods in the theory of quantum spin systems. Phys. Rep., 1992, 216, P. 179-251.
[30] Chudnovsky E. M., Garanin D. A. First- and Second-Order Transitions between Quantum and Classical Regimes for the Escape Rate of a Spin System. Phys. Rev. Lett., 1997, 79, P. 4469-4472.
[31] Wernsdorfer W., Aliaga-Alcalde N., Hendrickson D. N., Christou G. Exchange-biased quantum tunnelling in a supramolecular dimer of single-molecule magnets. Nature, 2002, 416, P. 406-409.
[32] Wernsdorfer W., Murugesu M. and Christou G. Resonant Tunneling in Truly Axial Symmetry Mni2 Single-Molecule Magnets: Sharp Crossover between Thermally Assisted and Pure Quantum Tunneling. Phys. Rev. Lett., 2006, 96, P. 057208(4 pp.)
[33] Owerre S. A. and Paranjape M. B. Phase transition between quantum and classical regimes for the escape rate of dimeric molecular nanomagnets in a staggered magnetic field. Phys. Lett. A, 2014, 378, P. 1407-1412.
[34] von Delft J. and Henley C. L. Destructive quantum interference in spin tunneling problems. Phys. Rev. Lett., 1992, 69, P. 3236-3239.
[35] Loss D., DiVincenzo D. P. and Grinstein G. Suppression of tunneling by interference in half-integer-spin particles. Phys. Rev. Lett., 1992, 69, P. 3232-3236.
[36] Garg A. Topologically Quenched Tunnel Splitting in Spin Systems without Kramers Degeneracy. Europhys. Lett. 1993. 22(3), P. 205-210.
[37] Vlasov S., Bessarab P. F., Uzdin V. M. and Jonsson H. Classical to quantum mechanical tunneling mechanism crossover in thermal transitions between magnetic states. Faraday Discuss, 2016, 195, P. 93-109.
[38] Vlasov S. M., Bessarab P. F., Uzdin V. M. and Jonsson H. Calculations of the onset temperature for tunneling in multispin systems. Nanosystems: Physics, Chemistry, Mathematics, 2017, 8, P. 454-461.
[39] Fradkin E. Field Theories of Condensed Matter Physics. Cambridge University Press, New York, 2013, 190 p.
[40] Kochetov E. SU(2) coherentstate path integral. J. Math. Phys., 1995, 36, P. 4667-4679.
[41] Berry M. V. Quantal Phase Factors Accompanying Adiabatic Changes. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1984, 392(1802), P. 45-57.
[42] Henkelman G., Jonsson H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys., 2000, 113, P. 9978-9985.
[43] Vlasov S. M., Bessarab P. F., Uzdin V. M. and Jonsson H. Crossover temperature for quantum tunnelling in spin systems. J. Phys. Conf. Ser., 2016, 741, P. 012183(5 pp.).
[44] Bessarab P. F., Uzdin V. M. Jonsson H. Calculations of magnetic states and minimum energy paths of transitions using a noncollinear extension of the Alexander-Anderson model and a magnetic force theorem. Phys. Rev. B, 2014, 89(21), P. 214424(12 pp.).
