Lattice solutions in a Ginzburg-Landau model for a chiral magnet

We examine micromagnetic pattern formation in chiral magnets, driven by the competition of Heisenberg exchange, Dzyaloshinskii-Moriya interaction, easy-plane anisotropy and thermodynamic Landau potentials. Based on equivariant bifurcation theory we prove existence of lattice solutions branching off the zero magnetization state and investigate their stability. We observe in particular the stabilization of quadratic vortex-antivortex lattice configurations and instability of hexagonal skyrmion lattice configurations, and we illustrate our findings by numerical studies.


Introduction and main results
Dzyaloshinskii-Moriya interaction (DMI) is the antisymmetric counterpart of Heisenberg exchange. It arises from the lack of inversion symmetry in certain magnetic system, induced by the underlying crystal structures or by the system geometry in the presence of interfaces. In mircromagnetic models DMI arises in form of linear combinations of the so-called Lifshitz invariants, i.e. the components of the chirality tensor ∇m × m, and is therefore sensitive with respect to reflections and independent rotation in the domain and the target space, respectively. It is well-known that DMI gives rise to modulated phases. The basic phenomenon is that the energy of the homogeneous magnetization state m = const. can be reduced below zero by means of spiralization in form of periodic domain wall arrays, the helical phase. A prominent form of doubly periodic lattice state arises from the stabilization of topological structures, so-called chiral skyrmions, in two-dimensional chiral ferromagnets. Chiral skyrmions are localized structures which are topologically characterized by a unit S 2 degree and a well-defined helicity depending on the specific form of DMI. In the presence of sufficiently strong perpendicular anisotropy and/or Zeeman field interaction, chiral skyrmions occur as local energy minimiziers in form of isolated topological solitons. Zeeman field interaction enables the possibility of an intermediate regime where chiral skyrmion embedded into a hexagonal lattice are expected to be globally energy minimizing. Micromagnetic theories including DMI have been proposed in [11] with the idea that skyrmion lattice configurations represent a magnetic analogue of the mixed state in type-II superconductors. Corresponding phase diagrams and stability questions have been examined analytically and numerically in the seminal work [9,10], see [22] for a recent review. A fully rigorous functional analytic theory on the existence, stability, asymptotics, internal structure and exact solvability of isolated chiral skyrmions has recently started to emerge [25,16,23,20,21,7,8]. Lattice states in two-dimensional Ginzburg-Landau models for chiral magnets including thermodynamic effects have been proposed theoretically in [29], see also [26,34]. Mathematically one may expect a close analogy with Abrikosov's vortex lattice solutions in Ginzburg-Landau models for superconductors [1] with a well-established theory in mathematical analysis. The occurrence of Abrikosov lattices in the framework of gauge-periodic solutions of Ginzburg-Landau equations and the optimality of hexagonal lattices have been thoroughly investigated by means of variational methods and bifurcation theory [6,28,32,5,2,30].
Ginzburg-Landau models in micromagnetics are, in contrast to superconductivity, directly formulated in terms of the physically observable magnetization field. A class of such models has been proposed and examined computationally in [29,26]. Compared to the purely ferromagnetic case |m| = const., Ginzburg-Landau models offer a larger variety of patterns, including vortex and half-skyrmion arrays of opposite in-plane winding on square lattices, and skyrmions on hexagonal lattices, see Figure 1.
We shall examine the occurrence of periodic solutions near the paramagnetic state of zero magnetization, i.e., in a high temperature regime. Our starting point are magnetization fields m : R 2 → R 3 governed by energy densities of The Dirichlet term with A > 0 is referred to as Heisenberg exchange interaction. The helicity term with D = 0 is a prototypical form of DMI arising in the context of noncentrosymmetric cubic crystals. Note that for m = (m, m 3 ) defined on a two-dimensional domain A key analytical feature induced by DMI is the loss of independent rotational symmetry in magnetization space R 3 and the domain R 2 . Finally, the Landau term f is an even polynomial in the modulus |m| where T − T C is the deviation from the Curie temperature. We shall focus on a minimal model with For stability reasons we also include the easy-plane anisotropy K(m ·ê 3 ) 2 with K > 0, which typically emerges as a reduced form of magnetostatic stray-field interaction in thin-film geometries, see e.g. [8,15,24].
We are interested in magnetization fields m which are periodic with respect to a two-dimensional lattice rΛ (r > 0) with where τ is a complex number in the fundamental domain of the modular group, referred to as the lattice shape parameter, see Section 2.1. Rescaling space we may assume r = 1. The rescaled energy density reads The energy density is invariant with respect to translations and joint rota- . Since a sign reversal of the DM density can be achieved by reflections such as m 3 → −m 3 , we may assume w.l.o.g that κ > 0.
We shall investigate the occurrence and stability of non-trivial periodic solutions m to the Euler-Lagrange equation associated to (1) We first discuss energy minimizing solutions.
(i) If λ > κ 2 , then m ≡ 0 is the unique energy minimizer on every lattice.
(ii) If λ < κ 2 and β = 0, then the helix (see Figure 2) is up to a joint rotation the unique energy minimizer on suitable lattices.
In the isotropic case β = 0, Theorem 1 indicates the existence of only two phases, paramagnetic or helical, while the picture in the anisotropic case β > 0 is incomplete. The occurrence of helical phases is common to other mathematically related theories for condensed matter such as the Oseen-Frank model for chiral liquid crystals (see e.g. [33]) or the Gross-Pitaevskii model for spin-orbit coupled Bose-Einstein condensates (see e.g. [3]). Helical structures in chiral ferromagnets are also discussed in [14,27].
Here we are interested in doubly periodic solutions. Given a lattice Λ, we aim to find λ and a non-trivial Λ-periodic solution m of (3) at λ. We call such pairs (m, λ) Λ-lattice solutions and will prove the following: has a branch of Λ-lattice solution (m s , λ s ) in a neighbourhood of m 0 ≡ 0 and satisfying the non-resonances condition for any ω in the dual lattice Λ * , see Sec. 2.3, with |ω| = 1 Im τ . The branch (m s , λ s ) is analytically parameterized by a real parameter s near 0 and has the form as s → 0 with λ 2 < 0 and ϕ 1 explicitly determined.
The morphology of bifurcation solutions is related to symmetry properties of the underlying lattice. Depending on this, the first-order bifurcation solution ϕ 1 , arising from the first critical wave number, indicates a threefold pattern formation.
1. helical pattern: on non-equilateral lattices, the first-order bifurcation solution (26) is given by a single helical mode (see Figure 2); 2. vortex-antivortex pattern: on square or rhombic lattices, the first-order bifurcation solution (27) is a superposition of two helices propagating in different directions, see Figure 1(a); 3. skyrmionic pattern: on hexagonal lattices, the first-order bifurcation solution (28) is a superposition of three helices propagating in distinct directions, see Figure 1  Reducing the domain of the bifurcation parameter s if necessary we shall prove the following stability properties of the bifurcation solutions obtained in Theorem 2: Theorem 3. Let α > 0, β ≥ 0 and κ > 0 satisfying and where γ is the quotient of the second and first critical wave number, i.e., depending on the lattice shape τ = |τ |e iθ .
(iii) The bifurcation solution on hexagonal lattices is unstable for any β ≥ 0.
Corollary 1. The quadratic vortex-antivortex lattice configuration exists and is stable if The admissible set of (κ, β) is not empty, see Figure 3. The existence and stability results are only an initial step towards understanding the stabilization of two-dimensional lattice solutions in chiral magnets. In particular, stability of lattice solutions is only examined under the simplest perturbations which preserve lattice periodicity. A more general stability result in the style of [31] is beyond the scope of this work and requires a different approach.
The mathematical framework for our construction of lattice solutions is the equivariant branching lemma [12,17], a concept of symmetry-breaking bifurcation based on a particular type of (axial) symmetry group. For each lattice we identify an isotropy subgroup Σ ⊂ SO(3) so that the fixed subspace of Σ in the kernel of linearized operator D m F (m 0 , λ 0 ) is one-dimensional. By means of an equivariant Lyapunov-Schmidt procedure, (3) reduces to a onedimensional bifurcation equation. The implicit function theorem provides a solutions to the bifurcation equation in the one-dimensional fixed subspace of Σ, from which a solution to (3) can be reconstructed. This solution is the bifurcation solution and inherits the symmetries featured by Σ.
In Section 2 we shall briefly recall the representation of lattices in the plane with an emphasis on symmetry and Fourier series which are key to our bifurcation argument. In Section 3 we shall derive energy bounds proving Theorem 1. Solving the linearized equation of (3) explicitly by Fourier methods is the key ingredient to the proof of Theorem 2 in Section 4 and provides insight about the morphology and topology of bifurcation solutions. In Section 5 we investigate the stability of bifurcation solutions under Λperiodic perturbations proving Theorem 3 . Finally in Section 6 we validate our analytical results by a series of numerical simulations of gradient flows using a modified Crank-Nicolson scheme.

Representation of lattices
Recall that a planar lattice Λ is the integer span of two linearly independent vectors t 1 , t 2 ∈ R 2 , i.e., Given x ∈ R 2 , a primitive cell of Λ is a set of the form The lattice basis {t 1 , t 2 } is clearly non-unique. Identifying R 2 ∼ = C, however, the complex ratio τ ∈ C of two basis vectors of Λ contained in the fundamental domain parametrizes the lattice shape uniquely, see e.g. [4]. Writing τ = |τ |e iθ , the range of |τ | and θ corresponding to fundamental domain (9) is A lattice is called equilateral if |τ | = 1, where the borderline cases θ = π/2 and θ = π/3 are referred to as square and hexagonal lattice, respectively; other equilateral lattices are called rhombic.
There are two distinct types of symmetries preserving the lattice: the translation R 2 /Λ and the holohedry, which is a finite subgroup of O(2). For planar lattices there are four basic holohedries: non-equilateral lattices have olohedry D 1 (oblique) or D 2 (rectangular); rhombic lattices have holohedry D 2 ; square lattices have holohedry D 4 ; hexagonal lattices have holohedry D 6 , see e.g. [12].

Function spaces on lattices
As the governing energy densities and the Euler-Lagrange equations are invariant under translation and joint rotation we can fix one basis vector of the lattice Λ as 2πrê 1 so that Λ = 2πr(Z ⊕ τ Z), is uniquely characterized by r and τ . Upon rescaling we can arrange Λ = 2π(Z ⊕ τ Z) spanned by {2πê 1 , 2πτ } and consider the rescaled density (1) containing only dimensionless parameters.
For Λ-periodic functions or fields f, g on R 2 we denote the average by the L 2 scalar product by and the L 2 norm by f := f, f once existent. Accordingly we define , respectively. Thanks to Sobolev embedding the average energy (11) defines an analytic functional on

Dual lattice and Fourier series
Fourier expansion on Λ requires the notion of dual lattice given by In particular Λ * = A −T Z 2 for Λ = 2πAZ 2 where in our setting In the equilateral case |τ | = 1, dual lattices remain square if θ = π/2 and hexagonal if θ = π/3. For f ∈ L 2 Λ and ν ∈ Λ * Fourier coefficients are defined as and the following Fourier expansion holds true in the L 2 sense along with Parseval's identity

Equivariance and lattice symmetry
We consider the following action of σ ∈ SO(3) on a field m : under which (12) is equivariant in the sense that In a periodic setting we need to restrict to certain discrete subgroups Σ of the holohedry that preserve the lattice Λ = 2π(Z ⊕ τ Z) and suitable subspaces of the kernel of the linearized operator D m F (0, λ 0 ). Three cases are to be distinguished: For |τ | > 1 we consider the symmetry group consisting of the identity and rotation through π where Rhombic and square lattices (|τ | = 1 and π 3 < θ ≤ π 2 ) have an additional symmetry given by reflections across the diagonals of the lattice cell. Therefore, in this case, we consider the symmetry group where For hexagonal lattices |τ | = 1 and θ = π 3 , the symmetry group considered is the subgroup of D 6 containing the rotations by an angle of kπ where In each of the three cases Σ = Σ i , the group action is an isometry on H k Λ for every k ∈ N 0 , and the operator (12) is Σ-equivariant in the sense that Λ orthogonal projections on such invariant subspaces of H 2 Λ are equivariant. Equivariance therefore propagates to the Lyapunov-Schmidt decomposition enabling a reduction of the bifurcation equation to Fix X 0 (Σ) where X 0 is the kernel of the linearization of F at a bifurcation point (0, λ 0 ), see e.g. [12]. The symmetry of a field φ ∈ H 2 Λ is given in terms of the isotropy subgroup i.e., the largest subgroup of SO(3) which fixes φ. Bifurcation solutions arising from the equivariant branching lemma turn out to have full Σ symmetry and are unique in this class.

Energy bounds on lattices
For a lattice Λ and e(m) given by (1), we examine ansatz-free lower bounds We start by expressing the total exchange density e 0 (m) = 1 2 as a sum of sign definite terms and a null Lagrangian also know as Frank's formula in the theory of liquid crystals, see e.g. [33] Chapter 3.
in the sense of distributions. In particular for m ∈ H 1 with equality if and only if ∇ × m + κm = 0 and β = 0.
From the lemma we obtain immediately: Proposition 1. If λ > κ 2 , then m ≡ 0 is the unique energy minimizer on every lattice. If λ < κ 2 the energy admits a lower bound which for β = 0 is precisely attained for m of constant modulus |m| = κ 2 − λ α and such that ∇ × m + κm = 0.
The unimodular Beltrami fields being parallel to their curl have been classified by Ericksen within the variational theory of liquid crystals, see e.g. [33] and references therein. For the present case of constant κ those are helices of pitch 2π/|κ|, i.e., (4). For the convenience of the reader we present this fundamental result in the Appendix Lemma 6. Thus in order to realize the lower energy bound, the underlying lattice Λ is required to accommodate such a helix. In this case the zero state losses its linear stability at λ = κ 2 .
In fact, the Hessian By the preceding arguments it satisfies and has a helical instability at λ = κ 2 .
A straightforward analysis of the minimizing problem (Im τ k 1 ) 2 + (Re τ k 1 − k 2 ) 2 yields (see Figure 5): Hence, in this case, we shall consider the bifurcation occurring at satisfying (6), which guarantees that only the first non-trivial wave number contributes to the kernel of the linearization.

Proposition 2. Under the assumptions of Theorem 2, the linearization
is a self-adjoint Fredholm operator with dim ker L 0 = N , where N = 2 on non-equilateral lattices, N = 4 on rhombic and square lattices and N = 6 on hexagonal lattices. The fixed subspace of the corresponding symmetry group Σ (see section 2.4) in X 0 = ker L 0 is one-dimensional, Proof. The normalized Fourier coefficients (24) obtained in the proof of Lemma 2 are The corresponding real-space solutions of (22) to the wave vector ν are We need to find all solutions of L 0 φ = 0 for a given wave number |ν|. According to Lemma 3, we consider following cases separately.
Inserting the analytic expansion of bifurcation solutions into (18) and (17) we obtain (19)- (21). Explicit calculations are carried out in Appendix B.
Topology of bifurcation solutions On hexagonal lattices, the bifurcation solution is nowhere vanishing and features in every primitive cell a skyrmion, i.e. a vortex-like structure with the magnetizations pointing upwards at the core and downwards at the perimeter, see Figure 6(a).
On rhombic and square lattices, the horizontal component of ϕ 1 has a finite number of isolated zeros in a primitive cell and forms a vortex or an antivortex around each zero. The antivortices are half-skyrmions (sometimes referred to as merons, see e.g. [34]) and have magnetizations pointing upwards or downwards at the core; while the center of vortices are singularity points (m = 0) due to the continuity and the Σ-invariance of bifurcation solutions, see Figure 6(b).

Linear stability of the bifurcation solutions
In this section we discuss the stability of bifurcation solutions following the standard perturbation argument, see e.g. [19]. Suppose (m s , λ s ) is a We focus on the larger root as positivity turns out to be necessary for the stability of bifurcation solutions. The linearization of (3) at (m s , λ s ) is given by We first investigate the spectrum of L 0 .
Lemma 4. The spectrum of L 0 is discrete and non-negative if (7) and (8) hold. Otherwise L 0 has negative eigenvalues.
From now on we focus on the case λ 0 > 0. It follows from Lemma 4 and the standard perturbation theory of eigenvalue [18] that the spectrum of L s consists of eigenvalues of the same multiplicities in an neighbourhood of the eigenvalues of L 0 . Thus the stability of (m s , λ s ) depends on the perturbation of the critical eigenvalue 0. It follows from Proposition 2 that there exist the following topological decompositions where X 1 = {φ ∈ H 2 Λ : φ ⊥ ker L 0 in L 2 }. We first consider the perturbation of the zero eigenvalue corresponding to the eigenvector ϕ 1 spanning Fix X 0 (Σ).
It remains to examine the properties of µ(s). Differentiating (29) with respect to s in s = 0 yieldṡ Calculating the second derivative of µ(s) at s = 0, we obtain
The helical solution on non-equilateral lattices is stable.
Proposition 3. For β ≥ 0 and every nonzero s ∈ (−δ, δ), the bifurcation solution on non-equilateral lattices is linearly stable in the sense that L s ≥ 0 with ker L s = span{∂ 2 m s }.
The hexagonal skyrmion lattice is unstable independently of any additional easy-plane anisotropy: For example, for any nonzero s ∈ (−δ, δ) We aim to find equilibria of the energy functional E Λ by solving (32) numerically on the primitive cell Ω τ induced by lattice spanned by {2πê 1 , 2πτ }. Eq. (32) is discretized by a modified Crank-Nicolson approximation for the time variable and a Fourier collocation method for the space variable. We denote m N the trigonometric interpolation function of m on the discretized grid by N 2 collocation points for N N = {0, . . . , N − 1}, N ∈ N and odd. For continuous fields u, v on Λ τ , we define the discrete L 2 scalar product the associated norm · N , defined by u 2 N := u, u N and the discrete energy where I N denotes the trigonometric interpolation operator and At each time iteration, in order to find the solution m n+1 N , we use a fixed point iteration for some time step ∆t. Well-posedness and a-priori error bounds of this numerical scheme follow analogously as in [13]. Given initial data m 0 N , the corresponding sequence (m n N ) n∈N 0 satisfies the following energy law We have implemented this numerical scheme in MATLAB. At each timestep the iteration process stops if a certain norm of the difference of two successive iterations becomes smaller than a chosen stopping tolerance. In our case we choose the L ∞ -norm and set the stopping tolerance to 10 −8 . The discrete energy is evaluated in each time step, and the terminal time is controlled through a smallness condition for the discrete energy gradient, i.e.
After the termination of the scheme, an equilibrium configuration is reached approximately.

Numerical experiments
We have implemented the method on a lattice of 275 × 275 grid points and with a time increment ∆t = 0.1 for different parameters and a randomly distributed initial field with modulus between 0 and 0.1 as initial condition.
Bifurcation without anisotropy We first performed the simulations at β = 0 on equilateral lattices , i.e. |τ | = 1 with different angle θ between π 3 and π 2 by setting λ = λ 0 + δλ 2 , where δ = 0.001, λ 0 and λ 2 are calculated according to (7) and (20), respectively. Despite our instability results at β = 0 (see Section 5) doubly periodic solutions emerge numerically. As displayed in Figure 7, vortex-antivortex lattice configurations were observed for different θ ∈ ( π 3 , π 2 ], while the skyrmion lattice was observed for θ = π 3 . Then we implemented the method on lattices |τ |e i π 2 with different |τ | ≥ 1. For |τ | ≈ 1 we obtained the vortex-antivortex lattice configuration, indicating its stability under perturbations of the lattice shape, while for |τ | large enough only the helix state as shown in Figure 8 was observable. Bifurcation with easy-plane anisotropy Next we included easy-plane anisotropy and implemented the simulations on a square lattice for different parameters κ and β near the bifurcation point λ 0 . For κ ∈ (0, 0.5), the bifurcation point λ 0 is negative for any β ≥ 0. In this case, we obtained an almost homogeneous field.
When the value of κ was increased over 0.5, vortex-antivortex lattice configurations were observed for β ≥ 0 small enough so that λ 0 > 0. For β large enough so that λ 0 < 0, the vortex-antivortex lattice configuration decayed to an almost homogenous field, as shown in Figure 9. Other patterns emerged for κ > 1.2 and small β ≥ 0. For example, at κ = 1.4 and small β ≥ 0 the solution converged to a stripe pattern, i.e. helices with a pitch smaller than 2π. Vortex-antivortex lattice configurations were observed for β larger than the stability threshold (about 2.3κ) and λ 0 > 0. Increasing β further so that λ 0 < 0, we obtained the almost homogeneous field again, as shown in Figure 10.

Acknowledgements
This work is support by Deutsche Forschungsgemeinschaft (DFG grant no. ME 2273/3-1). Proof. Upon rescaling one may assume M = 1 and κ = 1. Taking the divergence it follows that ∇ · m = 0 and hence ∆m + m = 0 in the sense of distributions by taking the curl. So m is smooth by virtue of standard elliptic regularity theory. In particular, it is enough to prove the claim locally. Denoting the componentwise gradient by (∇m) jk := (∂ j m k ) we claim that rank(∇m) = 1 and (∇m) 2 = 0.
Fixing a point x we may assume after rotation m(x) =ê 3 , so that by (35) the matrix ∇m(x) is given by a 2 × 2 matrix A such that tr A = 0 and tr A 2 = 0 by (36). It is easy to see that A 2 = 0 which implies det A = 0.
According to (33) there exist local smooth unit vector fields X and Y and a function λ such that ∇m = λX ⊗ Y and m = X × Y . Now ∇ × m = λ X × Y = λ m, hence λ = −1 and ∇m = −X ⊗ Y .
Assuming X is constant so that after rotation X =ê 1 , it follows that m = m(x 1 ) and m 1 = 0. Now the equation implies for the remaining components m 2 = −m 3 and m 3 = m 2 , so that after a rotation around thê e 1 axis, m = h.
To show that X = const. one may invoke the spectral theorem. In fact, symmetry of ∇X follows from the symmetry of ∇ 2 m = (∇ ⊗ ∇)m since −∇ 2 m = ∇X ⊗ Y + X ⊗ ∇Y so that − ∇ 2 m · Y = ∇X.

B Higher order terms in the bifurcation solution
We have the following analytic expansion for the bifurcation solution Inserted into (18), we obtain a hierarchy of equations in orders of s that can be solved successively. Solutions, which are fixed under the symmetry group, can be found by the Fourier method used for solving (22).