Phase transitions in chiral magnets from Monte Carlo simulations

TL;DR

This study uses Monte Carlo simulations to explore phase transitions in a model of chiral magnets, revealing how competing interactions cause a crossover from second-order to first-order transitions, explaining experimental observations in MnSi.

Contribution

It demonstrates how the ratio of Dzyaloshinskii-Moriya to exchange interactions influences the order of phase transitions in chiral magnets, linking microscopic interactions to macroscopic phase behavior.

Findings

The specific heat hump in MnSi is due to smearing of a second-order transition by helical fluctuations.

Increasing D/J causes a crossover from second-order to first-order phase transition.

Competing interactions are key to the first-order transitions in helical magnets.

Abstract

Motivated by the unusual temperature dependence of the specific heat in MnSi, comprising a combination of a sharp first-order feature accompanied by a broad hump, we study the extended Heisenberg model with competing exchange $J$ and anisotropic Dzyaloshinskii-Moriya $D$ interactions in a broad range of ratio $D/J$. Utilizing classical Monte Carlo simulations we find an evolution of the temperature dependence of the specific heat and magnetic susceptibility with variation of $D/J$. Combined with an analysis of the Bragg intensity patterns, we clearly demonstrate that the observed puzzling hump in the specific heat of MnSi originates from smearing out of the virtual ferromagnetic second order phase transition by helical fluctuations, which manifest themselves in the transient multiple spiral state. These fluctuations finally condense into the helical ordered phase via a first order phase transition as is indicated by the specific heat peak. Thus the model demonstrates a crossover from a second-order to a first-order transition with increasing $D/J$. Upon further increasing $D/J$ another crossover from a first-order to a second-order transition takes place in the system. Moreover, the results of the calculations clearly indicate that these competing interactions are the primary factor responsible for the appearance of first order phase transitions in helical magnets with the Dzyaloshinskii-Moriya (DM) interaction.

Equations

$$H_J=-J\sum _i{\mathbf{S}}_i\cdot ({\mathbf{S}}_{i+\hat{x}}+{\mathbf{S}}_{i+\hat{y}}+{\mathbf{S}}_{i+\hat{z}})-J^{\prime }\sum _i{\mathbf{S}}_i\cdot ({\mathbf{S}}_{i+2\hat{x}}+{\mathbf{S}}_{i+2\hat{y}}+{\mathbf{S}}_{i+2\hat{z}}),$$

$$H_D=-D\sum _i\left({\mathbf{S}}_i\times {\mathbf{S}}_{i+\hat{x}}\cdot \hat{x}+{\mathbf{S}}_i\times {\mathbf{S}}_{i+\hat{y}}\cdot \hat{y}+{\mathbf{S}}_i\times {\mathbf{S}}_{i+\hat{z}}\cdot \hat{z}\right)-D^{\prime }\sum _i\left({\mathbf{S}}_i\times {\mathbf{S}}_{i+2\hat{x}}\cdot \hat{x}+{\mathbf{S}}_i\times {\mathbf{S}}_{i+2\hat{y}}\cdot \hat{y}+{\mathbf{S}}_i\times {\mathbf{S}}_{i+2\hat{z}}\cdot \hat{z}\right).$$

Full text

Introduction

The experimental and theoretical study of the model itinerant helimagnet MnSi is a dynamic area of active research which still exhibits many unsolved problems [1, 2]. In particular, a tiny first order phase transition in MnSi and its accompanying phenomena present a real challenge that prevent a full understanding of the material. In this Letter we perform classical Monte Carlo simulations to get new insights on phase transitions in chiral spin systems. As is known the crystal structure of MnSi belongs to the cubic space group P$2_{1}3$ , which does not contain a center of inversion. This results in a non-zero value of the antisymmetric DM interaction term $D{\bf S}_{i}\times{\bf S}_{j}$. While the Heisenberg ferromagnetic interaction $J{\bf S}_{i}\cdot{\bf S}_{j}$ favors homogeneously aligned configurations of spins, the DM interaction leads to a long-wave modulation of the magnetic spin structure in a spiral type, appearing below a certain temperature $T_{c}$ ($T_{c}\approx 29$ K in MnSi) with a pitch wave vector ${\bf Q}$ directed along the $[111]$ crystallographic direction [2, 3]. The magnetic moments are ferromagnetically aligned in planes perpendicular to $[111]$, and the planes are consecutively turned by some angle $\theta$ relative to each other, with $\tan\theta=D/J$.

There are strong indications that the magnetic transition in helimagnets is of first order [3, 4, 5, 6, 7]. However, the nature of the phase transition in MnSi is still not well understood. This transition shows some remarkable features seen in the specific heat, thermal expansion coefficient, temperature coefficient of the resistivity and sound absorption [3, 4, 5, 8], As can be seen in Fig. 1 these quantities display sharp peaks at the phase transition temperature $T_{c}$ followed by well-defined rounded humps on the high-temperature sides of the curves.

The origin of these enigmatic second maxima or humps is also poorly understood. Analyses of small angle neutron scattering data strongly advanced our knowledge, but still did not deliver completely satisfactory explanations [7, 9, 10]. What seems certain by now is that the domain of the secondary maxima is characterised by strong helical fluctuations evidenced by diffuse neutron scattering, with the intensity distributed on a sphere of radius ${\bf q}={\bf Q}$ in momentum space [7, 9].

Inspired by the successful simulation of a phase diagram of MnSi [11], we address the question on the nature of the phase transition in a chiral helimagnet by using a classical Monte Carlo technique. The chiral helimagnet is modelled by a lattice spin Hamiltonian consisting of Heisenberg exchange and varying DM interaction terms extended to include interactions between next-nearest neighbours. This simple model leads to unusual magnetic properties based on the competing nature of the exchange and DM interactions. In particular, we show that the hump in the physical properties at the phase transition in MnSi originates in the perturbation of a virtual ferromagnetic second order phase transition by helical fluctuations that arise due to the DM interaction. In the other words the ”hump” may be interpreted as a ferromagnetic second order phase transition smeared out by the helical fluctuations, which eventually condense into the helically ordered phase. This conclusion agrees completely with the analysis of the experimental thermodynamic data on the magnetic phase transition in MnSi, performed in [12].

Model and simulation

For a description of the itinerant helimagnet MnSi we adopt a local-moment model [10, 11, 13, 14], whose Hamiltonian $H=H_{J}+H_{D}$ consists of two terms, describing the exchange (J) and DM (D) interactions between nearest-neighbor spins, with DM vector directed along the corresponding bond direction,

[TABLE]

We consider classical Heisenberg spins ${\bf S}_{i}=(S_{i}^{x},S_{i}^{y},S_{i}^{z})$ of unit length, $|{\bf S}_{i}|=1$ placed on a simple cubic lattice $L\times L\times L$ with periodic boundary conditions. The lattice spacing is taken to be unity. The summation is over the sites of the cubic lattice spanned by the vectors $\hat{x}$, $\hat{y}$ and $\hat{z}$. We also supplement exchange and DM interactions with next-nearest neighbor interaction terms $J^{\prime}$ and $D^{\prime}$ to compensate induced anisotropies originating from the discretization of the corresponding continuum spin model, as was first proposed by Buhrandt and Fritz [11]. Setting $J^{\prime}=-J/16$ and $D^{\prime}=-D/8$ allows to compensate these anisotropies to leading order [11]. In all calculations, except the one presented in Fig. 2 below, we held the parameter $J$ fixed at $J=1$, serving as a unit of temperature.

To find thermodynamic properties we perform Monte Carlo (MC) simulations using a standard single-site Metropolis algorithm. The simulation starts at some temperature well above $T_{c}$ from some ordered state. During the equilibration stage of simulation it melts into paramagnetic disordered state. We then gradually decrease the temperature by small steps $\Delta T=10^{-2}$. This enables us to obtain the dependence of specific heat as a function of temperature for varying DM strengths. Calculations were carried out with L= 30 and $5\times 10^{5}$ MC steps per spin at each temperature to acquire the statistics. The averaged lattice spin configuration $\langle{\bf S}_{i}\rangle$ at each $T$ is used to find Fourier components, $\langle{\bf S}_{q}\rangle=1/N\sum_{i}\langle{\bf S}_{i}\rangle e^{-i{\bf q}\cdot{\bf R}_{i}}$, where $N=L^{3}$, and, subsequently, the Bragg intensity profile $I({\bf q})\propto|\langle{\bf S}_{q}\rangle|^{2}$. For ease of presentation, we also introduce the Bragg intensity profiles projected onto the ($q_{x},q_{y}$) plane, $I^{*}(q_{x},q_{y})=\sum_{q_{z}}I(q_{x},q_{y},q_{z})$. Besides the Bragg intensity we have also studied the specific heat $C(T)$ and magnetic susceptibility $\chi(T)=N(\langle{\bf M}^{2}\rangle-\langle{\bf M}\rangle^{2})/T$, where ${\bf M}=(1/N)\sum_{i}{\bf S}_{i}$. We found it most convenient to calculate the specific heat from direct differentiation of the energy $E(T)=\langle H\rangle$. We set $2\pi/L$ as a unit length in the reciprocal space and vectors ${\bf q}$ discussed below are scaled accordingly.

Results and discussion

The Hamiltonian $H$ contains two main terms, the Heisenberg interaction term $H_{J}$, and the antisymmetric DM intraction $H_{D}$. Separately, each of them results in a phase transition of second order at some $T_{c}$ from an ordered phase to a disordered paramagnetic phase. The behavior of the specific heat and magnetic susceptibility for these two transitions is illustrated in Fig. 2 (in both cases the compensation terms $J^{\prime}$ and $D^{\prime}$ are included). The corresponding susceptibility for the Hamiltonian $H_{D}$ (curve 1 in the inset of Fig. 2) shows an antiferromagnetic signature, whereas Hamiltonian $H_{J}$ leads to a singularity in the susceptibility at $T=T_{c}$ (curve 2 in the inset of Fig. 2) as it expected for ferromagnets.

The specific heat shown in Fig. 2 (curve 1 and 2) displays the singular behavior at $T_{c}$ typical of second order transitions. In case of the Heisenberg interaction the ordered ferromagnetic (FM) state is characterized by a single peak in Bragg intensity at zero wave vector. With increasing temperature this peak gradually vanishes signalizing the disappearance of the ordered FM state, at a transition temperature $T_{c}/J\simeq 1.44$ [15]. In our simulation the transition occurs at slightly lower temperature, $T_{c}\simeq 1.28$ due to finite size effects.

In the case of DM interactions the lower-temperature phase is characterized by a spiral magnetic order ${\bf S}({\bf R}_{i})=S_{Q}\cos({\bf Q}\cdot{\bf R}_{i}){\bf e}_{1}+S_{Q}\sin({\bf Q}\cdot{\bf R}_{i}){\bf e}_{2}$, with unit vectors ${\bf e}_{1}$ and ${\bf e}_{2}$ being perpendicular to the spiral vector ${\bf Q}$. Such a single spiral results in two Bragg peaks at ${\bf Q}$ and $-{\bf Q}$. In our simulation the corresponding spiral wave vector turns out to be ${\bf Q}=(2\pi/L)(9,9,9)$, which corresponds to relative rotations of spins in neighboring planes perpendicular to $[111]$ by an angle approximately equal to $\pi/2$ (generally, the magnitude of the spiral vector is affected by the periodic boundary conditions and the competing term $D^{\prime}$). The eventual transition to the paramagnetic state, accompanied by simultaneous diminishing of the Bragg peaks’ intensity, occurs at $T_{c}/D\simeq 0.56$.

For the combined Hamiltonian $H=H_{J}+H_{D}$ the competing nature of the terms $H_{J}$ and $H_{D}$, which have different ground states, results in strong fluctuations that drive the system from a second-order to a first order transition at some values of $D/J$. The hump or shoulder in the specific heat at the high temperature side of the phase transition appears to result from a degradation of the ferromagnetic second order phase transition caused by the intense helical fluctuations. To reach this conclusion we performed a number of runs, varying the parameter $D$, while $J=1$ being fixed. The evolution of the specific heat and susceptibility for small and large values of $D$ are presented in panels (a) and (b) of Fig. 3. At $D<0.3$, the specific heat exhibits a clear maximum, resembling the singular behavior of $C(T)$ at $T_{c}$ for $D=0$. With increasing $D$ this maximum evolves into a round hump, which eventually disappears for $D>1$. Concurrently, a first order sharp peak of $C(T)$ develops at the low temperature side of the hump. Upon further increasing $D$ beyond $D=1$ the phase transition temperature $T_{c}$ also increases while simultaneous changing the shape of the $C(T)$-curve. The transition becomes second order again. Examining the corresponding dependencies of the susceptibility $\chi(T)$ reveals that with increasing $D$ the maximum of $\chi(T)$ evolves into a step-like feature with considerably diminished magnitude of $\chi$. At $D\geq 3$ the susceptibility is almost featureless with a nearly vanishing magnitude.

To examine the corresponding real-space spin structures and their counterparts in reciprocal space we show in Fig. 4 changes of the averaged spin configuration $\langle{\bf S}_{i}\rangle$ and Bragg intensity $I^{*}(q_{x},q_{y})$ as function of temperature. For illustration we chose the value of $D=1.1$ where the transition temperature $T_{c}\simeq 0.86$ (for $T>T_{c}$ the averaged spin at a given site $\langle{\bf S}_{i}\rangle$ becomes very small, and we have to re-scale spin lengths by some factor for presentation purposes).

For temperatures below the transition temperature one can clearly see a spiral state (upper panel), with a two-peak structure of the Bragg intensity profile at (two-dimensional) vectors ${\bf Q}_{1}=(2\pi/L)(3,2)$ and ${\bf Q}_{2}=(2\pi/L)(-3,-2)$. The corresponding spiral vector is ${\bf Q}=(2\pi/L)(3,2,3)$. When increasing $T$ just above $T_{c}$ the spin pattern changes drastically (middle panel): one can see an emergence of the four-peak structure. This four-peak structure gradually transforms into a ring-shape form (bottom panel). A similar intensity patterns is also observed when projecting onto the $(q_{x},q_{z})$ plane. A close examination of the spin configuration in real space shows that there are sites with vanishing values of averaged spin $\langle{\bf S}_{i}\rangle$. Configurations with $\langle{\bf S}_{i}\rangle=0$ arise in multi-spiral states predicted for 2D models with DM interaction [14, 16, 17, 18], or for competing frustrated Heisenberg models on triangular lattices [19].

To take a close look at this multi-spiral spiral state we present $I({\bf q})$ as function of $q_{z}$ in Fig. 5. The data are shown for $D=1.1$ at three temperatures around the transition temperature $T_{c}=0.86$. The one-spiral state is easily identified at $T=T_{c}$ (upper panel), where peaks at $q_{z}=3$ and $q_{z}=-3$ correspond to ${\bf Q}$ and $-{\bf Q}$. Slightly above $T_{c}$, (middle and bottom panels) the system starts to develop two main peaks at ${\bf Q}_{1}$ and ${\bf Q}_{2}$, which signals the emergence of a quasi two-spiral state. With increasing $T$ other peaks emerge and this quasi two-spiral state gradually transforms into the ring-shape structure seen in Fig. 4 (bottom panel).

With increasing $D$ the multi-spiral state transforms into a ring-shape state of very small intensity $I({\bf q})$, which subsequently fades away. This nicely correlates with a disappearance of the hump in $C(T)$ as $D$ increases. So it is important to emphasize that the spin configuration in the hump temperature region corresponds to the above-mentioned multi-spiral structure, which melts away with an increase in temperature or $D$.

To better understand the phase transition studied here we show $I(\bf{q})$ in Fig. 6, which illustrates the drastic change at the transition temperature for the case $D=1.1$. The left panel shows $I({\bf q})$ at $T_{c}$, and the right panel just above it. The intensity profile at $T_{c}$ displays the spiral state with peaks at the spiral vector ${\bf Q}=2\pi/L(3,2,3)$ and at $-{\bf Q}$. It is seen that the intensity drops by nearly three orders of magnitude just above $T_{c}$. It results from the corresponding drop in the magnitude of $\langle{\bf S}_{i}\rangle$ and is a clear signature of a first order transition. A similar behavior of $I({\bf q})$ is observed over a large $D$ range, $0.3<D<2$. This range for $D$ also covers the region where $T_{c}(D)$ has its minimum, as shown in Fig. 7. Upon further increase of $D$, the discontinuous jump in $I({\bf q})$ at $T_{c}$ gradually vanishes, signaling that the transition becomes second order again.

The full dependence $T_{c}(D)$ is shown in Fig. 7. The data points around the minimum of $T_{c}(D)$ correspond to the first order transition in the system. Away from this region, when $T_{c}$ gradually grows, the transition becomes second order. Looking again at Fig. 6 it is tempting to consider the ratio $D/J$ as a parameter that could control the evolution of the phase transition in the helical magnets at high pressures. This supposition can probably be valid for systems with local magnetic moments (for instance, in the case of Cu2OSeO3 [20]), but not directly applicable to an itinerant system like MnSi. However, a variation of the ratio $D/J$ could be an important factor influencing the phase diagrams of various materials. Another factor that can govern the phase diagram is a magnitude of compensating terms $J^{\prime}$ and $D^{\prime}$ in the Hamiltonian $H$. In our case these terms slightly shift a position of the first-order peak in the dependence of $C(T)$ to a lower temperature, with the spin structure of the hump being not affected.

In Fig. 7 we also show the $D$-dependence of the magnitude $Q$ of the spiral vector in the ordered low-temperature phase. With increasing $D$ the spiral wave vector continuously increases, while simultaneously rotating from direction $[100]$ via $[110]$ to the body diagonal $[111]$ to better adjust the spiral to the structure of the lattice.

Finally, classical Monte Carlo simulations of a spin system with the competing interactions $J$ and $D$ reproduce very well the experimental situation observed in helical magnets MnSi and Cu2OSeO3. The hump or shoulder in the specific heat on the high temperature side of the first order peak arises as a consequence of a perturbation of the virtual second order ferromagnetic phase transition by helical fluctuations. The hump domain therefore has a very complicated spin structure (for illustration see Fig.  4) stipulated by the competing ferromagnetic and helical fluctuations. As a result, at lower temperature the system cannot proceed to an ordered state in a continuous way and takes on helical order through a first order transition. We therefore propose that the competing interactions are the primary factor responsible for an occurrence of first order phase transitions in helical magnets with the Dzyaloshinskii-Moriya interaction. Indeed, with two competing interactions and two interacting order parameters in a system one may expect a renormalisation and a sign change of the forth order term in the Landau expansion, which would lead to a first order phase transition. The arising spin fluctuations (see [7, 21]) are a direct consequence of this situation.

Conclusion

Summarizing our results we found strong evidence for a continuous transformation of a second-order transition into a first-order transition with varying $(D/J)$. At low and high $D$ the system is perturbed weakly from the respective ground state of $H_{J}$ or $H_{D}$. At intermediate values of $(D/J)\sim 1$ the system exhibits strong competing interactions resulting in a continuous change of the nature of the transition from second to first-order, with subsequent reentrant behavior when at large $D\gtrsim 2$ the system again exhibits second-order behavior. We show that the hump in the physical properties at the phase transition in MnSi and Cu2OSeO3 originates from a perturbation of a virtual ferromagnetic second order phase transition by helical fluctuations arising due to the DM interaction. In the other words the ”hump” may be viewed as a smeared out ferromagnetic second order phase transition by the helical fluctuations which eventually condense into the helical ordered phase. This conclusion completely agrees with the assumption formulated in [12].

Acknowledgements

AMB (calculations, data analysis, writing the paper) and SMS (data analysis,writing the paper) greatly appreciate financial support of the Russian Foundation for Basic Research (grant 15-02-02040, 17-52-53014), Program of the Physics Department of RAS on Strongly Correlated Electron Systems and Program of the Presidium of RAS on Strongly Compressed Matter.