B–T phase diagram of Pd/Fe/Ir(111) computed with parallel tempering Monte Carlo

We consider an ultra-thin film built from a monolayer of Pd in fcc-stacking on a monolayer of Fe in fcc-stacking on an Ir(111) surface. For this system, DFT calculations were carried out to determine the magnetic interactions. The parameters were obtained with the FLEUR ab initio package⁷ [45, 46]. Our spin models are based on the extended Heisenberg model:

$$\begin{eqnarray}&&H=-\displaystyle \sum _{{ij}}{J}_{{ij}}({{\bf{m}}}_{i}\cdot {{\bf{m}}}_{j})-\displaystyle \sum _{{ij}}{{\bf{D}}}_{{ij}}\cdot ({{\bf{m}}}_{i}\times {{\bf{m}}}_{j})+\displaystyle \sum _{i}K{({m}_{i}^{z})}^{2}-\displaystyle \sum _{i}{\mu }_{s}{\bf{B}}\cdot {{\bf{m}}}_{i},\end{eqnarray} \tag{ 2 }$$

where ${{\bf{m}}}_{i}$ is the unit vector of the magnetization of the Fe atoms at sites ${{\bf{R}}}_{i}$, J_ij are the magnetic exchange coefficients, ${{\bf{D}}}_{{ij}}$ is the DMI-vector, K the coefficient of the magnetocrystalline anisotropy and B the magnetic field perpendicular to the film. The absolute value of the magnetic moment of one Fe atom is μ_s = 3 μ_B, where μ_B is the Bohr magneton. The exchange constants ${J}_{\mathrm{DFT}}=\{{J}_{n}\}$ between the nth nearest-neighbors in the Fe layer ${J}_{1}\cdots {J}_{9}$ are given, in meV, by 14.40, −2.48, −2.69, 0.52, 0.74, 0.28, 0.16, −0.57, −0.21, respectively. The anisotropy constant is given by K = −0.7 meV, which corresponds to an out-of-plane easy axis and the DMI is restricted to nearest-neighbors with a value given by D = 1.0 meV and its direction is obtained in agreement with the model in [18]. We use these DFT parameters unless stated otherwise. When only an effective nearest neighbor exchange parameter is taken into account J_eff = 3.68 meV, K = −0.7 meV and D = 1.39 meV. Note, that these values are very close to those obtained based on experiments [31, 47]. For further details see [9, 13]. The atomistic spin simulations were performed in a super cell with N = (100 × 100) magnetic moments, where we used periodic boundary conditions in the horizontal directions. This corresponds to a super cell of (27 × 23.4 nm²) where 6 spin spiral periods are stabilized.

To study the temperature dependence of Pd/Fe/Ir(111), we perform PTMC simulations to handle the occurrence of metastable states. PTMC makes use of the temperature to exit metastable states that occur when frustration of exchange is present and isolated skyrmions or antiskyrmions appear [14, 48]. We used 240 spin configurations (replicas), which are simultaneously simulated at different temperatures distributed in a geometric temperature set [49]. The simulations were initialized with the ground state (T ≈ 0 K) configurations. We calculate the total energy E (according to equation (2)) and the magnetization density M ($M=| \sum {{\bf{m}}}_{i}| /N$) of the whole super cell of one configuration. The replicas of adjacent temperatures were swapped 10⁴ times (average steps) from which we calculate the average total energy $\langle E\rangle$ and the average magnetization $\langle M\rangle$. In the following, the braket notation is used for all averaged values. In between these swapping steps, the spin configurations were thermalized with 10⁶ Metropolis MC steps. In addition, we calculate the topological charge density $q$ by discretizing equation (1) as explained in [50]. In our atomistic model $\parallel {\bf{m}}\parallel =1$ therefore we calculate $q$ in a unit cell as:

$$\begin{eqnarray}&&\tan \left(\displaystyle \frac{1}{2}q\right)=\displaystyle \frac{[{{\bf{m}}}_{1}{{\bf{m}}}_{2}{{\bf{m}}}_{3}]}{1+{{\bf{m}}}_{1}{{\bf{m}}}_{2}+{{\bf{m}}}_{1}{{\bf{m}}}_{3}+{{\bf{m}}}_{2}{{\bf{m}}}_{3}},\end{eqnarray} \tag{ 3 }$$

where m_i is the magnetization of site i, $[{{\bf{m}}}_{1}{{\bf{m}}}_{2}{{\bf{m}}}_{3}]$ is the triple scalar product between magnetic moments at sites 1–3. We can then define unit cells in which $q$ is positive or negative. The topological charge of the whole super cell Q is then defined by the difference between the absolute values of the SkD and ASkD contributions ${Q}_{+}$ and ${Q}_{-}$, which are the integral of the positive and negative parts of q, respectively.

To distinguish between the different phases, we define the susceptibilities associated with the above mentioned quantities:

We define the heat capacity C as

$$\begin{eqnarray}&&C=\displaystyle \frac{(\langle {E}^{2}\rangle -{\langle E\rangle }^{2})}{{k}_{{\rm{B}}}{T}^{2}}\end{eqnarray} \tag{ 4 }$$

and the magnetic susceptibility as:

$$\begin{eqnarray}&&{\chi }_{M}=\displaystyle \frac{(\langle {M}^{2}\rangle -{\langle M\rangle }^{2})}{{k}_{{\rm{B}}}T},\end{eqnarray} \tag{ 5 }$$

where ${k}_{{\rm{B}}}$ is the Boltzmann constant and T is the temperature.

In analogy with the number of particles in the grand canonical ensemble [51], we define the topological susceptibility as

$$\begin{eqnarray}&&{\chi }_{\langle {Q}_{\pm }\rangle }=\displaystyle \frac{(\langle {{Q}_{\pm }}^{2}\rangle -{\langle {Q}_{\pm }\rangle }^{2})}{\langle {Q}_{\pm }\rangle {k}_{{\rm{B}}}T}\end{eqnarray} \tag{ 6 }$$

which describes the fluctuations and correlations of the SkD and ASkD.

Each order parameter is modeled by an arctangent function with a linear noise:

$$\begin{eqnarray}&&f(T)=\displaystyle \frac{I}{\pi }\arctan \left(\displaystyle \frac{2(T-{T}_{{\rm{c}}})}{w}\right),\end{eqnarray} \tag{ 7 }$$

where I is the intensity (area underneath the peak), ${T}_{{\rm{c}}}$ the critical temperature and w the mean height width. This simple model allows a precise evaluation of the critical temperature by fitting the different susceptibilites with the derivative of the arctan function i.e. a Lorentzian function of the type:

$$\begin{eqnarray}&&f^{\prime} (T)=\displaystyle \frac{I}{\pi }\displaystyle \frac{2w}{4{(T-{T}_{{\rm{c}}})}^{2}+{w}^{2}}.\end{eqnarray} \tag{ 8 }$$

This fitting procedure leads to a precision of around 10% on T_c.

First, we have performed standard Metropolis MC simulations of a spin spiral, a SkX containing 21 skyrmions and a FM state at T = 0.001 K to obtain the stability diagram in figure 1. This stability diagram shows the most stable phase for different magnetic fields at 0 K. The magnetic structures can be seen in the plots above the corresponding ground state area in figure 1. At low magnetic fields up to ${B}_{{\rm{s}},1}=1.6\,{\rm{T}}$, the system has a spin spiral ground state in agreement with DFT. The spin spiral propagates along the $\bar{{\rm{\Gamma }}}\mbox{--}\bar{{\rm{K}}}$ direction with a wavelength of λ = 4.5 nm which compares well with λ ≈ 5–7 nm obtained experimentally [30, 52]. Between 1.6 T and ${B}_{{\rm{s}},2}=2.6\,{\rm{T}}$ we find a stable SkX. The energy gain of the SkX state compared to the spin spiral or FM state is up to 0.04 meV per spin. The ordering of the skyrmions in the SkX along lines tilted to the lattice vectors was found with a first initial PTMC simulation at a magnetic field with SkX ground state. This previous PTMC simulation finds the most stable SkX, which then can be used as the starting configuration for the PTMC simulation to obtain the thermodynamical quantities. At magnetic fields above 2.6 T the FM phase is the ground state. These field values are in good agreement with previous MC simulations [32] and experimental measurements [30].

Zoom In Zoom Out Reset image size

To identify the transition temperature from the ordered, long or short range, to the PM phase we use the peak in the heat capacity C as in the work of Buhrandt et al [24]. Therefore, the temperature dependence of the total energy $\langle E\rangle$ and the heat capacity C for selected magnetic fields with different ground states (according to the stability diagram in figure 1) is shown in figures 2(a) and (b), respectively. For all magnetic fields the energy $\langle E\rangle$ in figure 2(a) increases with temperature as expected. To obtain the inflection point of this increase, we calculate the heat capacity C, which shows a peak for all magnetic fields (see figure 2(b)). The field dependence of the energy and the heat capacity is minimal, as seen in figures 1(a) and (b). We do not observe the first order peak in C that is characteristic for the Brazovskii scenario in MnSi as we have a two dimensional magnet with strong anisotropy and not a three dimensional magnet without anisotropy [24, 53]. We determine a critical temperature at about ${T}_{{\rm{c}}}\approx 214\,{\rm{K}}$.

Zoom In Zoom Out Reset image size

In a second step, we analyze the magnetization $\langle M\rangle$ in figure 3(a) and magnetic susceptibility χ_M in figures 3(b)–(e) for selected magnetic fields with different ground states according to the stability diagram in figure 1. For comparison, the magnetic susceptibility without magnetic (B = 0 T) field is replotted in figures 3(b)–(e).

Zoom In Zoom Out Reset image size

The behavior of the magnetization depends on the ground state (see figure 3(a)). Without a magnetic field (B = 0 T), the magnetization $\langle M\rangle$ is zero due to the spin spiral structure. For that case, the magnetic susceptibility shows a broad peak similar to the one reported in Rózsa et al for Pd/Fe/Ir(111) [32] and Buhrandt and Fritz for MnSi [24].

At low magnetic field in the spin spiral ground state regime (B = 1 T), the magnetization (dotted blue line in figure 3(a)) is non-zero even at T = 0 K, since the spin spiral rotates inhomogeneously due to the applied magnetic field. At T ≈ 81 K, the magnetization increases steeply and then decreases with temperature. In the magnetic susceptibility (figure 3(c)), we find an additional sharp peak at the position of the increase of magnetization.

At magnetic fields with a SkX ground state, the magnetization is non-zero at low temperature and decreases smoothly with temperature (chained red line in figure 3(a)). For this case, we find a broad peak in the magnetic susceptibility (figures 3(c) and (e)).

In the case of a FM ground state, at low temperature the system is fully magnetized (dashed green line in figure 3(a)). At about T ≈ 71 K the magnetization decreases steeply. In the magnetic susceptibility, we also find an additional peak at the low temperature side of the broad peak as in the spin spiral case with magnetic field (see inset figure 3(d)).

The above analysis of the thermodynamical quantities indicate an additional intermediate region within the ordered phase. We shade this region in gray in figures 3(a)–(e) for guidance.

At first, we analyze the topological charge $\langle Q\rangle$ in figure 4(a). Without magnetic field, $\langle Q\rangle$ remains zero. This means, the DMI itself does not create the net topological charge and it needs the magnetic field to break the symmetry in agreement with Hou et al [29]. In the case of a SkX ground state, $\langle Q\rangle$ is maximum at low temperature and decreases smoothly with temperature. When the magnetic fields stabilize the spin spiral or the FM ground states, $\langle Q\rangle$ shows a similar trend when the temperature increases. In both cases, the topological charge increases steeply at a certain temperature and then decreases smoothly again. However, the net topological charge does not indicate the possible coexistence of skyrmions or antiskyrmions which may appear in Pd/Fe/Ir(111) ultra-thin film [13, 14]. We therefore analyze in a second step the temperature dependence of both SkD and ASkD contributions.

Zoom In Zoom Out Reset image size

The SkD and ASkD contributions $\langle {Q}_{\pm }\rangle$ at B = 1 T (spin spiral ground state) are shown in figure 4(b). At low temperature, the increase of the SkD (full red line) and ASkD (dashed–dotted blue line) contributions do not differ which means that the topological charge $\langle Q\rangle$ remains constant. At high temperature, the SkD and ASkD are equal which is characteristic of the PM phase. However, in a certain temperature range (shaded gray area), the SkD is larger contribution to $\langle Q\rangle$ than the ASkD which means that on average, more skyrmions than antiskyrmions are created. Due to this difference, the topological charge is non-zero in this temperature range.

At magnetic fields with a FM ground state (not shown), the curves for the SkD and ASkD look similar to the ones of the spin spiral ground state. In the case of a SkX ground state (not shown), both densities already differ at T = 0 K, since the SkD is non-zero due to the SkX.

To visualize the intermediate region, we show in figure 5(a) snapshot of the magnetic structure and the distribution of the topological charge at the steep increase of the net topological charge at T = 72 K for B = 3 T. The magnetic structure shows 13 skyrmion-like patches (see figure 5(a)) in a magnetized background. All patches are more or less spherically symmetric and have similar sizes. The calculated net topological charge for this structure is Q = 13. The spacial distribution of the corresponding SkD and ASkD can be seen in figure 5(b). At the same position as the skyrmion-like patches in the magnetic structure, an accumulation of highter SkD is visible. We also find small ASkD contributions in the region of the patches, and both SkD and ASkD in the background, which would not be the case for condensed skyrmions in a FM background at very low temperature. In supplementary videos available online at stacks.iop.org/NJP/20/103014/mmedia, we show the evolution of these magnetic structures and spacial distributions of the SkD and ASkD as a function of MC steps for selected magnetic fields and temperatures to show that their modification changes Q.

Zoom In Zoom Out Reset image size

To analyze the fluctuations of the topological charge densities SkD and ASkD, we calculate the topological susceptibilities ${\chi }_{{Q}_{\pm }}$ which are shown in figure 6. We identify the transition temperatures ${\chi }_{\langle {Q}_{\pm }\rangle ,1}$ and ${\chi }_{\langle {Q}_{\pm }\rangle ,2}$, as the peak positions of the topological susceptibilities at fixed magnetic fields. As an example, we plot the topological susceptibilities ${\chi }_{\langle {Q}_{\pm }\rangle }$ for B = 1 T in figure 6(a). We find a broad peak with its maximum at about ${T}_{\langle {Q}_{\pm }\rangle ,2}\approx 250\,{\rm{K}}$. This peak indicates the transition from the intermediate region to the PM phase, as in the simple Heisenberg model in [42]. However, at ${T}_{\langle {Q}_{\pm }\rangle ,1}$ ≈81 K an additional sharper peak appears, which indicates the transition from the ordered phase to the intermediate region. In figures 6(b) and (c) we plot a heat map of these topological susceptibilities at different temperatures and B fields. At high temperature, both topological susceptibilities have a maximum at all magnetic fields, indicating the transition to the PM phase. At low temperature, both susceptibilities show an asymmetric behavior. At magnetic fields with a spin spiral ground state, SkD and ASkD susceptibility are increasing symmetrical with magnetic field. This does not mean, that the same amount of SkD and ASkD is created, but that the fluctuations have the same correlation length and that both skyrmions and antiskyrmions may appear.

Zoom In Zoom Out Reset image size

We now analyze the effect of the frustration of exchange interaction on the topological charge and the topological susceptibility. We therefore performed PTMC simulations at different magnetic fields with only an effective exchange coefficient J_eff with 192 temperature steps. These parameters were derived in [13] and give a similar stability diagram as in figure 1. They are also very similar to those obtained based on fitting experimental data to the micromagnetic model [31].

In figures 7(a) and (c) is plotted the topological charge $\langle Q\rangle$ for J_DFT and J_eff, respectively, for B = 2 T (magnetic field with SkX ground state) and B = 4.5 T (magnetic field with FM ground state). It shows, that in the case of a frustration of exchange (J_DFT), the maximum is in the SkX, whereas for the case of only an effective exchange parameter (J_DFT), the maximum of topological charge is in the intermediate region at a magnetic field with FM ground state.

Zoom In Zoom Out Reset image size

The figures 7(b) and (d) show the positive and negative topological charge contributions $\langle {Q}_{\pm }\rangle$ for J_DFT and J_eff, respectively, for B = 4.5 T (magnetic field with FM ground state). In both cases, the topological charges do not differ at low and high temperature, but in the intermediate region the positive charge exceeds the negative one. The inset plots show the corresponding topological susceptibilities ${\chi }_{\langle {Q}_{\pm }\rangle }$. In the case of frustrated exchange (J_DFT), we identify a broad peak indicating a transition to the PM phase for both charges at about ${T}_{\langle {Q}_{\pm }\rangle ,2}=260\,{\rm{K}}$ as in [42] and for the positive topological charge also an additional peak indicating the transition to the intermediate region at about ${T}_{\langle {Q}_{+}\rangle ,1}=100\,{\rm{K}}$. However, there is no peak indicating a transition to the interVonMalottkimediate region for the negative charge, of which the slope is smooth as a function of temperature. In the case of an effective exchange parameter (J_eff), both charges show a broad peak at about ${T}_{\langle {Q}_{\pm }\rangle ,2}=115\,{\rm{K}}$, and here also both charges show a peak indicating the transition to the intermediate region at about ${T}_{\langle {Q}_{\pm }\rangle ,1}=45\,{\rm{K}}$, whereby the peak of the positive charge is much larger than the negative one. Both peaks indicate the transition to the intermediate region (${T}_{\langle {Q}_{\pm }\rangle ,1}$).

The above results could indicate, that the frustration of exchange yield different chiral excitations. In the case of a frustration of exchange, the magnons are driven by exchange, and in the case of an effective exchange parameter driven by DMI, which would favor shorter wavelength excitations in the intermediate region and unfavor the creation of negative topological charge in agreement with the calculated energy barriers for antiskyrmion collapse of [13] when using J_DFT, the magnetic exchange interactions favor magnons of 6–7 nm periodicity and the SkX is therefore the maximum of $\langle Q\rangle$.

We summarize the information of the transition temperatures T_c, ${T}_{\langle {Q}_{\pm }\rangle ,1}$ and ${T}_{\langle {Q}_{\pm }\rangle ,2}$ in a B–T phase diagram which is shown in figure 8(a) for the case of frustration of exchange (J_DFT), where we plot the mean topological charge $\langle Q\rangle$ as contour lines. In total we identified five different phases. At low temperature three ordered phases occur with increasing magnetic field: spin spiral, SkX and FM phase. At higher temperatures (${T}_{\langle {Q}_{\pm }\rangle ,1}$) the system shows a transition into the intermediate region, with a non-zero topological charge $\langle Q\rangle$. In contrast to a previously reported phase diagram [32], we find critical temperatures which depend on the magnetic field. Here, we find a decrease of the critical temperature with increasing magnetic field in the region with a spin spiral ground state. When the FM state is lowest in energy, the critical temperatures increase with increasing magnetic field. With further increase of the temperature, the topological charge vanishes and the system becomes PM. The critical temperature of this transition to the PM state can be defined as the position of the peak of the heat capacity (${T}_{{\rm{c}}}$) or the broad peak of the topological susceptibility (${T}_{\langle {Q}_{\pm }\rangle ,2}$).

Zoom In Zoom Out Reset image size

As a last step, we calculate the B–T phase diagram when only an effective exchange parameter is taken into account (J_eff) which is shown in figure 8(b). The magnetic structures of the different ground states do not differ from the ones in figure 1, the different regimes overlap at low temperature and the B–T phase diagram shows again five phases (spin spiral, SkX, FM phase, intermediate region, PM phase). However the critical temperatures are reduced by a factor of roughly 2, compared to the case, when frustration of exchange is included (see figure 8(a)). Furthermore, the maximum of topological charge is not in the SkX phase as in the case for a full set of exchange coefficients, but in the intermediate region with a magnetic field which stabilizes a FM ground state.