Skyrmion Crystal and Phase Transition in Magneto-Ferroelectric Superlattices: Dzyaloshinskii-Moriya Interaction in a Frustrated $J_1-J_2$ Model

The formation of a skyrmion crystal and its phase transition are studied taking into account the Dzyaloshinskii-Moriya (DM) interaction at the interface between a ferroelectric layer and a magnetic layer in a superlattice. The frustration is introduced in both magnetic and ferroelectric films. The films have the simple cubic lattice structure. Spins inside magnetic layers are Heisenberg spins interacting with each other via nearest-neighbor (NN) exchange $J^m$ and next-nearest-neighbor (NNN) exchange $J^{2m}$. Polarizations in the ferroelectric layers are assumed to be of Ising type with NN and NNN interactions $J^f$ and $J^{2f}$. At the magnetoelectric interface, a DM interaction $J^{mf}$ between spins and polarizations is supposed. The spin configuration in the ground state is calculated by the steepest descent method. In an applied magnetic field $\mathbf H$ perpendicular to the layers, we show that the formation of skyrmions at the magnetoelectric interface is strongly enhanced by the frustration brought about by the NNN antiferromagnetic interactions $J^{2m}$ and $J^{2f}$. Various physical quantities at finite temperatures are obtained by Monte Carlo simulations. We show the critical temperature, the order parameters of magnetic and ferroelectric layers as functions of the interface DM coupling, the applied magnetic field and $J^{2m}$ and $J^{2f}$. The phase transition to the disordered phase is studied in details.\\ Keywords: kyrmions; phase transition; frustration; superlattice; magneto-ferroelectric coupling; Dzyaloshinskii-Moriya interaction; $J_1-J_2$ model; Monte Carlo simulation.

Real magnetic materials have complicated interactions; there are large families of frustrated systems, such as the heavy lanthanides metals (holmium, terbium, and dysprosium) [34,35] and helical MnSi [36]. Other interesting properties of skyrmions in magnetically frustrated systems have also been investigated [31,[37][38][39][40][41][42][43][44][45][46]. Multiferroics and superlattices of multiferroics (for example, PZT/LSMO and BTO/LSMO) currently attract many research activities on these materials because of the coexistence of ferroelectric and magnetic orderings. A large number of investigations was devoted to the non-uniform states in magneto-ferroelectric superlattices both theoretically [47] and experimentally [48][49][50][51][52][53][54]. In Ref. [55,56], Janssen et al. proposed a new model for the interaction between polarizations and spins in magneto-ferroelectric superlattices. Using this model, the dynamics and configurations of domain walls for the unidimensional case were simulated. Li et al. [57] proposed an algorithm based on the Monte Carlo method for a two-dimensional (2D) lattice. Recently, magneto-ferroelectric superlattices have drawn as much attention as magneto-electric (ME) materials. This is due to the intrinsic magneto-electric effects stemming from the spin-orbit interaction [52] as well as from the spin charge-orbital coupling [58]. It has been shown that surface ME effects appear due to the charge and spin accumulation [53,59]. The enhancement of magnetoelectric effects due to phase separation was shown in Ref. [60]. Many microscopic interaction mechanisms for different materials have been suggested. Among these, we can mention the lone skyrmion pair mechanism [61], the ferroelectricity in manganites [62], the multiferroicity induced by the spiral spin ordering [63], and the off-center shifts in geometrically frustrated magnets [64]. In Ref. [65], it was shown that magnetic frustration results in a phase competition, which is the origin of the magnetoelectric response. A possible experimental realization of an isolated skyrmion as well as a skyrmion lattice has been suggested [66]. In Ref. [15], the 2D Heisenberg exchange model with Dzyaloshinskii-Moriya (DM) interaction was used to study the lifetimes of antiferromagnetic skyrmions. Spin waves and skyrmions in magnetoelectric superlattices with a DM interface coupling have been studied [47]. Yadav et al. [67] have produced complex topologies of electrical polarizations-namely, nanometer-scale vortex-antivortex structures-using the competition between charge, orbital, and lattice degrees of freedom in superlattices of alternate lead titanate and strontium titanate films. They showed that the vortex state is the low-energy state for various superlattice periods. In Ref. [19], the authors explored skyrmions and antiskyrmions in a 2D frustrated ferromagnetic system with competing exchange interactions based on the J 1 − J 2 classical Heisenberg model on a simple square lattice [38] with dipole-dipole interaction, which was neglected in previous works [37,68]. Dipole-dipole interaction has been shown to cause frustration in systems of skyrmions. The interface-induced skyrmions have been investigated. The superstructures are obviously the best for creating interactions of skyrmions on different interfaces causing very particular dynamics compared to interactions between skyrmions of the same interface. We can mention a theoretical study of two skyrmions on two-layer systems [69] using micromagnetic modeling and an analysis based on the Thiele equation. It has been found that there is a reaction between them, such as the collision and bound state formation. The dynamics of these skyrmions depends on the sign of the DM interaction; the number of two skyrmions is well described by the Thiele equation. Another interesting aspect is a colossal spin-transfer torque effect of bound skyrmion pairs, revealed in antiferromagnetically coupled bilayer systems. Note that the study of current-induced motion using the Thiele equation was carried out in a skyrmion lattice using two models of the pinning potential [70]. Shi-Zeng Lin et al. [71] studied the dynamics of skyrmions in chiral magnets in the presence of a spin-polarized current. They also showed that skyrmions can be created by increasing the current in the magnetic spiral state. Numerical simulations of the Landau-Lifshitz-Gilbert equation were performed in Ref. [72], which reveals a remarkably robust and universal current-velocity relation of the skyrmion motion driven by the spin-transfer torque. This is unaffected by impurities and non-adiabatic effects, in sharp contrast to the case of the domain wall or spin helix.
Note that in Ref. [47], we studied the effects of Dzyaloshinskii-Moriya (DM) magneto-ferroelectric coupling in a superlattice formed by "unfrustrated" magnetic and ferroelectric films. In a zero field, we showed that the ground state (GS) spin configuration is periodically non-collinear. Through the use of a Green's function technique, we calculated the spin wave spectrum in a monolayer and in a bilayer sandwiched between ferroelectric films. In particular, we showed that the DM interaction strongly affects the long-wavelength spin wave mode. In a magnetic field H applied in the z direction perpendicular to the layers, we showed that skyrmions are arranged to form a crystalline structure at the interface.
In this paper, we study a superlattice composed of alternating "frustrated" magnetic films and "frustrated" ferroelectric films. The frustration due to competing interactions has been extensively investigated over the last four decades. The reader is referred to Ref. [73] for reviews on theories, simulations, and experiments on various frustrated systems. Our present aim is to investigate the effect of the frustration in the presence of the DM interaction at the magnetoelectric interface. It turns out that the frustration gives rise to an enhancement of skyrmions created by the DM interaction in a field H applied perpendicularly to the films. Monte Carlo simulations are carried out to study the skyrmion phase transition in the superlattice as a function of the frustration strength.
The paper is presented as follows. In Section 2, we describe our model and show how to determine the ground-state spin structure. The skyrmion crystal is shown with varying frustration parameters. Section 3 is devoted to the presentation of the Monte Carlo results for the phase transition in the system as a function of the frustration, in the presence of the interface DM coupling. These results show that the skyrmion crystalline structure is stable up to a transition temperature T c . Concluding remarks are given in Section 4.

Model
The superlattice that we study here is shown in Figure 1a. It is composed of magnetic and ferroelectric films with a simple cubic lattice. The Hamiltonian of this magneto-ferroelectric superlattice is given by where H m and H f indicate the magnetic and ferroelectric parts, respectively. H m f denotes the Hamiltonian of magnetoelectric interaction at the interface of two adjacent films. We are interested in the frustrated regime. Therefore, we consider the following magnetic Hamiltonian with the Heisenberg spin model: where S i is the spin occupying the i-th lattice site, H denotes the magnetic field applied along the z axis, and J m ij denotes the magnetic interaction between two spins S i and S j . We shall take into account both the nearest neighbor (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-as the same everywhere. The interactions between spins and polarizations at the interface are given below.
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 [74] and the Ising model [75], where J 1 and J 2 are both antiferromagnetic (< 0). The critical value |J c 2 | = 0.25|J 1 |; above this, the bipartite antiferromagnetic ordering changes into a frustrated ordering, where a line of spins up is surrounded by its neighboring lines of spins down. In the cases 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.
For the ferroelectric film, we assume that electric polarizations are of an Ising model of magnitude 1, pointing in the ±z direction. The Hamiltonian reads where P i is the polarization on the i-th lattice site and J f ij is the interaction parameter between polarizations at sites i and j. Similarly to with 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.
We know that the DM interaction is written as where S i is the spin at the i-th magnetic site, while D i,j denotes the Dzyaloshinskii-Moriya vector which is defined as R × r i,j , where R is the displacement of the ligand ion (oxygen) and r i,j is the unit vector along the axis connecting S i and S j (see Figure 1b). One then has We define where D is a constant, z is 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 Ref. [47]: 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 (see Figure 1c). In this expression, J m f i,j e i,j P k is defined as the DM vector which is along the z axis, given by Equation (6). When summing the neighboring pairs (i, j), attention should be paid to the signs of e i,j and S i × S j (see example in Ref. [47]).
, 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).
In Equation (7), the magnetoelectric interaction J m f favors a non-collinear spin structure in competition with the exchange interactions J m and J 2m , which favor collinear (ferro and antiferro) spin configurations. In ferroelectric layers, only collinear, ferro-, or antiferromagnetic ordering is possible because of the assumed Ising model for the polarizations. In historical demonstrations, the DM interaction was supposed to be small with respect to the exchange terms in the Hamiltonian. However, in superlattices, the magnetoelectric interaction is necessary to create non-collinear spin ordering. It has been shown that Rashba spin-orbit coupling can lead to a strong DM interaction at the interface [76,77], where the broken inversion symmetry at the interface can change the magnetic states. The DM interaction has been identified as a key ingredient in the creation and stabilization of skyrmions and chiral domain walls.

Ground State
From Equation (7), we see that the interface interaction is at its minimum when S i and S j lie in the xy interface plane and are perpendicular to each other in the absence of exchange interactions. When the exchange interactions are turned on, the collinear configuration will compete with the DM perpendicular configuration. This results in a non-collinear configuration, as will be shown below.
We note that when the magnetic film has more than one layer, the angle between NN spins in each magnetic layer is different. The determination of the angles is analytically difficult. We have to recourse to the numerical method known as the "steepest descent method" to minimize the energy to get the ground state (GS): We calculate the local field acting on a spin and align it in the direction of the local field. We go to another spin and do the same thing until all spins are visited. We repeat the operation many times until the total energy reaches its minimum.
In the simulations, a sample size of N × N × L was used, with the linear lateral size N = 60 and thickness L = L m + L f , where L m = L f = 4 (L m is the magnetic layer's thickness and L f is the ferroelectric layer's thickness). We use the periodic boundary conditions in the xy plane.
For simplicity, we take exchange parameters between NN spins and NN polarizations equal to 1-namely, J m = J f = 1-for the simulations. We investigate the effects of the interaction parameters (J 2m , J 2 f ) and J m f . We note that the steepest descent method calculates the GS down to the value J m f = −1.25. For values lower than this, the DM interaction is so strong that the spin-spin angle θ tends towards π/2 so that magnetic exchange terms are zero.
We now consider a case with the frustrated regime for (J 2 f , J 2m ) ∈ (−0.4, 0), namely above the critical value of -0.5, as mentioned above.
The spin configuration in the case where H = 0 is shown in Figure 2 for the interface magnetic layer. We observe here a stripe phase with long islands and domain walls. The inside magnetic layers have the same texture.
When H is increased, we observe the skyrmion crystal, as seen in Figure 3: The GS configuration of the interface and beneath the magnetic layers was obtained for J m f = −1.25, with J 2m = J 2 f = −0.2 and external magnetic field H = 0.25. A zoom of the skyrmion shown in Figure 3c and the z-components across the skyrmion shown in Figure 3d indicate that the skyrmion is of a Bloch type. At this field strength of H = 0.25, if we increase the frustration-for example, J 2m = J 2 f = −0.3-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. This is shown in Figure 4, where the interface and the second layer are displayed. The highest value of frustration where the skyrmion structure can be observed is when J 2m = J 2 f = −0.4 close to the critical value -0.5. We show this case in Figure 5: The GS configuration of the interface (a) and second (interior) (b) magnetic layers are presented. Other parameters are the same as in the previous figures, namely J m f = −1.25 and H = 0.25. We can observe a clear 3D skyrmion crystal structure in all magnetic layers, not only near the interface layer. Unlike in the case where we did not take into account the interaction between NNN [47], in the present case, where the frustration is very strong, we see that a large number of skyrmions are distributed over the all magnetic layers with a certain periodicity close to that of a perfect crystal. Though we take the same values for J 2m and J 2 f in the figures shown above, it is obvious that only the magnetic frustration J 2m is important for the skyrmion structure. The ferroelectric frustration affects only the stability of the polarizations at the interface. As long as J 2 f does not exceed -0.5, the skyrmions are not affected by J 2 f . In Figure 6 Figure 5. We conclude here that, when the magnetic frustration is strong enough, the ferroelectric frustration does not affect the skyrmion structure. , where skyrmions are clearly formed, we conclude that while magnetic frustration J 2m enhances the formation of skyrmions, the ferroelectric frustration J 2 f in the weak magnetic frustration tends to suppress skyrmions. The mechanism of these parameters when acting together seems to be very complicated. Let us now show the effect of H. In the case of zero frustration, skyrmions disappear with an external field larger than H = 0.4 [47]. In the case in which we take into account the negative interaction between NNN neighbors, the skyrmion structure is stable with an external field of up to H = 1.0 (see Figure 9). The spins are almost aligned in the direction of the field.  Figure 10b shows the dependence of the ratio of the number of skyrmions on the interior layer N 2 to that on the interface layer N 1 . We see that, as the frustration becomes stronger, the ratio N 2 /N 1 tends towards 1. Let us discuss about some theoretical observations of skyrmions in frustrated magnets [31,[37][38][39][40][41][42][43]46]. Each of these works used a different model, so comparison is impossible. However, all show very similar skyrmion textures. For experiments, many observations have been made in various magnetic materials [2,4,5,7,[19][20][21][22][23][24][25][26][27], in multiferroic materials [28], in ferroelectric materials [29], in semiconductors [30], and in helimagnets [2,20]. Again, as each real material corresponds to a particular microscopic mechanism, the comparison is not simple. However, one can observe many similar topological textures.

Skyrmion Phase Transition
The magnetic transition is driven by the competition between T, the DM interaction (namely < P k > ), the field H, and the magnetic texture (skyrmions). The stronger < P k > and/or J m f is, the higher the transition temperature T c of the skyrmion structure will be. As mentioned above, strong DM interaction helps stabilize the skyrmion crystal [76,77] at the superlattice interface. We used a strong J m f , as in the previous section.
We used the Metropolis algorithm [78,79] to simulate the system at T = 0. We performed calculations for systems with different sizes N × N × L, where N varied from 40 to 100 and the thickness L varied from 2 to 16. It should be noted that changing the lateral size of N does not affect the results on skyrmions shown in the article. However, the influence of the total thickness L of the magnetic and ferroelectric layers is very significant: With an increase in the thickness of the magnetic and ferroelectric layers from L m = L f = 4 to 8, skyrmions were formed only near the interface, not on layers far inside. In most calculations, we used N = 60 and L m = L f = 4. With this thickness, skyrmions were observed in the two interior layers, as seen in the previous section. Usually, we discard 10 6 Monte Carlo steps (MCS) per spin in order to equilibrate the system and average physical quantities over the next 10 6 MCS/spin. Such long Monte Carlo times are needed, since it has been tested for the skyrmion crystal similar to that of the present model [31].
For the ferroelectric layers, the order parameter M f (n) of layer n is given by where .. indicates the thermal average. For the magnetic layers where the spin configuration is not collinear, the definition of an order parameter is not easy. One way to proceed is to heat the system from a selected GS configuration. At a given T, we compare the actual spin configuration observed at the time t with its GS. The comparison is made by projecting that configuration onto the GS. The order parameter of layer n can be thus calculated by where S i (T, t) is the i-th spin at the time t and at T, and S 0 i (T = 0) denotes its orientation at T = 0. We see that M m (n) is close to 1 at very low T, where each spin is close to its orientation in the GS. At high T, where every spin strongly fluctuates, M m (n) becomes zero.
Note that M m (n) is defined in a similar way to the Edwards-Anderson (EA) order parameter in spin glasses (SG) [80]. The EA order parameter is calculated by taking the time average of each spin. When it is frozen at low T, its time average is not zero. At high T, it fluctuates greatly with time so that its time average is zero. The EA order parameter is just the sum of the squares of each spin average. It expresses the degree of freezing, but does not express the kind of ordering.
Note that if the system makes a global rotation during the simulation, then ∑ t a t=t 0 S i (T, t) · S 0 i (T = 0) = 0 for a long time average. To check this, the most efficient way to proceed is to calculate the relaxation time to obtain properties at the infinite time, in the same spirit as the finite-size scaling used to calculate properties at the infinite system size. We have calculated the relaxation time of the 2D skyrmion crystal [31] using the order parameter defined by Equation (9). We have found that skyrmions need more than 10 6 MCS/spin to relax to equilibrium, and the order parameter follows a stretched exponential law as in SG for T < T c .
Another way to check the stability of the skyrmion crystal is to numerically count the topological charges around each skyrmion [81]. If there is a phase transition, the charge number, which plays the role of an order parameter, evolves with T and goes to zero at the phase transition.
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 11, we display the magnetic energy and the magnetic order parameter vs. T 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 11a where Z = 6 (simple cubic lattice), S = 1 (spin magnitude), and k B = 1.3807 × 10 −23 Joules/Kelvin have been used. J e f f is a combination of J m , J 2m , and J m f . Knowing the GS, we can deduce these interactions. We can then calculate the energy in units of Joules by multiplying the value of E m in Figure 11a by the value of J m . Unfortunately, at the time being, there is no experimental energy measurement for comparison. The magnetic order parameters shown in Figure 11b confirm the skyrmion transition temperatures seen by the curvature change of the energy in Figure 11a.  Figure 12b shows the ferroelectric order parameters for the three sets of NNN interactions shown in Figure 12a. These curves confirm the transition temperatures given above.

Conclusions
In this paper, we have studied the effect of the NNN interactions in both magnetic and ferroelectric layers of a magneto-ferroelectric superlattice. A Dzyaloshinskii-Moriya (DM) interaction was assumed for the magneto-ferroelectric interface coupling.
We found the formation of a skyrmion crystal in the GS under an applied magnetic field in a large region of parameters in the space (J 2m , J m f ). As expected, the magnetic frustration enhances the creation of skyrmions while, when strong enough, the ferroelectric frustration destabilizes skyrmions if the magnetic frustration is weak.
We have studied the phase transition of the skyrmion crystal through the use of Monte Carlo method. Skyrmions have been shown to be stable at finite temperatures. While the magnetic frustration helps enhance the creation of skyrmions, it reduces the transition temperature considerably.
The existence of very stable skyrmions confined at the magneto-ferroelectric interface at finite T is very interesting and may have potential applications in spintronics. Many applications using skyrmions have been mentioned in the Introduction. As a last remark, let us mention that the present magneto-ferroelectric superlattice model can be used in the cases of magnetic monolayers or bilayers to study the dynamics of the skyrmions driven by a spin-polarized current or by a spin-transfer torque. Due to the small thickness, the skyrmions created by the interface are as well confined as in 2D. Our model is therefore suitable for creating skyrmion pinning by using an electric field acting on the ferroelectric polarizations. This is the subject of our future investigations.