Spin Waves and Skyrmions in Magneto-Ferroelectric Superlattices: Theory and Simulation

: We present in this paper the effects of Dzyaloshinskii–Moriya (DM) magnetoelectric coupling between ferroelectric and magnetic layers in a superlattice formed by alternate magnetic and ferroelectric ﬁlms. Magnetic ﬁlms are ﬁlms of simple cubic lattice with Heisenberg spins interacting with each other via an exchange J and a DM interaction with the ferroelectric interface. Electrical polarizations of ± 1 are assigned at simple cubic lattice sites in the ferroelectric ﬁlms. We determine the ground-state (GS) spin conﬁguration in the magnetic ﬁlm. In zero ﬁeld, the GS is periodically non-collinear (helical structure) and in an applied ﬁeld H perpendicular to the layers, it shows the existence of skyrmions at the interface. Using the Green’s function method we study the spin waves (SW) excited in a monolayer and also in a bilayer sandwiched between ferroelectric ﬁlms, in zero ﬁeld. We show that the DM interaction strongly affects the long-wave length SW mode. We calculate also the magnetization at low temperatures. We use next Monte Carlo simulations to calculate various physical quantities at ﬁnite temperatures such as the critical temperature, the layer magnetization and the layer polarization, as functions of the magneto-electric DM coupling and the applied magnetic ﬁeld. Phase transition to the disordered phase is studied.

Multiferroics and superlattices of multiferroics (for example, PZT/LSMO and BTO/LSMO) are attracting increasing interest and provoking much research activity driven by the profound physics of these materials, the coexistence and coupling of ferroelectric and magnetic orders. Enormous attention was paid to studying the non-uniform states in magnetic/ferroelectric superlattices both theoretically [32] and experimentally [33][34][35].
The purpose of this paper is to show the main features of our results obtained for a magneto-ferroelectric superlattice where the interface coupling is a DM interaction. The effect of the frustration due to the competing next-nearest-neighbor (NNN) interaction is emphasized.
The paper is organized as follows. Section 2 is devoted to the description of our model and the determination of the ground-state spin configuration with and without applied magnetic field. Skyrmion crystal is shown with varying interface parameters. Section 3 shows the results for spin-wave spectrum obtained by the Green's function method. Results of Monte Carlo simulations for the phase transition in the system is shown in Section 4. These results show that the skyrmion crystal is stable at finite temperatures below a transition temperature T c . Concluding remarks are given in Section 5.

Model: Ground State
Consider a superlattice composed of alternate magnetic and ferroelectric films. Both have the structure of simple cubic lattice of the same lattice constant, for simplicity. The Hamiltonian of this multiferroic superlattice is expressed as: where H m and H f are the Hamiltonians of the magnetic and ferroelectric subsystems, respectively, H m f denotes Hamiltonian of magnetoelectric interaction at the interface of two adjacent films. We are interested in the frustrated regime. Therefore we describe the Hamiltonian of the magnetic film with the Heisenberg spin model on a simple cubic lattice as follows: where S i is the spin on the i-th site, H the external magnetic field, J m ij the magnetic interaction between two spins at i and j sites. We shall take into account both the nearest neighbors (NN) interaction, denoted by J m , and the next-nearest neighbor (NNN) interaction denoted by J 2m . We consider J m > 0 to be the same everywhere in the magnetic film. To introduce the frustration we shall consider J 2m < 0, namely antiferromagnetic interaction, the same everywhere. The external magnetic field H is applied along the z-axis which is perpendicular to the plane of the layers. The interaction of the spins at the interface will be given below.
For the ferroelectric film, we suppose for simplicity that electric polarizations are Ising-like vectors of magnitude 1, pointing in the ±z direction. The Hamiltonian is given by where P i is the polarization on the i-th lattice site, J f ij the interaction parameter between polarizations at sites i and j. Similar to the magnetic subsystem we will take the same J f ij = J f > 0 for all NN pairs, and J ij = J 2 f < 0 for NNN pairs.
Before introducing the DM interface interaction, let us emphasize that the bulk J 1 − J 2 model on the simple cubic lattice has been studied with Heisenberg spins [36] and the Ising model [37] where J 1 and J 2 are both antiferromagnetic (< 0). The critical value |J c 2 | = 0.25|J 1 | above which the bipartite antiferromagnetic ordering changes into a frustrated ordering where a line is with spins up, and its neighboring lines are with spins down. In the case of J m > 0 (ferro), and J 2m < 0, it is easy to show that the critical value where the ferromagnetic becomes antiferromagnetic is J 2m c = −0.5J m . Below this value, the antiferromagnetic ordering replaces the ferromagnetic ordering.
We know that the DM interaction is written as where S i is the spin of the i-th magnetic ion, and D i,j is the Dzyaloshinskii-Moriya vector. The vector D i,j is proportional to the vector product R × r i,j of the vector R which specifies the displacement of the ligand (for example, oxygen) and the unit vector r i,j along the axis connecting the magnetic ions i and j (see Figure 1). One then has We define thus where D is a constant, z the unit vector on the z axis, and e i,j = −e j,i = 1. For the magnetoelectric interaction at the interface, we choose the interface Hamiltonian following reference [32]: where P k is the polarization at the site k of the ferroelectric interface layer, while S i and S j belong to the interface magnetic layer. In this expression J m f i,j e i,j P k plays the role of the DM vector perpendicular to the xy plane, given by Equation (6). When summing the neighboring pairs (i, j), attention should be paid on the signs of e i,j and S i × S j (see example in reference [32]).
Hereafter, we suppose J m f i,j = J m f independent of (i, j). Since P k is in the z direction, i. e. the DM vector is in the z direction, in the absence of an applied field the spins in the magnetic layers will lie in the xy plane to minimize the interface interaction energy, according to Equation (7).
From Equation (7), we see that the magnetoelectric interaction J m f favors a canted (non collinear) spin structure. It competes with the exchange interactions J m and J 2m of H m which favor collinear (ferro and antiferro) spin configurations. In ferroelectric layers, there is just ferro-or antiferromagnetic ordering due to the Ising nature. Usually the magnetic or ferroelectric exchange interaction is the leading term in the Hamiltonian, so that in many situations the magnetoelectric effect is negligible. However, in nanofilms of superlattices the magnetoelectric interaction is crucial for the creation of non-collinear long-range spin order. It has been shown that Rashba spin-orbit coupling can lead to a strong DM interaction at the interface [38,39], where the broken inversion symmetry at the interface can change the magnetic states.
Since the polarizations are along the z axis, the interface DM interaction is minimum when S i and S j lie in the xy interface plane and perpendicular to each other. However the collinear exchange interactions among the spins will compete with the DM perpendicular configuration. The resulting configuration is non-collinear. In the general case where a magnetic field is applied and there is also the frustration, we use the steepest descent method [31,32] to determine the ground state.
We show in Figure 2 an example where H = 0 with no NNN interaction. The magnetic film has a single layer. The GS spin configurations have periodic structures. Several remarks are in order: (i) Each spin has the same turning angle θ with its NN in both x and y direction. The schematic zoom in Figure 2c shows that the spins on the same diagonal (spins 2 and 3, spins 4 and 5) are parallel. This explains the structures shown in Figure 2a   At this field strength H = 0.25, if we increase the frustration, for example, J 2m = J 2 f = −0.3, then the skyrmion structure is enhanced: we can observe a clear 3D skyrmion crystal structure not only in the interface layer but also in the interior layers (not shown).
Many other cases will be shown in a full paper.

Spin Waves
Here let us show theoretically spin-waves (SW) excited in the magnetic film in zero field, in some simple cases. The method we employ is the Green's function technique for non-collinear spin configurations which has been shown to be efficient for studying low-T properties of quantum spin systems such as helimagnets [40] and systems with a DM interaction [41]. This technique is rather lengthy to describe here. The reader is referred to these papers for mathematical details.
Let us show here only the resulting spin-wave energy where the reduced anisotropy is d = K/J m and γ = (cos k x a + cos k y a)/2, k x and k y being the wave-vector components in the xy planes, a the lattice constant, θ the angle between two NN spins given by θ = arctan(− 2J m f J m ). We show in Figure 4 the spin-wave energy E versus the wave vector k = k x = k y , for weak and strong values of the DM interaction strength D. For the weak value of θ (Figure 4a), namely weak J m f , the long-wave length mode energy (k → 0) is proportional to k 2 as in ferromagnets, while for strong θ (Figure 4b), E is proportional to k as in antiferromagnets. The effect of the DM interaction is thus very strong on the spin-wave behavior.
The case of a bilayer as well as the effects of other parameters on the spin-wave spectrum will be shown in the full paper.

Phase Transition
We have seen the skyrmion crystal in Figure 3 at T = 0. Shall this structure survive at finite T? To answer this question we have performed Monte Carlo simulations using the Metropolis algorithm [42]. We have to define an order parameter for the skyrmion crystal. In complicated spin orderings such as spin glasses, we have to follow each individual spin during the time. If it is frozen then its time average is not zero. This is called Edwards-Anderson order parameter for spin glasses [43].
For the magnetic layers, the definition of an order parameter for a skyrmion crystal is not obvious. Taking advantage of the fact that we know the GS, we define the order parameter as the projection of an actual spin configuration at a given T on its GS and we take the time average. This order parameter of layer n is thus defined as where S i (T, t) is the i-th spin at the time t, at temperature T, and S i (T = 0) is its state in the GS. The order parameter M m (n) is close to 1 at very low T where each spin is only weakly deviated from its state in the GS. M m (n) is zero when every spin strongly fluctuates in the paramagnetic state. For the ferroelectric layers, the order parameter M f (n) of layer n is defined as the magnetization where ... denotes the time average. The total order parameters M m and M f are the sum of the layer order parameters, namely M m = ∑ n M m (n) and M f = ∑ n M f (n). In Figure 5 we show the magnetic energy and magnetic order parameter versus temperature in an external magnetic field, for various sets of NNN interaction. Note that the phase transition occurs at the energy curvature changes, namely at the maximum of the specific heat. The red curve in Figure 5a

Conclusions
We have briefly shown in this paper some striking features of the properties of a magnetoferroelectric superlattice. We have taken into account a DM interaction at the interface and the frustration due to the NNN interactions. The skyrmion crystal is shown to be stable at finite temperatures. This may have interesting applications in transport by skyrmions.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.