CC BY
4Russian Journal of Nonlinear Dynamics, 2018, vol. 14, no. 2, pp. 217-234. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd180206
MSC 2010: 34C23, 34D20
Analysis of a Mathematical Model for Nuclear Spins
in an Antiferromagnet
This paper is concerned with a system of three nonlinear differential equations, which is a mathematical model for a system of nuclear spins in an antiferromagnet. The model has arisen in recent physical studies and differs from the well-known and well-understood Landau-Lifshitz and Bloch models in the manner of incorporating dissipation effects. It is established that the system under consideration is related to the Landau-Lifshitz system by the passage to the limit only on one invariant sphere. The initial equations contain three dimensionless parameters. Equilibrium points and their stability are examined depending on these parameters. The position of the bifurcation surface is found in the parameter space. It is proved that the corresponding equilibrium is of saddle-node type. Exact statements are illustrated by results of numerical experiments.
Keywords: nonlinear equations, equilibrium, stability, bifurcation
Received December 14, 2017 Accepted January 31, 2018
Leonid A. Kalyakin klenru@mail.ru
Institute of mathematics USC RAS
ul. Chernyshevskogo 112, Ufa, 450008 Russia
L. A. Kalyakin
1. Introduction
Formulation of the problem. The purpose of this paper is to analyze equations for the model of nuclear magnetization of an antiferromagnet in the form that has recently been presented in [5]. These equations for three sought-for functions (x, y, z)(t) with values in R3 can be brought to the form
^ = -z)y-Bzx, ^ = —(A — z)x — Bzy — Az, = B(1 - z2) + Ay. (1.1)
Here, A, B, A ^ 0 are given parameters. They are regarded as nonnegative both in the physical sense and from formal considerations taking into account the symmetries of the equations.
The main goal is to investigate the structure of solutions in the case where the parameters are constant. The primary interest is in equilibrium points, their stability, and bifurcations. The paper also discusses the relation of the system (1.1) to the well-known Landau-Lifshitz model. This work continues the investigation of magnetodynamics models in a rigorous mathematical setting [7-9].
2. The origin of the equations
The parameters of the initial physical model are defined by given magnetic fields. One of them is taken to be constant and homogeneous; the other field, which is orthogonal to the first, rotates fast with frequency w close to the frequency of nuclear magnetic resonance, wn. The initial equations are written in coordinates which rotate with the frequency of the external field, w. This implies representation of the magnetization vector m(T) = (mx,my,mz) in the form mx = mX cos ^(T) + my sin p(T), my = -mX sin p(T) + my cos ^>(T), mz = mZ. The frequency is defined by the derivative of the phase with respect to time, w = ^(T). The equations for the amplitudes (mX ,my ,mZ) are written in the form [5]:
^ = [{un-u)-^mz]my-B0mzmx,
dmY = -[(con -co)- %-mz]mx ~ B0m,zm,y - A0mz, (2.1)
dT ~ j mQ
7<m r, 9 9
= B0(m5 - mZ) + A0my,
dT
where wn, A0, B0, wp > 0 are known coefficients; they are defined by external fields. The constant m0 > 0 represents the equilibrium value of magnetization; the sign of difference of frequencies, wn - w, can be arbitrary. Equations (2.1) yield the relation
+ my + m|) = 2Bomz[m% - (mx + my + mz)],
which implies the existence of the partial integral m\ + my + m?Z = m9. Thus, trajectories of solutions which start on a sphere of radius m0 stay on it for ever, and trajectories starting outside this sphere never evolve onto it. We will call such a sphere invariant.
The number of independent initial parameters can be decreased by renormalizing the variables (by rescaling). In this case, one of the coefficients can always be reduced to unity. When
rescaling time, it is desirable to use a quantity that does not vanish in various physical situations. From this point of view, the difference of frequencies, wn — w, and the amplitude of external excitation, A0, look unsuitable. In general, the scale of time can be chosen in different ways depending on what effects should be detected on times of order unity. We use a scale showing nonlinearity, which determines (by means of the coefficient wp) the dynamical shift of frequency. Therefore, we rescale time as
t
t = j up dT 0
and renormalize the dependent variables
mx(T) = mo • x(t), my(T) = mo • y(t), mx(T) = mo • x(t).
As a result, the equations are brought to the form (1.1), in which three independent dimensionless parameters A = (wn — w)/wp, A = A0/wp, B = B0m0/wp are contained. Now the partial integral takes the form x2 + y2 + z2 = 1.
To analyze this model, it is sometimes more convenient to use cylindrical coordinates r, z, which are related to Cartesian coordinates by x = r cos y = —r sin In these coordinates the equations become
dr
dt
d^
dz 2
^ = z[A sin 0 - Br], r-^--{A-z) = Az cos ip, ^ = B{1 - z2) - At sin tp. (2.2)
dt
In view of the existence of a partial integral, the system on the sphere r2 + z2 = 1 is reduced to two equations:
%=B(1- z2) - Ay/1 - z2 sin ^ = A-z+A—J—cos'0. (2.3)
dt dt 1 z2
In the case where B = 0, the reduced system will be Hamiltonian with the Hamiltonian H = = Az — z2/2 — Ay/l — z2 cos ip. The system describes a dissipation-free oscillator such that almost all its phase trajectories are closed and correspond to periodic solutions. However, when B = 0, there exists a general integral x2 + y2 + z2 = const for the initial equations (1.1). In this case, any sphere turns out to be invariant and reduction to a Hamiltonian system takes place on it.
The positive coefficient B > 0 corresponds to representing relaxation (dissipation) in a specific form proposed in [5]. There are other well-known means for representing relaxation which lead to other forms of the mathematical model: in the form of Bloch or Landau-Lifshitz equations. Such models are well-understood [6, 10] and are therefore not discussed here. The possibility of further simplification of the equations using the smallness of coefficients is not considered; it is assumed that A,B, A w 1.
We elucidate the meaning of the coefficients in relation to the physical interpretion. The controllable parameters of the initial physical system are defined by given (external) magnetic fields. One field is defined by the vector of length H. The other field with amplitude 2H\ is orthogonal to the first and rotates fast with frequency w = ^'(T). Following [5], we show the dependence of the coefficients A, B, A on the initial parameters H, Hi, w of the physical model.
The expression for the coefficient A = (wn — w)/wp contains wn, the frequency of nuclear magnetic resonance (NMR), which is defined by the initial physical system and possibly
by the field H. The sign A depends on the difference wn — w, i.e., on the choice of the pumping frequency w. The field H is involved in determining wp, the coefficient of dynamical frequency shift:
Ci
Uv = ----, Ci, C-2 = const. ^ 0.
p H(H + c2y
In the case of a dissipation-free oscillator, i.e., (2.1) with B0 = 0, the coefficient wp determines the shift (due to nonlinearity) of the eigenfrequency relative to the frequency of the linearized system, wn — w.
Formally the coefficient A = (wn — w)/wp arises when time is normalized. It determines the time scale of the processes described by equations in the form (1.1). Near a resonance, when the frequency of the external field is close to the NMR frequency \wn — w\ ^ wn, the components (x,y,z)(t) represent a slow evolution (on the scale of time t) of the amplitudes of the fast rotating (on the scale T) magnetization vector. In this case, the system (1.1) can be interpreted as the result of averaging (simplification) of more complex equations, which in explicit form are not even discussed. Traces of this simplification are contained in the specific structure of nonlinearities, in which only quadratic terms are retained. Some details of such a simplification are discussed in [5]. Systems like (1.1) are sometimes called equations of the main resonance.
We note the invariance of Eqs. (1.1) under the change of sign (A, z, t) ^ (—A, —z, —t). This allows us to restrict our analysis to the case of a nonnegative coefficient A ^ 0 and to obtain statements for the case A < 0 from symmetry. Of course, in statements about stability we will have to take into account the change of the sign of time t.
The coefficient
H
A = 7 nHn—^-, 7 nHn = const Hw
p
represents, up to a multiplier, the amplitudes of external fields, H\(H + c2). It characterizes the value of pumping with a fast oscillating field. This parameter is regarded as nonnegative, A ^ 0, in view of the invariance of the equations under the transformation x,y, A ^ —x, —y, —A.
The parameter B = B0m0/Hwp is a coefficient related to the rate of nuclear-spin relaxation induced by electron-spin relaxation; an expression for B0 is presented in [5]. This parameter is regarded as nonnegative, B ^ 0, both in its physical sense and from formal considerations due to invariance of the equations under the transformation x, z, A,B, A ^ —x, —z, —A, —B, —A.
Thus, the parameters A and B determine the influence of pumping and relaxation on the nonlinear system (taking into account the dynamical frequency shift). The parameter A characterizes the mismatch of the NMR frequency and the frequency of the external field relative to the dynamical frequency shift.
The amplitudes of the fields H, H\ and the frequency w can depend on time. The dependence on slow time et, 0 < e ^ 1, is a case of practical importance. The small parameter e allows one to use the asymptotics for an approximate analysis of solutions of such weakly nonautonomous systems. It is well known that slow deformations in nonlinear systems sometimes lead to unexpected phenomena whose essence, within the framework of mathematical models, is associated with the loss of stability and rapid rearrangements of solutions. A key role in investigating such problems is played by analysis of a "frozen" system, i.e., at various (constant) parameter values. Such an analysis is made in this paper for Eqs. (1.1).
3. Connection with the Landau — Lifshitz equations
Consider the dynamics of the magnetization vector M(T) for a ferromagnet at given external magnetic fields in the Landau-Lifshitz model. Initial equations for the magnetization vector M are usually written in terms of vector products and involve an effective magnetic field Heff. For a spatially homogeneous system these equations will be ordinary [6, 10]:
^ = -7[M, Heff] - JL[M, [M, Heff}}.
The coefficients y, b = const ^ 0 correspond to a gyromagnetic relation and a damping (relaxation) parameter; they are taken to be constant.
With the general first integral \M(T)| = M0, VM0 = const > 0 in mind, these equations can be rewritten on a sphere of radius M0 for two angles 9, p, which are spherical coordinates of the magnetization vector M = M0(sin 9 cos p, sin 9 sin p, cos 9). The structure of the equations thus obtained depends on the choice of the vector Heff. Below we consider a model incorporating a given external magnetic field defined by the vector h = (hi(cos$, sin$),h0). The transverse and longitudinal amplitudes of the field are given by hi,h0. The axis z is chosen along the direction of the field, so that h0 > 0. The phase $ = $(T) defines the rotation of the transverse component with frequency w = $'(T). In addition, in an effective field one takes into account uniaxial anisotropy with the coefficient K = const > 0, so that Heff = h + 2K/M0(0,0, cos 9). It is the last term that determines the dynamical frequency shift, or, to put it simply, the nonlinearity of the dissipation-free part of the system.
The equations take the simplest form in a coordinate system rotating with the frequency of the external field. Transformation to such a coordinate system corresponds to replacement of the azimuth angle p by the variable ^ = p — $. It is convenient to use, instead of the angle 9, a variable corresponding to one of the Cartesian coordinates, z = cos 9. Thus, the equations are actually considered in cylindrical coordinates. Precisely in these variables, z,^, the Landau-Lifshitz equations take the form [9]:
= 7/74 v/l Sin 0 + [/7.0(1 - z2) - hizVl-z2 cos + ||(1 - z2)-
d^ 2K 1 z b i 1 (31)
W = 7/?'° - w + 7Moz - lhlVT=^cos 0 + M~ohlVT=ism
Under the assumption of smallness of the relaxation coefficient b this system can be simplified. In addition to the trivial passage from b = 0, it is possible to partially take relaxation into account in the form of the equations
^ = 7/Wl-^sin 0 + ^/7,0(1 - z2),
d^ 2K 1 z (32)
W = 7/?'° - w + 7Moz - lhlVT=^cos
Conditions for such a passage are described by the inequalities
¿/7 - §/7,0 - hi/ho « 1.
After rescaling time as t = T ■ 2Kj/M0 and making the change of variable p ^ —p, Eqs. (3.2) simplified in this way transform exactly into the reduced (on the unit sphere) system (2.3) with the coefficients
A = hi M0/2Kj, B = bh0/2KY, A = (7^ — w)M0/2Ky.
Conclusion. In contrast to the Landau - Lifshitz equations, for the system (1.1) there exists a unique sphere on which it reduces to a pair of equations. The reduced system (2.3) is an approximation of the Landau - Lifshitz equations at a small relaxation and a strong longitudinal
4. Stability of the invariant sphere
A consequence of the initial equations in the form (2.2) is the relation
J-(r2 + z2) = 2Bz[l - (r2 + z2)], £>0, (4.1)
which implies the invariance of the unit sphere. It is seen that in the case B > 0 the upper hemisphere r2 + z2 = 1, z > 0, attracts the trajectories and the lower hemisphere r2 + z2 = 1, z < 0, repels them. However, one should not think that, when B > 0, all trajectories are attracted to the upper hemisphere.
Theorem 1. For any constants A, B > 0 and 1 < 5 < 2 there exists a number C > 1 such that in the area of a circular cone:
K = {(z, r, p): 0 < r < (—z)s — C6, z< —C, 0 < p< 2n}
any phase trajectory of the system (1.1) goes to infinity.
The proof involves analysis of the vector field of the system on the boundary of region K. A key role in the proof is played by the velocity components corresponding to the polar radius r and the coordinate z.
On the boundary of the cone r = (—z)6 + C6 we calculate the normal vector n = (1,5(—z)6-1, 0). Obviously, it has zero projection onto the direction of the angle p and is directed "outwards" from the cone. In considering the velocity vector along the trajectory v = (r,z,pp) in polar coordinates (2.2), we calculate its projection onto the normal (v, n):
(v, n) = z[A sin p — Br] + [B (1 — z2) — Ar sin p]5(—z)6-1 =
= zA sin p — Bz((—z)6 — C6) + B5(1 — z2)(—z)6-1 — 5(—z)6-1 [(—z)6 — C6]A sin p.
Rough estimates from above make it possible to get rid of the angle p and to single out the term that is the principal one as z ^to:
(v,n) < (—z)A — B(5 — 1)(—z)6+1 + B5(—z)6-1 + A5(—z)6-1 [(—z)6 + C6] < < —B(5 — 1)(—z)6+1 + (—z)A + B5(—z)6-1 + 2A5(—z)26-1.
Since 1 < 5 < 2, it follows that 5 + 1 > 25 — 1, and in the last expression the term with the highest power (—z)6+1 appears with the negative multiplier —B(5 — 1) < 0. Therefore, at all sufficiently large —z>C = C(A,B,5) the expression on the right is negative and hence (v, n) < 0. This implies that the trajectory that crosses the boundary evolves into the cone.
A separate analysis should be made of the trajectory passing through the vertex of the cone, where r = 0, 2 = —C and a description in polar coordinates does not hold good. A formal proof follows from a comparison of the velocity vector in Cartesian coordinates v = (x, y, z)\r=0z=-c = = (0,AC,B(1 — C2)) with a cone generatrix inclined at the vertex, r'z = —5C6-1. It is seen that, for 1 < 5 < 2 and for a sufficiently large value of C, the inclination of the trajectory will be larger than that of the generatrix: B(C2 — 1)/AC > 5C5-1.
Thus, all trajectories starting from region K do not leave it. Then, in view of Eq. (4.1) and by the theorem on the continuation of a solution [11, p. 13], trajectories must go to infinity. This proves the theorem.
Remark. The rough estimates of the field of directions use only equations for r, z. The parameter A and the equation for — do not appear either in these estimates or in the statement of the theorem. The exact boundary of the domain of attraction to the upper hemisphere depends, of course, on the angle —, and remains indefinite.
Figures 1 and 2 give examples for two solutions with close initial data below the invariant sphere. The trajectories diverge with time: one of them is attracted by the upper hemisphere and the other goes to infinity.
5. Fixed points (equilibrium points)
The fixed points of the system (1.1) are defined by
(A — z)y — Bzx = 0, (A — z)x + Bzy + Az = 0, B (1 — z2) + Ay = 0.
It is easily seen from the resulting relation 2Bz[1 — (x2 + y2 + z2)] = 0 that, when B, A = 0, all fixed points are located on the unit sphere.
We note that the situation is somewhat different for B = 0. In this case, the differential equations (1.1) have the general first integral x2 + y2 + z2 = R2, V R = const. The phase space is foliated by invariant spheres and there exist two to four fixed points on each of them. Such a situation was analyzed in [9].
R
Fig. 1. Numerical solution of Eqs. (1.1). Parameters A = 1, B = 0.19, A = 0.8. Left: trajectory of a solution with initial point (0,0, —1.16) below the invariant sphere. The trajectory is stabilized to a stable equilibrium on the upper hemisphere. Right: stabilization of the spherical radius with time to unity.
R
10
J_i_i_i_I_i_i_i_I_i_i_i_L
2 4 6
8
Fig. 2. Numerical solution of Eqs. (1.1). Parameters A = 1, B = 0.19, A = 0.8. Left: trajectory of a solution with initial point (0, 0, —1.17) below the invariant sphere. The trajectory goes to infinity. Right: an unlimited increase in the spherical radius with time on the trajectory.
Next, analysis is made of equations for fixed points on the unit sphere in the general case for B ^ 0 and A, A > 0. In the cylindrical coordinates z, p this reduces to a pair of equations
It is seen that there are no roots for z2 > 1, and if A = 0, there are no roots at the poles z = ±1. If A = 0, there is no zero root. Thus, Eq. (5.2) should be considered for 0 < z2 < 1.
Lemma 1. If z is the root of the algebraic equation (5.2), then at the equilibrium point the angular coordinate p <E [0,n] is defined uniquely from Eqs. (5.1).
Proof. It follows from the equation for z that on the root
Hence, By/l — z2/A ^ 1. Therefore, the first trigonometric equation is solvable on any root. Moreover, sin p ^ 0; therefore, p G [0,^]. Two trigonometric equations are compatible by virtue of (5.2). Their common solution satisfies the equation tgp = —Bz/A — z. The solution is obviously unique. This proves the lemma.
Let us analyze the position of the fixed points.
Theorem 2. Let A, B, A > 0. Then the system (2.2) has two to four fixed points. All of them lie on the unit sphere r2 + z2 = 1. On the lower hemisphere z < 0 there always exists a unique equilibrium. On the upper hemisphere z > 0 there are one to three fixed points.
Proof. To abbreviate computations, we write Eq. (5.2) in the form
B( 1 - z2) - Ay/1 - z2 sin p = 0, A - z + A z cos p = 0.
y/1 — z2
For the coordinate z we obtain an algebraic fourth-degree equation
P(z; A,B, A) = (z2 — 1)[B2z2 + (A — z)2]+ A2z2 = 0.
(5.1)
(5.2)
z2[A2 — (1 — z2)B2] = (1 — z2)(A — z)2 ^ 0.
P (z; A, A, B) = /3z4 — 2Az3 + [a — 3]z2 + 2Az — A2 = 0
(5.3)
with three independent parameters a = A2 + A2 > 0, 3 = B2 + 1 > 1, A > 0. Since P\z=0 = = —A2 < 0, P\z=1 = A2 > 0, Eq. (5.3) always has a pair of roots, one of which is negative and the other positive.
For the case of four roots we show that there are three roots on the right (when z > 0). In this case, the graph of the fourth-degree polynomial P4(z) = P(z; A, A, B) has three extrema. They are the roots of the equation
P4 = 4/z3 — 6Az2 + 2[a — 3]z + 2A = 0.
The three roots of this equation indicate that the cubic parabola has two extrema. These extrema are the zeros of the second derivative
P" = 12/z2 — 12Az + 2[a — 3].
Since the coefficient A > 0, the largest root of this equation is positive. It corresponds to the minimum of the cubic parabola. In the situation of three roots this minimum is negative. On the other hand, the value in the zero is positive P4(0) = 2A > 0. In this case, two roots of the cubic parabola are on the right when z > 0 and one is on the left when z < 0. These roots correspond to the extrema of the polynomial P4 (z), so that two extrema of this fourth-degree polynomial are on the right. Since this polynomial is negative in the zero, P4(0) = —A < 0, we conclude that three of four zeros of this polynomial must be on the right (when z > 0). This proves the theorem.
Corollary 1. Three fixed points on the upper hemisphere exist if and only if the parameters A, B, A > 0 lie at the intersection D+ n D- of regions each of which is defined by an inequality given in parametric form
D+ = {A,B, A > 0: P'z(z; A , A, B) = 0 , P(z; A , A, B) > 0 , P'z (z; A , A, B) < 0} , D- = {A,B, A > 0: P'z(z; A , A, B) = 0 , P(z; A , A, B) < 0 , P'z (z; A , A, B) > 0}. Here, the variable z £ (0,1) is a parameter.
Proof. The conditions for the existence of four roots for Eqs. (5.3) require that, when z > 0, the polynomial P4(z) = P(z;A,A,B) have one positive maximum and one negative minimum. Each of these conditions defines region D+ or D- in the three-dimensional parameter space (A,B, A) £ R3. The intersection of the regions, D+ n D-, defines the set of parameters for which there exist three equilibria on the upper hemisphere. This proves the corollary.
6. Bifurcation set
One of the important problems in magnetodynamics is to describe collisions (abrupt rearrangements of solutions) which occur under slow deformations of parameters. Such problems have been examined in detail for oscillating systems in the absence of dissipation. The well-known results describe the asymptotics of a solution when the system undergoes a bifurcation of saddle-center type [1-3]. Incorporation of dissipation leads to a bifurcation of saddle-node type. This considerably complicates the analysis of the solution, and there are no results in this direction. Below we analyze a particular system (1.1) to describe the boundary in the parameter space A, B, A on which the number of fixed points changes. This so-called bifurcation set is defined by the presence of a multiple root, i.e., by the pair of equations P(z; A,B, A) = 0, P'z(z; A,B, A) = 0. Here, the variable z £ (0,1) plays the role of a parameter, and it can be eliminated in some way.
6.1. Equation of the bifurcation surface
We first specify the parameter space which can have multiple roots. On such roots the derivative
P'z(z; A, B, A) = B2z(4z2 — 2) + 2(A — z)[—z2 + 1 + z(A — z)] + 2A2z = 0
vanishes. A combination of this relation with the initial equation (5.2) in the form zPZ/2 — P gives the equality
^zP^-P = B2z4-(A-z)(z3-A) = 0. (6.1)
Since a multiple root appears only in the interval z £ (0,1), the last equality can be satisfied only for A ^ 1.
Theorem 3. A bifurcation set occurs only at A ^ 1. It can be represented by the intersection of a pair of surfaces given by
P (zm(A,B, A); A,B, A) = 0, P'z (zm(A,B, A); A,B, A) = 0, (6.2)
where zm(A,B, A) = —L/K is defined by
K = A(a + 5/3 )((a — /)[2/ (a — /3) — 3A2]+A2(a + 5/3)) — — [2/(a — /3) — 3A2] (3A[2/(a — /3) — 3A2] + A3(1 — 4/)), L = A2(1 — 4/3) ((a — /3 )[2/(a — /3) — 3A2] + A2(a + 5/3)) + 2A2[2/(a — /3) — 3A2]2
with a = A2 + A2, /3 = B2 + 1, 0 < A < 1.
Proof. Consider the following equations for the multiple root:
P(z; A, B, A) = /3z4 — 2Az3 + [a — /3]z2 + 2Az — A2 = 0,
(6 3)
i P'z{z; A, B, A) = 2pz3 - 3Az2 + [a - ¡3}z + A = 0.
By successively eliminating terms with the highest power of z these equations can be reduced to one first-degree equation. For example, in the first step after multiplication of the second equation by z and a combination with the first, one obtains the third-degree equation
—Az3 + [a — /]z2 + 3Az — 2A2 = 0.
After that the combination of two third-degree equations leads to the second-degree equation
[2/(a — /3) — 3A2]z2 + A(a + 5/3 )z + A2(1 — 4/) = 0.
After multiplication of the latter by z the combination with one of the third-degree equations leads to another second-degree equation:
((a — / )[2/(a — /) — 3A2 ] + A2(a + 5/3)) z2 + ^3A[2/(a — /) — 3A2] + A3(1 — 4/^ z —
— 2A2[2/(a — /) — 3A2 ] =0.
This yields a first-degree equation in the form K z + L = 0. Formulae for the coefficients L = L(a,3, A), K = K (a,/, A) are presented in the formulation of the theorem. They depend on the method of reduction, but their relationship does not depend on the reduction.
According to the above reasoning, the sought-for multiple root, if any, is defined by the expression zm(A,B, A) = —L/K(a,3, A) with a = A2 + A2, 3 = B2 + 1. It should be emphasized that reduction may result in "redundant" values of zm, which are not the roots of the corresponding initial equation. The formula presented above defines the multiple root zm only for those values of parameters (A,B, A) on which two polynomials P = 0 and P'z = 0 vanish simultaneously. Each of these equations defines the surface in the parameter space (A, B, A). Their intersection (and not each single surface) gives a bifurcation set. This proves the theorem.
Remark. In reducing to the first-order equation it was implicitly assumed that the coefficients of the leading terms do not vanish. It is easy to verify that such exceptional values do not affect the formula for the multiple root and the bifurcation set.
6.2. Bifurcation line at zero relaxation
In order to have an idea of the position of a bifurcation set, it is useful to analyze its limiting position as B — 0.
Theorem 4. The limit of a bifurcation set at the zero coefficient of relaxation B — 0 is a line that consists of two branches
A = 0, A2/3 +A2/3 = 1, VA £ (0,1).
Proof. In the case B = 0 the equations for the multiple root (6.3) take the form
(z2 — 1)(A — z)2 + A2z2 = 0, (A — z)(z3 — A) = 0. (6.4)
For the bifurcation line, this yields two equations in parametric form: A = z, A = 0 and A = z3, A = (1 — z2 )3/2. Eliminating the parameter z, we obtain the required equalities. This proves the theorem.
Theorem 5. When B = 0, three fixed points on the upper hemisphere exist only for parameters A and A from the region between the branches of the bifurcation line: 0 < A < (1 — A2/3)3/2, 0 < A < 1.
Proof. Consider the second derivative of the polynomial
i P'JZ{Z] A, 0, A) = - A) + A2 + A2
on the interval 0 < z < 1. It is easy to see that this expression will be positive for sufficiently large A2 + A2. The monotonicity of this derivative implies that the cubic polynomial P'z(z; A, 0, A) has at most one root for 0 < z < 1 and large A2 + A2. Then the fourth-degree polynomial P(z; A, 0, A) has at most one extremum and hence has only one zero for 0 < z < 1. The simple zeros of the polynomial continuously depend on the parameters A and A, therefore there is only one zero above the bifurcation line A > (1 — A2/3)3/2.
The bifurcation curve can be represented by the equation A2/3 + A2/3 = 1, 0 < A < 1. The properties 0 < A < 1,0 < A < 1 take place in the region below the curve A2/3 + A2/3 < 1. Since the power 2/3 < 1, for points A and A inside this region we have A + A < A2/3 + A2/3 < 1. Consider the value of the polynomial for different values of z £ [0,1]:
P(0; A, 0, A) = —A2 < 0, P(A; A, 0, A) = A2A2 > 0, P(1; A, 0, A) = A2 > 0,
P(z; A, 0, A) = (z2 - 1)(A - z)2 + A2z2\z=a+a = [2(A + A)2 - 1]A2 <0, A + A < 1/^2.
It follows from these inequalities that the polynomial has three zeros on the interval z £ (0,1). By virtue of the continuous dependence of simple zeros on the parameters the number of zeros remains unchanged in the entire region bounded by the bifurcation line. This proves the theorem.
Explanation. If we fix the plane B = const > 0 in the parameter space, the trace of the bifurcation set gives a line separated from the boundary A = 0. Indeed, when A = 0, B > 0, only z = 1 can be a nonnegative root of the algebraic equation. However, in this case the condition of a multiple root does not hold since P'z|z=i = 2B2 + 2(A — 1)2 > 0. Thus, as the relaxation parameter B increases, the bifurcation line moves away from the boundary A = 0.
When using Eqs. (6.2) in practice, it is useful to find the bifurcation line (trace) on the fixed plane B = const. After that one can easily choose curves in the parameter space A, A which cross the region of three roots. Figure 3 gives examples of traces of the bifurcation surface for different values of B. Figure 4 shows the motion of the roots (the coordinate z of the equilibrium points) as the parameter A is deformed. The region of the fold corresponds to the presence of three roots. In the variables A, A the deformation of the root is described by the surface with assembly; Figures 4 correspond to the section of this surface for A = const.
A A
Fig. 3. The bifurcation line on the parameter plane B = const consists of two branches. Their intersection corresponds to a point on the critical line. The value of the relaxation parameter is: on the left B = 0, on the right B = 0.2.
A A
Fig. 4. Change of the coordinates z of the fixed point as the parameter A is deformed. The other parameters are fixed: on the left B = 0.2, A = 0.5, on the right B = 0.2, A = 0.2.
6.3. Critical line
It can be seen from the expression for the second derivative
P'z'z(z; A, B, A) = 12(B2 + 1)z2 — 12Az + 2[A2 +A2 — (B2 + 1)],
taken for the values B = 0, A = 0, A = 1, that Pz'z|z=i = 0. This means that at such parameter values (at the intersection of two branches of the bifurcation line) the corresponding root of the polynomial z = 1 will have multiplicity 3.
In the general case, there is a critical line in the parameter space with B > 0 on the bifurcation surface. At values of the parameters A, B, A from this line the root of the algebraic equation has multiplicity 3. Such a set of values of A, B, A is defined by the equations
P(z; A,B, A) = 0, Pz(z; A,B, A) = 0, P'zz(z; A, B, A) = 0.
Here the variable z £ (0,1) plays the role of a parameter. The equations can be transformed to a simpler form:
Theorem 6. In the region {(A,B, A) £ R3: A> 0, B > 0, 0 < A < 1} on the bifurcation surface there exists a unique critical line which is smooth and is described in parametric form by the relations
A2 = (1 — z2 )3(z2 + 3)/16z2, B2 =3(1 — z2)2(z2 + 3)/16z2, A = z(z2 +3)/4, (6.5) z £ (0,1).
Proof. The equations of the bifurcation set, taken in the form of Eqs. (5.2) and (6.1), are supplemented by the condition for the second derivative:
(B2 + 1)z4 — 2Az3 + [A2 + A2 — (B2 + 1)]z2 — Az — A2 = 0,
B 2z4 — (A — z)(z3 — A) = 0, (6.6)
P'z'z = 12(B2 + 1)z2 — 12Az + 2 [A2 + A2 — (B2 + 1)] = 0.
The last equation in combination with the first in the form P"z2/2 — P = 0 gives
5(B2 + 1)z4 — 4Az3 — 2Az + A2 = 0.
Using the second equation, one can eliminate the term with B2. As a result, one obtains the expression A = z(z2 + 3)/4. After that, from the second equation one writes an expression for B2 in terms of z and then, from the first equation, an expression for A2 in the form (6.5). The smoothness of the line follows from the fact that the tangent vector, calculated by virtue of (6.5), does not vanish and does not go to infinity. This proves the theorem.
Corollary 2. In the region A > 0, B > 0, 0 < A < 1 the bifurcation surface is smooth everywhere except for the critical line, which divides the surface into two smooth leaves.
By the implicit function theorem for the equation Pz(z; A,B, A) = 0 the simple root z = z(A,B, A) with values from (0,1) is a smooth function in (A,B, A) everywhere except for the critical line. Substituting this root yields an equation for the bifurcation surface P(z(A, B, A); A, B, A) = 0. The normal to this surface is defined at all points except for the critical line and is calculated taking into account the vanishing of the derivative with respect to z. For example,
P(z(A,B, A); A,B, A) = cjAP(z; A,B, A)Iz=z(a,b,a) =2Az(A,B, A) = 0.
The bifurcation surface is smooth by virtue of smoothness of the derivatives and and the nonzero normal vector. This proves the corollary.
6.4. Parametric description of the bifurcation surface
The equations of the bifurcation surface in the form (6.2) are convenient for analysis at given parameters A, B, A. However, for a graphic representation of this surface using standard software packages it is more convenient to make a parametric description. It is obtained from the initial equations P(z; A,B, A) = 0, P'z(z; A,B, A) = 0 by solving them for A, B and by using A as an additional parameter. As a result, one obtains a parametric description of the surface in the form
Here the parameters (n,z) G R2 change in the region 0 <z3 < n < z < 1, which corresponds to conditions for the existence of a multiple root.
Figure 5 gives an idea of a bifurcation surface in different projections.
Fig. 5. The bifurcation surface in the parameter space A, B, A consists of two leaves. Their intersection corresponds to the critical line.
7. Stability of simple fixed points
To find fixed points, we solve Eqs. (5.1). They can be considered on an invariant sphere as fixed points for the two-dimensional dynamical system (2.3). The following statement applies to stability on the sphere. We examine the situation for different values of the parameters A > 0, B> 0, A > 0.
Theorem 7. For the two-dimensional system (5.1) on the lower hemisphere —1 ^ z < 0 there always exists a unique fixed point (z-,^-) which is an unstable focus. On the upper hemisphere there always exists a fixed point (z+ which is a stable focus (or node). In the case where on the upper hemisphere there exist three fixed points (zi,^i), (zs), (z2,^2) with coordinates 0 < zi < zs < z2 < 1, two of them (zi(z2,^2) are stable foci (or nodes), and the third point (zs,^s) is a saddle.
Proof. The stability of the equilibrium points is analyzed by the first Lyapunov method. It is convenient to calculate the eigenvalues for the matrix of a system linearized at a fixed point
A = 1 y/r/(z - rj), B = \\/{z -rj){rj - z3), A = r?.
2 — 2 z2 z2
2
by using the Hamiltonian H(z, tp) = Az — z'2/2 — Ay/1 — z'2 cos tp introduced for a dissipation-free system. The matrix of the linearized system has the form
— 2Bz
Hz,z H^,z J
If the equations of the fixed point in the form Hz(z, p) = 0, H^(z, p) = B(1 — z2) are taken into account, the eigenvalues are written in terms of the root z of the algebraic equation (5.2) in the form
A±(z) = -Bz ± \y/(z3 - A)(A - z). (7.1)
Since B, A > 0, it immediately follows from the last formula that if z < 0, then Re \± = = —Bz > 0, Im \± = 0. Therefore, the fixed point on the lower hemisphere will be an unstable focus.
On the upper hemisphere z > 0, additional analysis is required of the roots of the equation P(z; A, B, A) = 0. On these roots, the value of the derivative is considered in the form (6.1)
\ zPl = lzP'z-P = B2Z4 - (A - z)(z3 - A).
Further, we need to take into account the signs of the derivative in simple roots of the polynomial. For roots in the order of increase, 0 < z1 < zs < z2 < 1, the derivative in the roots z1>2 is positive and that in the root zs is negative. If there is only one root, z1 = z+, the derivative in it is positive, and we may assume that zs and z2 are absent in the arguments that follow. This leads to the following inequalities:
(A — zi,2)(z3,2 — A) < B2z4,2, (A — zs)(z? — A) > B2z4s.
The use of these estimates in the expression for the eigenvalues allows the conclusion that for the root zs the numbers A±(zs) will be real and have different signs. Consequently, the corresponding fixed point will be an unstable saddle. In the roots z1>2 the real part of the eigenvalues is negative
ReA±(*i,2) = Bz\ 2 ± -^Re - A)(A - z1>2) < 0.
Therefore, the corresponding fixed points will be stable nodes or foci depending on the sign of the quantity (z32 — A)(A — z1>2). This proves the theorem.
The statement proved above applies to the two-dimensional system on the sphere. This implies that for the initial three-dimensional system (1.1) the analysis concerns stability with respect to perturbations on the sphere. Obviously, the property of instability of two fixed points with coordinates z_ and zs holds also for a complete system. The stability of equilibria with coordinates z1>2 is intuitively clear since these points are on the upper hemisphere, to which close trajectories "press themselves". However, an accurate proof needs to be formalized.
We now consider the initial system in cylindrical coordinates (2.2). In this system we make the following change of variable: the polar radius r is changed to the spherical radius R by the formula R2 = x2 +y2 + z2 = r2 + z2, or, which is the same, r = VR2 — z2. The equations take the form
^ = 2Bz(l- R), i±=B(l- z2) - AVR2 ~ z2 sin 0, at dt
dp (7.2)
—rr=A-z+ * A cos ih.
dt y/R2 _ ~2
The equations for the fixed points
R = 1, B( 1 - z2) - AsjR2 - z2 sin ip = 0, A -z+A— z ms )A = n
VR2 - z2
lead to the known answer: all such points lie on the unit sphere R = 1 and were found in Theorem 2.
We write the matrix of the system (7.2) linearized at the fixed points. Now it is a 3 x 3 matrix and has the form
I -2Bz 0 0 \
Asin'0/Vl — z2 —H^^—'IBz —H^.ip . yzAcos ^/(1 - z2)3/2 Hz,z H^z /
It is seen from the structure of the matrix that another eigenvalue, A3 = -2Bz, appears in addition to the eigenvalues A±, which were calculated above. For the upper hemisphere z > 0 this eigenvalue points to exponential stability in the direction of R, while for the lower hemisphere it points to instability.
Thus, the stable equilibrium points on the upper hemisphere are of focus-node or node-node type. This proves the theorem.
Corollary 3. If the parameter A < 0 is negative, then on the lower hemisphere z < 0 there are one to three equilibrium points; all of them are unstable even with respect to perturbations on the sphere. On the upper hemisphere z > 0 there is always only one equilibrium, which is a stable focus-node.
The proof follows from the previous statement taking into account the invariance of Eqs. (1.1) under the change of sign (A,z,t) ^ (—A, -z, -t).
Remark. In the absence of relaxation, when B = 0, there will be fixed points of center type instead of foci; they were investigated earlier in [9].
7.1. Analysis of a multiple fixed point
In the case where the parameters A, B, A are on the bifurcation surface, there is the root z = zm > 0 of the algebraic equation (5.2), which turns out to be multiple and for which the additional equation P'z = 0 holds. As shown above, the last relation on the multiple root can be reduced to the form B2zm — (A — zm)(zm — A) = 0. This implies that the eigenvalues of (7.1) will be real and one of them will be zero: A+(zm) = —2Bzm < 0, A_ = 0. In this case, the first Lyapunov method is unsuitable for characterization of the stability property. Of course, the negativity of one of the eigenvalues A+(zm) < 0 guarantees the presence of a one-parameter family of solutions which evolve into a multiple fixed point. However, more insight into the structure of the trajectories can be gained only by analyzing the highest terms in the expansion of the two-dimensional system on the sphere near the equilibrium.
The initial system on the invariant sphere is described by the pair of Eqs. (2.3). Considering z, ^ as Cartesian coordinates on the plane, we shift the origin into the fixed point (zm) and make the change of variables £ = z — zm, n = ^ — ^m. The right-hand sides are expanded by the Taylor formula for n — 0 and are written, up to quadratic terms, in the explicit form
4=0^ + bv + a20£2 + an^V + «W + 0((\£\ + |r?|)3),
dt (7.3)
^ = +drj+b20e+&n£r?+w+om + M)3).
The coefficients are calculated in terms of the derivatives of the right-hand sides of (2.3) at the fixed point and obviously depend on the parameters A, B, A.
Theorem 8. At values of the parameters A, B, A > 0 on the bifurcation surface outside the critical line, the multiple fixed point of the two-dimensional system (2.3) is of saddle-node
The proof reduces to verifying the conditions for the point of saddle-node type, as formulated in [4, p. 197]. The matrix of the linear part in (7.3) has a zero determinant A = ad — bc = 0, which corresponds to the presence of a zero eigenvalue (or, which is the same, to multiplicity of the root zm). The sum of diagonal elements is different from zero: a = a + d = —2B < 0. Then, according to [4], the fixed point will be a saddle-node at values of A, B, A for which the following combination of Taylor coefficients for expansion of the right-hand sides is different from zero:
I = 1 (ca-02 - ab02) + , 9 1 7 n9 [ab(abn ~ can) + 26c(a,602 - ca,02)] +
a(a2 + bc) a(a2 + bc)2
H——t^-—r[a262(ca2o - ab-2o) - ab2c(abn - can) - b2c2{ab02 - ca,02)]-
a(a2 + bc)3
This expression is easily brought to the more symmetric form
I = , 0 a—-T[a2(ca02 - a.602) + ab(abn - can) + &2(ca,20 - afoo)]-(a2 + bc)3
The coefficients in the expansion of the right-hand sides of (7.3) are calculated in terms of the derivatives and depend on the parameters (A,B, A) and the chosen point z, p. In view of two equations of the fixed point (5.1) and under the multiplicity condition B2z4 = (A — z)(z3 — A), in the expression for l we can eliminate three variables A,B,p. This yields an expression in terms of z and A:
(z2 - 1)A 16(z-A)(z3- A)
By virtue of the multiplicity condition the denominator does not vanish here for B > 0. Since fixed points occur only for z2 < 1, the numerator vanishes only on the critical line, where A = (3z + z3)/4.
Thus, the relation l = 0, which is considered on the bifurcation surface, is equivalent to equations of the critical line. Formulae for B,A are written in terms of z (i.e., in parametric form) from the multiplicity condition B2z4 — (A — z)(z3 — A) = 0 and from the equation for the fixed point in the form P(z; A, B, A) = 0. This line determines the set of parameters (A, B, A) on the bifurcation set at which alone the equilibrium will not be a saddle-node. This proves the theorem.
Explanation. Near a point of saddle-node type there is a known structure of the trajectories of a two-dimensional system, [4, p. 88]. In the case at hand, there is a stable node sector and two saddle sectors. For the initial three-dimensional system (1.1) this point will be of saddle-node-node type.
8. Conclusion
For the system (1.1), a relation with the Landau-Lifshitz equation has been established. Equilibrium points, their stability and bifurcations have been investigated. The type of equilibrium point at parameters from the bifurcation surface has been established. These results allow
l = 777^—TTrk-rT(3z + z3-4A), z = z„
further investigation of the magnetodynamics model (1.1), for example, in the situation of slowly varying parameters. In addition, the equations written out above for the bifurcation surface can be used for efficient numerical experiments. It is of particular interest to describe the asymp-totics of a solution as it passes slowly the bifurcation point, as was done in the dissipation-free system [1-3, 9]. This problem is a subject for further research.
References
[1] Haberman, R. and Ho, E.K., Boundary of the Basin of Attraction for Weakly Damped Primary Resonance, Trans. ASME J. Appl. Mech, 1995, vol.62, no. 4, pp. 941-946.
[2] Itin, A. P., Neishtadt, A. I., and Vasiliev, A. A., Capture into Resonance in Dynamics of a Charged Partice in Magnetic Field and Electrostatic Wave, Phys. D, 2000, vol. 141, nos.3-4, pp. 281-296.
[3] Kiselev, O.M. and Glebov, S.G., An Asymptotic Solution Slowly Crossing the Separatrix near a Saddle-Centre Bifurcation Point, Nonlinearity, 2003, vol. 16, no. 1, pp. 327-362.
[4] Bautin, N. N. and Leontovich, E. A., Methods and Ways of the Qualitative Analysis of Dynamical Systems in a Plane, 2nd ed., Moscow: Nauka, 1990 (Russian).
[5] Borich, M. A., Bunkov, M.Yu., Kurkin, M.I., and Tankeyev, A. P., Nuclear Magnetic Relaxation Induced by the Relaxation of Electron Spins, JETP Lett., 2017, vol. 105, no. 1, pp. 21-25; see also: Pis'ma v Zh. Eksper. Teoret. Fiz, 2017, vol. 105, no. 1, pp. 23-27.
[6] Gurevich, A. G. and Melkov, G. A., Magnetization Oscillations and Waves, New York: CRC, 1996.
[7] Kalyakin, L. A., Analysis of the Bloch Equations for the Nuclear Magnetization Model, Proc. Steklov Inst. Math., 2013, vol.281, suppl. 1, S64-S81.
[8] Kalyakin, L.A., Sultanov, O.A., and Shamsutdinov, M. A., Asymptotic Analysis of a Model of Nuclear Magnetic Autoresonance, Theoret. and Math. Phys., 2011, vol.167, no. 3, pp. 762-771; see also: Teoret. Mat. Fiz., 2011, vol.167, no. 3, pp. 420-431.
[9] Kalyakin, L.A. and Shamsutdinov, M.A., Adiabatic Approximations for Landau - Lifshitz Equations, Proc. Steklov Inst. Math, 2007, vol. 259, suppl. 2, S124-S140.
[10] Monosov, Ya. A., Nonlinear Ferromagnetic Resonance, Moscow: Nauka, 1971 (Russian).
[11] Nemytskii, V. V. and Stepanov, V. V. Qualitative Theory of Differential Equations, Princeton Math. Ser., vol.22, Princeton, N.J.: Princeton Univ. Press, 1960.
