Signatures of a liquid-crystal transition in spin-wave excitations of skyrmions

Understanding the spin-wave excitations of chiral magnetic order, such as the skyrmion crystal (SkX), is of fundamental interest to confirm such exotic magnetic order. The SkX is realized by competing Dzyaloshinskii-Moriya and ferromagnetic-exchange interactions with a magnetic field or anisotropy. Here, we compute the dynamical spin structure factor, using Monte Carlo and spin dynamics simulations, extracting the spin-wave spectrum in the SkX, in the vicinity of the paramagnet to SkX transition. Inside the SkX, we find six spin-wave modes, which are supplemented by another mode originating from the ferromagnetic background. Above the critical temperature Ts for the skyrmion crystallization, we find a diffusive regime, reminiscent of the liquid-to-crystal transition, revealing that topological spin texture of skyrmionic character starts to develop above Ts as the precursor of the SkX. We discuss the opportunities for the detection of the spin waves of the SkX using inelastic-neutron-scattering experiments in manganite-iridate heterostructures. https://doi.org/10.1038/s42005-020-00489-w OPEN

The SkX phase exhibits spin-wave excitations, which are Goldstone modes associated with the spontaneously-broken translational symmetry of the chiral spin arrangement. The dispersion of these spin-wave modes is determined by the properties of the SkX. Three SkX spin-wave modes-viz. clockwise, counterclockwise circulation and breathing modes-have been predicted in previous theoretical analysis [32][33][34] and subsequently reported in experiments using broadband microwave spectroscopy [35][36][37] . Some of the spin-wave modes exhibit non-trivial topological properties, revealing in the form of chiral edge states in the SkX 38 . The transition to the SkX phase takes place within a range of external magnetic field strengths, when the temperature is slowly lowered down 39 .
An important puzzle which remained unresolved for years is the role of spin fluctuations in the transition to the SkX phase. This is particularly interesting since specific heat measurements report a first-order phase transition 40 to the helimagnetic phase at zero magnetic field, while neutron-scattering experiments reveal a Landau soft-mode mechanism of weak crystallization to the SkX phase at a finite field [41][42][43] . These experiments suggest that the emergence of the SkX phase occurs via a precursor regime in which non-trivial topological spin textures of skyrmionic character are developed within a paramagnetic background with abundance of fluctuations, similar to the precursor phenomenon proposed in the context of liquid-crystal transitions 44 . At low finite temperatures, the positional correlations of the skyrmions decay with distance as power laws while the orientational correlations remain finite 45 . With increasing temperature, a transition takes place from the solid to a liquid phase of skyrmions with asymptotically vanishing long-range correlations.
Here, we theoretically investigate the spin-wave excitations of a two-dimensional SkX phase and its evolution at higher temperatures, in the vicinity of the paramagnetic-to-SkX phase transition, by computing the dynamical spin structure factor S (q, ω), expecting its future direct detection in inelastic SANS experiments. We use a spin Hamiltonian to first find out the stable spin configurations at low temperatures using the Metropolis Monte Carlo annealing method, and then numerically solve the Landau-Lifshitz equation of motion to capture the spin dynamics. The computed S(q, ω) reveals an exotic array of six gapless spin-wave modes in the SkX phase, stabilized within a range of the external perpendicular magnetic fields, accompanied by another gapless mode originating from the ferromagnetic background. In the diffusive regime for the skyrmion crystallization, S(q, ω = 0) profile reveals a clear transition from a ringshaped profile to a six-peak profile, characteristic of the SkX phase, confirming the previous experimental findings that the transition to the SkX phase upon cooling slowly from a higher temperature occurs via a precursor regime. A gradual decrease in the S(q = 0, ω) intensity with temperature, at the SkX spin-wave modes, also indicates a phase transition similar to the liquidcrystal transition.
The spin spiral (SS) phase, which is the natural solution at zero magnetic field, exhibits two gapless spin-wave modes from the spiral order and a soft ferromagnetic mode that appears at finite energies and goes away with increasing the DMI strength. The field-polarized ferromagnetic (FM) phase exhibits a spin-wave mode, which appears with an excitation gap due to the Larmor precession in a magnetic field.
A particularly interesting platform that hosts a twodimensional SkX is oxide interfaces, such as the interface between the ferromagnetic metal La 1−x Sr x MnO 3 (0.2 ⪅ x ⪅ 0.5) and the non-magnetic semimetal SrIrO 3 . The Dzyaloshinskii-Moriya interaction (DMI), which arises at this interface due to the strong spin-orbit coupling in SrIrO 3 and the broken inversion symmetry, stabilizes a Néel-type SkX 39 . We discuss the feasibility of the experimental detection of our results in this manganite/ iridate interfaces based on the state-of-the-art SANS spectrometers and the developments required for a successful future confirmation of our predictions.

Results
Model Hamiltonian. We consider a model Hamiltonian, relevant to two-dimensional material systems such as the La 1−x Sr x MnO 3 / SrIrO 3 , Pt/Co, and Fe/Ir interfaces. We consider the manganite/ iridate interface as our example system and choose parameters that are relevant to this interface. The manganite is considered to be doped suitably to be in the ferromagnetic phase in which the localized spins can be described effectively by a Heisenberg exchange interaction. The influence of the strong spin-orbit coupling of the iridate, in the absence of explicit structural inversion symmetry, leads to the DMI at the interface. The total Hamiltonian, in the presence of an external magnetic field, is given by where J is the nearest-neighbor ferromagnetic exchange parameter, D is the DMI strength, H z is the magnetic field applied perpendicular to the interfacial plane, and A is the easy-plane magnetic anisotropy originating from the combined interfacial strain and Rashba spin-orbit coupling 46 . We used A = 0.01 J throughout this paper. S i is the localized spin of amplitude S = 3/ 2 on the Mn t 2g orbitals at site i on the manganite side of the interface. By comparing the critical temperature for the zero-field FM phase (at D = 0) extracted from our Hamiltonian with the experimental findings 47 , it was found that the Heisenberg exchange energy constant at the manganite/iridate films varies within the range 1.2 meV ⪅ J ⪅ 3.9 meV, depending upon the thickness of the film 39 . In the description below, we present energies in units of J and reveal magnetic field and temperature scales, at selected places, using J = 1.2 meV. The dipole-dipole interaction usually has an important contribution in ferromagnetic materials. However, recent experiments suggest that the short-range exchange interaction dominates over the long-range dipole-dipole interaction in Sr-doped manganite compounds 48 . While the long-range dipole-dipole interaction may be sizeable and help in stabilizing the magnetic phases in bulk manganite systems, here we focus on the short-range exchange interactions which are more relevant to the two-dimensional geometry.
Low-temperature spin configurations. We present the spin-wave excitations, obtained using a 100 × 100 lattice at different magnetic fields, first at a low temperature T = 0.001 J and then discuss the high-temperature diffusive regime in which the liquid-crystal transition takes place. The SS phase, which is stabilized at zero magnetic field, evolves into a triangular SkX within a range of magnetic fields and finally becomes a fully-polarized FM phase at higher fields. A detailed analysis of the phase diagram, spanned by the magnetic field and the temperature, was done in our previous work 39 . In Fig. 1a, b, we show the spin configurations of the SS phase and the SkX phase, obtained in a typical annealing session of the Monte Carlo simulations, respectively at zero field (H z = 0) and at a finite field (H z = 0.18 J). Since the DM vectors lie predominantly at the interface plane, the SS and the SkX textures are Néel-type in nature, unlike in three-dimensional chiral magnets such as MnSi in which out-of-plane DM vectors lead to Blochtype textures. A change in the sign of D reverses the sense of rotation of the spins in going from the skyrmion center towards the radially outward direction. At the considered DMI strength D = 0.5 J, which is close to the value realizable at oxide interfaces [49][50][51] , the period of the SS and the skyrmion diameter are nearly 12 lattice spacings. At smaller values of D, these length scales increase almost exponentially, D = 0 (H z = 0) being the rotationally-invariant Heisenberg FM phase.
Spin waves in the skyrmion crystal. The SkX phase is formed within an intermediate magnetic-field range in which, upon cooling, the skyrmions establish a long-range order to form a triangular crystal. At the DMI strength D = 0.5 J, the SkX phase is realized within the field range 0.13 J ⪅ H z ⪅ 0.20 J. We focus on the SkX phase obtained at H z = 0.18 J, shown in Fig. 1b. In Fig. 2a-f, we show the constant-energy contours of the dynamical spin structure factor S(q, ω) at different energies. The plot at ω = 0 in Fig. 2a depicts the six-peak structure of the SkX phase, as revealed in elastic neutron-scattering experiments (see e.g., Fig. 2e in ref. 9 ). The three characteristic momenta of the triangular SkX phase are denoted by Q 1 , Q 2 , and Q 3 . The peak at the Γ point originates from the zero-momentum magnetic ordering in the ferromagnetic background of the SkX. The constant-energy contours of S(q, ω) at different finite energies are shown in Fig. 2b-f. At a finite energy ω = J, the seven elastic peaks evolve into circles as shown in Fig. 2b. With the increase in energy, these circles grow in size and deviate their shape to rounded squares near energy ω = 6 J ( Fig. 2d). At much higher energies, the topology of the energy contour changes (characteristic of the Lifshitz transition), seven pockets appear at the corners of the Brillouin zone ( Fig. 2e) and these modes disappear at energies above ω ≈ 13 J.
In Fig. 2g, we present the momentum-energy dependence of S (q, ω) of the SkX phase, the main result of the paper, along the high-symmetry path The plot unveils a complex structure for the spin-wave modes, six of which originate from the SkX and another from the ferromagnetic background. The spin-wave modes from the SkX appear to follow a parabolic dispersion relation (ω ∝ q 2 ) in the low-energy regime ω ⪅ 0.2 J, around each elastic peak. As shown in Fig. 2h, the spinwave modes of the SkX along the Γ − X direction are twofold degenerate. The rich features shown in Fig. 2g are a remarkable result and its possible experimental confirmation is addressed below.
In the presence of an additional AC magnetic field, three spinwave modes of the SkX phase were numerically identified in the ref. 32 . Two of these three modes are rotational modes (clockwise and counterclockwise) that appear with in-plane AC magnetic field. The third one which was found with out-of-plane AC magnetic field is called the "breathing" mode, where the skyrmion core expands and shrinks alternatively. These modes were subsequently reported in experiments using the broadband microwave spectroscopy in the skyrmion-host compounds MnSi, Fe 1−x Co x Si, Cu 2 OSeO 3 , and GaV 4 S 8 [35][36][37] . We speculate that these three modes correspond to the excitations at energy ω ≈ 0.5 J at the Γ point in Fig. 2g. Previous efforts 33,34,52 analytically obtained the spin-wave modes of the SkX that we explore here. Below, we discuss the feasibility of the detection of the spin-wave modes, obtained in our calculations at constant magnetic field, in inelastic neutron-scattering experiments.
Diffusive regime for skyrmion crystallization. We now focus on the high-temperature regime in which the crystallization of the skyrmions takes place as the temperature is lowered down. In our model, the triangular SkX phase is stabilized below a temperature T s ⪅ 0.8 J (which, in Kelvin unit, lies in the range 11 K ⪅ T s ⪅ 35 K, depending on the thickness of the manganite layer), as revealed by the temperature dependence of the skyrmion number ∂y Þdxdy, plotted in Fig. 3a. The average magnetization component 〈M z 〉 also becomes nearly constant below T s . However, the development of the skyrmionic correlation begins at temperatures below T ≈ 2.5 J, much above T s , in the annealing process when the temperature is reduced gradually. The temperature range 0.8J ⪅ T ⪅ 2.5 J is the weak crystallization regime at the considered values of D and H z . The snapshots of the spin configuration at two temperatures T = J and T = 2 J are shown, respectively in Fig. 3b, d. This diffusive regime may also involve an extended phase of nucleated skyrmions, which we explored in our previous study 39 . The constant-energy contour of the dynamical spin structure factor S(q, ω) at ω = 0 and at temperature T = 2 J, in Fig. 3e, shows a ring-shaped profile. At the smaller temperature T = J, in Fig. 3c, the six-peak structure, which is a reminiscent of the triangular crystal of skyrmions, appears with the ring in the background. Such a paramagnetic-toskyrmion lattice transition has been observed in the transition metal helimagnet MnSi using small-angle neutron scattering, neutron-resonance spin-echo spectroscopy, and microwave spectroscopy experiments 43 . Figure 3f, g show the energy-momentum dependence of S(q, ω) along the momentum path Γ − Q 1 − Q 2 − Q 3 − Γ at the two temperatures T = J and T = 2 J, respectively. At T = J, the spinwave modes of the SkX can be resolved well from the merged background of the diffusively-scattered intensities. On the contrary, at T = 2 J, the ring-shaped elastic profile do not reveal any clear mode of the spin-wave branches. Interestingly, however, the modes near ω = 0.51 J at q = 0 persist, despite a smearing in the intensities within a narrow energy range. In Fig. 4, we show the evolution of this ω = 0.51 J mode at different temperatures across the paramagnet-SkX transition and find that the intensity of S(q, ω) decreases linearly below the temperature T s ≤ 3J. The absence of a sharp transition in the intensity of S(q, ω) further indicates the occurrence of a liquid to crystal transition leading to the SkX phase. The measurement of these smeared modes of S (q, ω) at zero wave vector can be useful to distinguish the elastic neutron-scattering pattern of a SkX from the nearly-similar one of the hexagonal ferrofluids 53 , iron oxide nanoparticles 54 , hexagonal magnets 55 , and possibly ferromagnetic domains.
The transition from the ring shape to the six-peak structure of elastic S(q, ω) profile upon reducing temperature in the diffusive regime establishes the fact, that fluctuating skyrmionic correlation starts to develop at temperatures much above the SkX-ordering temperature T s as the precursor of the long-range triangular crystal of skyrmions. Such a precursor phenomenon has been reported in the cubic helimagnets FeGe and MnSi using ac-magnetic susceptibility and magneto-heat capacity measurements 41,42 .
Spin waves in the spin spiral. We now turn our attention to the SS phase that appears spontaneously in the absence of any external magnetic field at low temperatures in our model. The SS phase has a single characteristic momentum Q 1 that is governed by the period of the spiral texture. It is shown in the constantenergy contour of S(q, ω) at ω = 0 in Fig. 5a. With an increase in energy, the two elastic peaks at Q 1 and Q 2 = −Q 1 , similarly as in the SkX phase, evolve into circles around these characteristic momenta, as shown by the constant-energy contour of S(q, ω) at ω = J in Fig. 5b. In addition to the two spin-wave modes of the spiral texture, a soft ferromagnetic mode appears at finite energies, as shown by the circle around the Γ point in Fig. 5b. The intensity of this ferromagnetic mode decreases with a decrease in the period of the SS texture, i.e., with an increase in the DMI strength D. The constant-energy contours of S(q, ω) at different energies in Fig. 5c-f show the spin-wave modes behaving in a similar manner as in the FM or in the SkX phase.
We show the momentum-energy dependence of S(q, ω) of the SS phase along the momentum path Q 1 − Γ − Q 2 in Fig. 5g, and along Γ − X − M − Γ in Fig. 5h. The plots show two spin-wave modes from the spiral order and one mode from the ferromagnetic order with a bandwidth~13 J. The spin-wave modes from the SS phase follow nearly a parabolic dispersion relation (ω ∝ q 2 ) in the low-energy regime ω ⪅ J. The two-fold degenerate spin-wave excitation at ω ≈ 0.8 J at the Γ point (Fig. 5g) is detectable in microwave spectroscopy experiments.
Theoretical calculations 34 and inelastic neutron-scattering experiments in the context of chiral magnets 35,56-64 analyzed the "helimagnon" excitation modes of the SS phase.  A ferromagnetic spin-wave mode, corresponding to the zero wave-vector magnetic order, is expected to robustly emerge from ω = 0 in the case of the conical spin spiral phase, that appears in three-dimensional chiral magnets at zero magnetic field and in two-dimensional oxide interfaces with an in-plane magnetic field 34 .
Spin waves in the field-polarized ferromagnet. Having explored all other magnetic phases in our model, we now discuss, for concreteness, the spin-wave mode of the field-polarized FM phase. We realized this FM phase at a magnetic field H z = 0.5 J in the MC simulations. Figure 6a shows the dynamical spin structure factor S(q, ω) computed along a momentum path Γ(0, 0) − X (π, 0) − M(π, π) − Γ in the Brillouin zone.
The plot describes the well-known spin-wave dispersion relation ω ¼ H z þ 4JSðcos q x þ cos q y Þ with a bandwidth~12 J and an energy gap Δω = 0.5 J. Figure 6b-e show the constant-energy contours of the spin-wave dispersion at different energies. The point-like peak in S(q, ω) at the Γ point at energy ω = 0.5 J evolves into a circle with increasing ω and the circle takes the shape of a rounded square at ω ≈ 6 J. With further increase in ω and four pockets appear at the corners of the Brillouin zone, like in the cases of the SkX and the SS phases. These hole pockets disappear above the upper limit, ω ≈ 12.5 J, of the energy spectrum.
The sharpness of the S(q, ω) dispersion in the frequencymomentum space is controlled primarily by the lattice size and the time step. At smaller lattice sizes, the S(q, ω) peaks appear discontinuously along the dispersion curve. With larger lattice sizes, the increase in the number of peaks form a continuous sharply-defined curve. An analysis of finite-size effects on the spin-wave dispersion was done in ref. 65 . On the other hand, a larger time step introduces tails along the ω axis surrounding the peak positions with an amplitude, that decays in an oscillatory manner with distance from the ω of the peak. A sufficiently small time step reduces these tails significantly from the dispersion curve. In three-dimensional materials with a finite DMI, the fullypolarized FM phase at a finite magnetic field exhibits an anisotropy in the spin-wave dispersion, the anisotropy being proportional to the DMI strength D and the perpendicular momentum k z 34,57 . In two dimensions, the anisotropy goes away and the spin-wave dispersion is similar to a two-dimensional ferromagnet.

Discussion
The presented dynamical spin structure factor of the skyrmion crystal reveals a rich structure of the spin-wave excitations that can be achieved at the maganite/iridate or similar interfaces where a skyrmion crystal is expected to be formed. Its detection in inelastic neutron-scattering experiments can offer an opportunity to effectively probe the topological magnetic textures, which is often difficult to deduce conclusively from only the topological Hall effect measurements. The challenge for the successful detection of the SkX spin-wave modes in inelastic scattering experiments lies in obtaining a sufficiently high signalto-noise ratio from a single interface. For this purpose, multilayer heterostructures, composed of alternate manganite and iridate layers can provide a workable platform. We note that with the mass of a typical La 1−x Sr x MnO 3 /SrIrO 3 /SrTiO 3 superlattice (of thickness 1-10 μm deposited on a 1 × 1 cm 2 substrate) is 0.00064 g (with average density~6.4 g/cm 3 ), implying that a sample mass of 0.1 g can be achieved by stacking~15-150 of these superlattices. Considering that LaMnO 3 has a molar mass 239.3 g/mol, a 0.1 g sample contains 0.4 millimoles, while a generic criterion for inelastic neutron scattering experiments is to have samples of about 6.2 (3.5) millimoles for S = 3/2 (S = 2) Mn spins. Thus, measurements with these small samples are likely a factor 10 beyond current capabilities, such as the Spallation Neutron Source in the United States, the Institut Laue-Langevin in France, the ISIS neutron source in the United Kingdom, and J-PARC in Japan. There are additional considerations such as neutron absorption, background from the substrate scattering, bandwidth of the excitations, achieving appropriate coverage in reciprocal space, etc for a successful neutron-scattering experiment. However, such experiments are likely feasible at planned or under construction next generation neutron sources, such as the European Spallation Source or the Second Target Station at the Spallation Neutron Source, Oak Ridge National Laboratory, where improvements in instrumentation and neutron flux are anticipated to result in performance gains in excess of an order of magnitude [66][67][68] .
Although the challenges related to the thin-film geometry of oxide interfaces are considerable with regards to sample mass, if robust results could be experimentally gathered they will enrich the understanding of the dynamics of the skyrmions, which are of high technological interests owing to their potential for wide-spread applications in the high-density data storage 69 and magnonic devices 70 . Another merit of the presented results and the proposed inelastic measurements would be to distinguish the elastic neutron-scattering pattern of a skyrmion crystal from the nearly-similar one of the hexagonal ferrofluids 53 , iron oxide nanoparticles 54 , and hexagonal magnets 55 . To conclude, we theoretically investigated the spin-wave excitations of a two-dimensional Dzyaloshinskii-Moriya magnet, in the context of the two-dimensional interface between lanthanum manganite and strontium iridate, by computing the dynamical spin structure factor. We explore the excitation modes appearing in the SkX, SS, and FM phases obtained using Monte Carlo simulations and Landau-Lifshitz spin dynamics at different magnetic fields. In particular, we focus on the SkX phase which shows a complex interesting structure in the spin dynamical spin structure factor at low temperatures. We also present the spinwave excitations in the diffusive regime above the SkX-ordering temperature T s and show that fluctuating skyrmionic correlation develops at temperatures much above T s as the precursor of the long-range SkX phase, indicative of a liquid-crystal transition. We finally discussed the possible challenges in the experimental detection and propose a multilayer heterostructure geometry that might overcome some of the obstacles to successfully observe these spin-wave modes from the SkX phase at oxide interfaces.

Methods
Monte Carlo annealing. We obtain the spin configurations on a square lattice of size 100 × 100 with periodic boundary conditions using the Monte Carlo (MC) annealing procedure. For all considered magnetic fields, the annealing process was started at a high temperature T = 30J with a completely random spin configuration and the temperature was lowered slowly down to a low value T = 0.001J in 1000 temperature steps. At each temperature, 10 9 MC spin updates were performed. In each spin update step, a new spin direction was chosen randomly within a small cone spanned around the initial spin direction. The MC sampling method provides a full and uniform coverage of the phase space spanned by the spin angles, as illustrated in Fig. 7. The new spin configuration was accepted or rejected according to the standard Metropolis algorithm by comparing the total energies of the previous and the new trial spin configurations, calculated using the Hamiltonian in Eq. (1).
Landau-Lifshitz spin dynamics. Starting with the initial conditions obtained from MC simulations at a given temperature, the spins are integrated numerically and averaged over the initial values using the Landau-Lifshitz equations of motion where H is the Hamiltonian described in Eq. (1) and S r (t) is the spin vector at position r and time t.
The central object of investigation is the dynamical spin structure factor S(q, ω), expressed as Sðq; ωÞ ¼ where α and β represent the x, y, and z components. N is the total number of lattice sites and the parameter η = 0.01 was used to obtain S(q, ω) at q = 0. The factor ðδ αβ À q α q β q 2 þη 2 Þ was used in Eq. (3) to include the off-diagonal components of the correlation function C αβ ðr À r 0 ; tÞ (α ≠ β).
While solving Eq. (2), we express the spins using spherical polar coordinates and evaluate the evolution of the polar and the azimuthal angles θ and ϕ. The time integration in Eq. (2) was carried out using the fourth-order Runge-Kutta method to a maximum time t max = 1000J −1 with a time step Δt = 0.001J −1 . The finite time cutoff may introduce small oscillations in S(q, ω) along the ω axis, resulting from the time Fourier transform. These oscillations were minimized using the convolution of the spin-spin correlation function with a resolution function in frequency, equivalent to the energy resolution of the SANS spectrometer. A Gaussian broadening function with a broadening parameter σ ω yields the modulated dynamical spin structure factor 65 S 0 ðq; ω 0 Þexp À ðω À ω 0 Þ 2 2σ 2 where S 0 ðq; ω 0 Þ is the dynamical spin structure factor, calculated using a time cutoff. In the results presented in this paper, σ ω = 0.001J was sufficient to resolve different spin-wave modes.
The Gilbert damping may also introduce a broadening of the spin-wave bands. However, the Gilbert damping coefficient (α) for manganite systems is very low (~10 −4 ) [see, for example, ref. 71 ]. These small values of α will introduce only a tiny amount of broadening to the higher-energy spin-wave excitations, and therefore, it can be safely ignored in the context of manganite thin films.
The computed S(q, ω) was averaged using 50 independent thermalized realizations of the MC spin configurations to reduce the noise arising from small statistical imperfections in those configurations. Increasing this number did not produce noticeable modifications in the results. One advantage of the current approach over linear spin-wave theory techniques is the scope to examine spin dynamics at finite temperatures, where the spin system does not have a wellestablished long-range order.

Data availability
The datasets obtained during this work are available from the corresponding author upon reasonable request.