(In-)Stability of singular equivariant solutions to the Landau-Lifshitz-Gilbert equation

(In-)Stability of singular equivariant solutions to the Landau-Lifshitz-Gilbert equation

Jan Bouwe van den Berg Dept. of Mathematics, VU University Amsterdam, de Boelelaan 1081, 1081 HV Amsterdam, the Netherlands, janbouwe@math.vu.nl    J.F. Williams Dept. of Mathematics, Simon Fraser University, Burnaby, Canada, jfw@math.sfu.ca
July 26, 2019
Abstract

In this paper we use formal asymptotic arguments to understand the stability properties of equivariant solutions to the Landau-Lifshitz-Gilbert model for ferromagnets. We also analyze both the harmonic map heatflow and Schrödinger map flow limit cases. All asymptotic results are verified by detailed numerical experiments, as well as a robust topological argument. The key result of this paper is that blowup solutions to these problems are co-dimension one and hence both unstable and non-generic.

Finite time blowup solutions are thus far only known to arise in the harmonic map heatflow in the special case of radial symmetry. Solutions permitted to deviate from this symmetry remain global for all time but may, for suitable initial data, approach arbitrarily close to blowup. A careful asymptotic analysis of solutions near blowup shows that finite-time blowup corresponds to a saddle fixed point in a low dimensional dynamical system. Radial symmetry precludes motion anywhere but on the stable manifold towards blowup.

The Landau-Lifshitz-Gilbert problem is not invariant under radial symmetry. Nevertheless, a similar scenario emerges in the equivariant setting: blowup is unstable. To be more precise, blowup is co-dimension one both within the equivariant symmetry class and in the unrestricted class of initial data. The value of the parameter in the Landau-Lifshitz-Gilbert equation plays a very subdued role in the analysis of equivariant blowup, leading to identical blowup rates and spatial scales for all parameter values. One notable exception is the angle between solution in inner scale (which bubbles off) and outer scale (which remains), which does depend on parameter values.

Analyzing near-blowup solutions, we find that in the inner scale these solution quickly rotate over an angle . As a consequence, for the blowup solution it is natural to consider a continuation scenario after blowup where one immediately re-attaches a sphere (thus restoring the energy lost in blowup), yet rotated over an angle . This continuation is natural since it leads to continuous dependence on initial data.

Keywords: Landau-Lifshitz-Gilbert, Harmonic map heatflow, Schrödinger map flow, asymptotic analysis, blowup, numerical simulations, adaptive numerical methods

1 Introduction

In this paper we are interested in the existence and stability of finite-time singularities of the Landau-Lifshitz-Gilbert equation for maps from the unit disk (in the plane) to the surface of the unit sphere, :

\hb@xt@.01(1.1)

We will always require the damping term and take without loss of generality. This problem preserves the length of the vector , i.e., for all  implies that for all positive time (for all ).

In the case this equation arises as a model for the exchange interaction between magnetic moments in a magnetic spin system on a square lattice [20, 21]. Taking recovers the harmonic map heatflow which is a model in nematic liquid crystal flow [8]. It is also of much fundamental interest in differential geometry [27]. Finally, the conservative case is the Schrödinger map from the disk to the sphere, which is a model of current study in geometry [12, 16, 17].

Stationary solutions in all these cases are harmonic maps. This allows us to analyze singularity formation in a unified manner for all parameter values. Traditionally, the Landau-Lifshitz-Gilbert equation is posed with Neumann boundary conditions, but this does not affect the local structure of singularities, should they arise. We note that in the harmonic map literature the second term on the right-hand side of the differential equation (LABEL:LLG1) is often rewritten using the identities

where denotes the inner product in and , with the components of .

As discussed in much greater detail in Sections LABEL:sec:formulation and  LABEL:sec:previousresults, there are initial data for which the solution to (LABEL:LLG1) becomes singular in finite time. In this paper we analyze this blowup behaviour, in particular its stability properties under (small) perturbations of the initial data. Considering initial data that lead to blowup, the question is whether or not solutions starting from slightly different initial data also blowup. Our main conclusion is that blowup is an unstable co-dimension one scenario. With this in mind we also investigate the behavior of solutions in “near-miss” of blowup, and the consequences this has for the continuation of the blowup solution after its blowup time.

1.1 Problem formulation

We will consider two formulations for equation (LABEL:LLG1). The first is the so-called equivariant case: using polar coordinates on the unit disk , these are solutions of the form

\hb@xt@.01(1.2)

which have the (intertwining) symmetry property for all  and each fixed , where is a rotation over angle around the origin in the plane , while is a rotation over angle around the -axis in .

The components then satisfy the pointwise constraint , as well the differential equations

\hb@xt@.01(1.3)

where

We will take in what follows, except in Section LABEL:sec:AsympN.

Alternatively, we can parametrize the solutions on the sphere via the Euler angles:

\hb@xt@.01(1.4)

where the equations for and are given by

\hb@xt@.01(1.5)

We note that due to the splitting in (LABEL:LLG_2form) the system (LABEL:LLG_2comp) has one spatial variable. In this equivariant case the image of one radius in the disk thus fixes the entire map (through rotation) and we write .

In the special case and only, there are radially symmetric solutions of the form  constant, reducing the system to a single equation

\hb@xt@.01(1.6)

1.2 Previous results

It is well known that not all strong solutions to the radially symmetric harmonic map heatflow (LABEL:radHM) are global in time. Equation (LABEL:radHM) is -periodic in . Supplemented with the (finite energy) boundary condition , it only has stationary solutions of the form for . Hence with prescribed boundary data and there is no accessible stationary profile. However, there is an associated Lyapunov functional,

whose only stationary points are the family . It is this paradox that leads to blowup: there is a finite collection of (possibly finite) times at which “jumps” from to , losing of energy [27, 6]. The structure of the local solution close to the jumps (in time and space) is known, which allows us to analyze the stability of these solutions.

The fundamental result in this area is due to Struwe [27] who first showed that solutions of the harmonic map heatflow could exhibit the type of jumps described above and derived what the local structure of the blowup profile is. Chen, Ding and Ye [11] then used super- and sub-solution arguments, applicable only to (LABEL:radHM), to show that finite-time blowup must occur when and . The blowup rate and additional structural details were determined through a careful matched asymptotic analysis in [29].

The analysis is based on the original result of Struwe who showed that any solution which blows up in finite time must look locally (near the blowup point) like a rescaled harmonic map at the so-called quasi-stationary scale. That is, there is a scale on which the solution takes the form

\hb@xt@.01(1.7)

From [29] it is known that for generic initial data one has

\hb@xt@.01(1.8)

for some and blowup time , which both depend on the initial data. This result is intriguing as the blowup rate is very far from the similarity rate of  [3].

While the blowup rate (LABEL:RadRate) was derived in [29] for the harmonic map heatflow, i.e. , , in this paper we demonstrate that formal asymptotics imply that this rate is universal for all parameter values of and . Very recently, proofs of the blowup rate (LABEL:RadRate) have appeared for the harmonic map heatflow [24] as well as the Schrödinger map flow [23], i.e., the two limit cases of the Landau-Lifshitz-Gilbert problem.

It is common to consider radial symmetry when analyzing the blowup dynamics of many reaction-diffusion equations. Typically there one can show that there must be blowup using radially symmetric arguments. Moreover, numerical experiments generically show that rescaled solutions approach radial symmetry as the blowup time is approached.

For the harmonic map problem the proof of blowup solutions due to Cheng, Ding and Ye is completely dependent on the radial symmetry. Moreover, there are stationary solutions to the problem which are in the homotopy class of the initial data, but which are not reachable under the radial symmetry constraint. This begs the question: What happens when we relax the constraint of radially symmetric initial data?

The above description is mainly restricted to the harmonic map problem () which has received considerably more attention than the general Landau-Lifshitz-Gilbert equation. Before addressing the question of stability under non-radially symmetric constraints for the harmonic map heatflow problem, we first show that blowup solutions are still expected for the full Landau-Lifshitz-Gilbert equation with , see Section LABEL:sec:Top. We note that

is a Lyapunov functional for the equivariant problem (LABEL:LLG_2comp) as long as , whereas it is a conserved quantity for the Schrödinger map flow ().

We discuss the question of stability for the full problem in a uniform manner. The topological argument in Section LABEL:sec:Top suggests that blowup is co-dimension one, and this is indeed supported by the asymptotic analysis in Section LABEL:sec:Asymp and the numerics in Section LABEL:sec:Num. The main quantitative and qualitative properties turn out to be independent of the parameter values, except for the angle between sphere that bubbles off and the remaining part of the solution. In Section LABEL:sec:MoreAsymp we analyze near-blowup solutions. These solution rotate quickly over an angle in the inner scale. For the blowup solution this implies a natural continuation scenario (leading to continuous dependence on initial data) after the time of blowup: the lost energy is restored immediately by re-attaching a sphere, rotated over an angle . Finally, in Section LABEL:sec:AsympN we present the generalization to the case , followed by a succinct conclusion in Section LABEL:sec:Concl.

2 The global topological picture

We present a topological argument to corroborate that blowup is co-dimension one. It does not distinguish between finite and infinite time blowup. Since the argument relies on dissipation, it works for , but since the algebra is essentially uniform in and , as we shall see in Section LABEL:sec:Asymp, we would argue that the situation for the Schrödinger map flow is the same.

Let us first consider the equivariant case, where, as explained in Section LABEL:sec:formulation, the image of one radius in the disk fixes the entire map, and we write .

Let (the north pole) and . In the notation using Euler angles from Section LABEL:sec:formulation, by rotational symmetry we may assume that , . The only equilibrium configuration satisfying these boundary conditions is , where . Note that for there is no equilibrium, hence blowup must occur for all initial data in that case [2].

For , i.e. , we shall construct a one parameter family of initial data , and we argue that at least one of the corresponding solutions blows up. Since the presented argument is topological, it is robust under perturbations, hence it proves that blowup is (at most) co-dimension one. The matched asymptotic analysis in Section LABEL:sec:Asymp confirms this co-dimension one nature of blowup.

Figure 2.1: One parameter families of initial data for the equivariant case; several members of half of each family are shown (the other half lives on the hemisphere facing away from us). Left: for one may obtain such a family for example by stereographic projection (w.r.t. ) of all straight lines through the point in the plane corresponding to . Right: for one can choose the stereographic projection of parallel lines covering the plane.

We choose one-parameter families of initial data as follows. The family of initial data will be parametrized by , or with the end points identified. Let be a continuous map from to , such that

\hb@xt@.01(2.1)

We may then view as a map from by identifying and to points. Now choose any continuous family satisfying (LABEL:e:endpoints) such that it represent a degree 1 map from to itself. One such choice is obtained by using the stereographic projection

Let be such that , i.e. . Then we choose

see also Figure LABEL:f:cover1. For the special case that we choose

\hb@xt@.01(2.2)

This is just one explicit choice; any homotopy of this family of initial data that obeys the boundary conditions (LABEL:e:endpoints) also represents a degree 1 map on . Let be the collection of initial data obtained by taking all such homotopies. It is not hard to see that is the space of continuous functions (with the usual supremum norm) satisfying the boundary conditions. Let be the subset of initial data in for which the solution to the equivariant equation (LABEL:LLG_3comp) blows up in finite or infinite time. The following result states that the co-dimension of is at most one. In particular, each one parameter family of initial data that represent a degree 1 map from to itself has at least one member that blows up.

Proposition 2.1

Let . The blowup set for the equivariant equation (LABEL:LLG_3comp) has co-dimension at most 1.

Proof. Let be any family of initial data that, via the above identification, represent a degree 1 map from to itself. Let correspond to the solution with initial data . As explained above, we see from (LABEL:e:endpoints) that we may view as a map from by identifying and to points. Since the boundary points are fixed in time, we may by the same argument view as a map from to itself along the entire evolution. Note that since the energy is a Lyapunov functional for , any solution tends to an equilibrium as . If none of the solutions in the family would blow up (in finite or infinite time) along the evolution, then all solutions converge smoothly to the unique equilibrium. In particular, for large the map has its image in a small neighborhood of this equilibrium, hence it is contractible and thus has degree 0. Moreover, if there is no blowup then is continuous in , i.e. a homotopy. This is clearly contradictory, and we conclude blowup must occur for at least one solution in the one-parameter family . This is a topologically robust property in the sense that any small perturbation of also represents a degree 1 map, and the preceding arguments thus apply to such small perturbations of as well. This proves that the co-dimension of is at most 1.     

One may wonder what happens when (equivariant) symmetry is lost. Although a priori the co-dimension could be higher in that case, we will show that this is not so. For convenience, we only deal with boundary conditions for all , which simplifies the geometric picture, but the argument can be extended to more general boundary conditions.

The only equilibrium solution in this situation is  [22]. Let be the family of initial data for the equivariant case with boundary condition (see (LABEL:e:mN) and Figure LABEL:f:cover1). Consider now the following family of initial data for the general case:

We see that maps to , and , but also . We may thus identify to a point, and interpret as a map from to . In particular, represents an element in the homotopy group . Furthermore, upon inspection, represents the generator of the group, since it is (a deformation of) the Hopf map (see e.g. [18]). Let be the collection of initial data in one parameter families obtained from all homotopies of that obey the boundary conditions

\hb@xt@.01(2.3)

Let be the subset of initial data in for which the solution to the differential equation (LABEL:LLG1) blows up in finite or infinite time. As before, the co-dimension of is at most 1, showing that dropping the equivariant symmetry does not further increase the instability of the blowup scenario.

Proposition 2.2

Let and . The blowup set for the general equation (LABEL:LLG1) has co-dimension at most 1.

Proof. The proof is analogous to the one of Proposition LABEL:prop:one, but one uses instead of , i.e. the degree, to obtain the contradiction.     

As a final remark, even though blowup is co-dimension one, this does not mean it is irrelevant. Clearly, by changing the initial data slightly one may avoid blowup. On the other hand, the arguments above indicate that blowup is caused by the topology of the target manifold, and one can therefore not circumvent this type of singularity formation by simply adding additional terms to the equation (for example a physical effect that works on a smaller length scale), unless additional equilibria are introduced which reflect the pinning of a defect.

3 Asymptotic analysis

In this section we extend the results of [29], where the rate of blowup for radially symmetric solutions to the harmonic heat map problem (LABEL:radHM) was determined. We will consider both the extension to the full Landau-Lifshitz-Gilbert equation (i.e. ), as well as allowing a particular class of non-radial perturbations. We find that blowup solutions are always unstable in the equivariant regime. It can be understood that the blowup solutions are separatrices between two distinct global behaviors.

3.1 The inner region

We will proceed with an expansion motivated by two facts: (i) blowup in the harmonic map heatflow is a quasi-static modulated stationary solution; (ii) the full LLG problem has the same stationary profiles as the harmonic map heatflow.

Without specifying the rescaling factor yet we introduce the rescaled variable

which transforms (LABEL:LLG_2comp) to

\hb@xt@.01(3.1)

To solve this in the limit we pose the expansion

where

\hb@xt@.01(3.2)
\hb@xt@.01(3.3)

represent slow movement along the two-parameter family of equilibria and .

At the next order we have

\hb@xt@.01(3.4)
\hb@xt@.01(3.5)

These equations can both be solved exactly, but we omit the algebraic details since for the matching we only need the asymptotic behaviour as , viz.:

\hb@xt@.01(3.6)
\hb@xt@.01(3.7)
Figure 3.1: Regions for asymptotic analysis

At this stage both and are unspecified functions. They will be determined through the matching of the inner () and outer regions (), cf. Figure LABEL:fig:regions.

3.2 The outer region

To make the mechanics of the linearization and matching as transparent as possible we shall now change variables by linearizing around the south pole in the formulation (LABEL:LLG_3comp):

This recovers

\hb@xt@.01(3.8)
\hb@xt@.01(3.9)

To solve this we introduce , whence

\hb@xt@.01(3.10)

This is simply the projection of the flow from the sphere on to the tangent plane at the pole. Notice that in the respective limits we the recover modified linear heat equation () and Schrödinger equation ( on the tangent plane as appropriate.

To match the inner and outer regions we first define the similarity variables

\hb@xt@.01(3.11)

To reduce confusion in what follows we shall denote

Under this change of variables equation (LABEL:Zeq1) becomes

\hb@xt@.01(3.12)

A solution for this equation is

for any — this is just in (LABEL:Zeq1). This corresponds to the eigenfunction of the dominant eigenvalue of , which governs the generic long time behaviour of solutions of (LABEL:ZeqSimVars). When we allow to vary slowly with , we obtain a series expansion for the solution of the form

\hb@xt@.01(3.13)

We now see that the introduction of non-zero does not affect the procedure for the expansion.

Denoting , we introduce

\hb@xt@.01(3.14)

We recover and through , and expand for small :

\hb@xt@.01(3.15)
\hb@xt@.01(3.16)

Here and in what follows one should be slightly careful interpreting all formulae involving the because of multi-valuedness. For future reference we note that

\hb@xt@.01(3.17)
\hb@xt@.01(3.18)

3.3 The matching

In order to match the inner region to the outer we first write the inner solution in the similarity variables:

\hb@xt@.01(3.19)
\hb@xt@.01(3.20)

The matching procedure now involves setting and such that the expansions (LABEL:ThetaInAsy) and (LABEL:PhiInAsy) agree with (LABEL:ThetaOutAsy) and (LABEL:PhiOutAsy) respectively, to two orders in :

To solve this we set (with algebraic in ), and after rearranging terms we get

\hb@xt@.01(3.21)

Since is defined not to change exponentially fast in , we neglect the terms of . Using the definition of we may simplify (LABEL:FullSys1) to get

\hb@xt@.01(3.22)

Finally, we introduce , so that

\hb@xt@.01(3.23)

which is independent of and . This formulation strongly suggests that the case is not different from the dissipative case . Before solving and studying the system (LABEL:FullSys3) let us recall what its solutions tell us: gives an algebraic correction to the blowup rate, determines the local behaviour of near blowup and describes the amplitude and orientation of the solution in self-similar coordinates, see (LABEL:e:angle1),(LABEL:e:angle2). In order to fully understand blowup, we need to solve for the blowup coordinates and determine their stability.

The blowup solution is represented by for some constant (or , i.e. ), with . In particular, equations (LABEL:e:angle1),(LABEL:e:angle2) show that there is an angle between the sphere bubbling off and the solution remaining at/after blowup, see Figure LABEL:fig:angle.

Figure 3.2: The tangent plane at the south pole can be identified with the complex plane. The thick curve represents for a time near blowup. The angle between the solution in the inner scale (which bubbles off) and the remote scale (which remains) is indicated.

By rotating the sphere we may take without loss of generality, i.e.  and (note the sign), see (LABEL:FullSys3). The blowup dynamics is described by

\hb@xt@.01(3.24)

and . This equation has general solutions of the form

We can immediately set as this “instability” reflects shifts in the blowup time and hence is not a real instability — this is common to all blowup problems [25]. We note that so that indeed as , and . This implies that , hence for large (cf. Figure LABEL:fig:angle).

At this stage we have an asymptotic description of the blowup rate and its local structure. Unfortunately, we do not have enough information to determine stability. To more carefully understand the dynamics of this system we need to linearize about this leading order solution to find the subsequent corrections , , and in , , and , respectively. Taking , , , , we get as the system for the next order

\hb@xt@.01(3.25)

which separates into two systems. The first one is (using )

\hb@xt@.01(3.26)

which, using that solves (LABEL:e:sigma0), reduces to

the same equation as for , and provides no additional information. The other system is

\hb@xt@.01(3.27)

which can be rewritten as

or, again using that solves (LABEL:e:sigma0),

i.e., once again equation (LABEL:e:sigma0). The asymptotic behaviour of and is thus given by ()

where the exponentially growing terms show that blowup is unstable (the neutral mode corresponds to a (fixed, time-independent) rotation of the sphere).

4 Numerical computations

To supplement the formal analysis above we now present some numerical experiments in the radial, equivariant and fully two-dimensional cases.

4.1 Numerical methods

To reliably numerically simulate potentially singular solutions to (LABEL:LLG1) one needs to use adaptivity in both time and space as well satisfy the constraint . For the former, we use adaptive numerical methods as described in [10]. This approach is based on the moving mesh PDE approach of Huang and Russell [19] combined with scale-invariance and the Sundman transformation in time. The expository paper [10] provides many examples of this method being effective for computing blowup solutions to many different problems. For the latter we can either use formulation (LABEL:LLG_3comp) and use a projection step or regularize the Euler angle formulation (LABEL:LLG_2comp). We have implemented both and found little difference in efficiency or accuracy and hence will use formulation (LABEL:LLG_3comp) for all but Example 1 as it directly follows the above asymptotic analysis.

Full two-dimensional calculations have only been performed in the case of formulation (LABEL:LLG1) and on the unit disk. This latter fact is for numerical convenience and in no way affects the structure of local singularities (should they arise). Here adaptivity was performed using the parabolic Monge-Ampere equation as described in [9]

4.2 Numerical results

Figure 4.1: Left: Physical grid on which the solution was computed. Right: Computational grid in the region of . Notice there is a region of essentially constant in grid trajectories in this region. Some trajectories are leaving this region as and . Note, in both figures only every fifth grid trajectory is plotted.
Figure 4.2: Left: Solution on physical grid at selected times. Right: Solution on computational grid at same times. This clearly shows that the blowup region is very well resolved.

In this Section we present a sequence of numerical experiments to validate the results above. For examples 1-3 we used spatial points, the monitor function

and took