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

We use an atomistic spin model derived from density functional theory calculations for the ultra-thin film Pd/Fe/Ir(111) to show that temperature induces coexisting non-zero skyrmion and antiskyrmion densities. We apply the parallel tempering Monte Carlo method in order to reliably compute thermodynamical quantities and the B-T phase diagram in the presence of frustrated exchange interactions. We evaluate the critical temperatures using the topological susceptibility. We show that the critical temperatures depend on the magnetic field in contrast to previous work. In total, we identify five phases: spin spiral, skyrmion lattice, ferromagnetic phase, intermediate region with finite topological charge and paramagnetic phase. To explore the effect of frustrated exchange interactions, we calculate the B-T phase diagram, when only effective exchange parameters are taken into account.


Introduction
Chiral magnetic spin structures, such as magnetic skyrmions [1], have received a lot of interest since they are possible candidates as bits in data storage [2][3][4][5][6]. Magnetic skyrmions are localized non-collinear spin-textures with a unique rotational sense which defines their chirality. The winding of the magnetization can be described by an integer topological charge which is called skyrmion number Q [7]: where x and y are the spatial coordinates and m the unit vector of the magnetization. Micromagnetic models have predicted the occurrence of magnetic skyrmions in condensed matter [1,8]. Density functional theory (DFT) calculations have demonstrated that skyrmion stabilization can be explained by the competition between the Dzyaloshinskii-Moriya interaction (DMI) and the magnetic exchange interaction beyond the first nearest neighbor approximation [9][10][11]. The latter is enough to stabilize topologically protected states [12,13]. Then, the stability of skyrmions and antiskyrmions [14,15] is enhanced only by DMI [16,17], which occurs when structural inversion symmetry is broken such as in chiral magnets or at surfaces or interfaces [18][19][20].
In order to use magnetic skyrmions in data storage devices, the knowledge of their temperature (T) and magnetic field (B) dependence, which can be summarized in a B-T phase diagram, is of great importance. B-T phase diagrams including skyrmion magnetic textures were first measured for MnSi bulk [21], Fe 0.5 Co 0.5 Si thinfilms [22] and in the multiferroics Cu 2 OSeO 3 [23]. In the case of MnSi, the B-T phase diagram was reproduced via Monte Carlo (MC) simulations [24]. All these diagrams show long ranged ordered phases (spin spiral, ferromagnetic (FM), skyrmion lattice (SkX)), a short range ordered phase and a paramagnetic (PM) phase which Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. are separated by critical temperatures T c which depend on the external magnetic field B. In none of the above mentioned cases, the topological charge has been measured. Recently, several works have been devoted to the analysis of the modified topological charge of skyrmions as a function of temperature and magnetic field. The topological properties of skyrmions have been analyzed at T=0 K in the context of electron transport [25]. A minimum model containing frustration of exchange interaction has been considered, showing a liquid phase of skyrmions and antiskyrmions [26]. The temperature dependence of the topological charge of skyrmions was explored via MC simulations with and without impurities [27,28], and shows decreasing charge at high temperature due to the melting of the skyrmion phase. A general B-T phase diagram of a two dimensional film of a chiral magnet obtained by MC simulations showed the occurrence of topological charge at high magnetic field in the intermediate region far from the spin spiral ground state [29].
It was shown, that magnetic exchange interaction beyond the first nearest neighbor approximation without dipole-dipole interaction [9,10,13] or first nearest exchange magnetic interaction with dipole-dipole interaction [39] was necessary to model its magnetic states. In particular, we have shown that the presence of frustrated exchange strongly enhance the barrier for skyrmion collapse into the FM state [13]. A B-T phase diagram of Pd/Fe/Ir(111) based on a Metropolis MC simulation has been reported by Rózsa et al [32]. At low temperature, the system transits from a spin spiral phase to a SkX phase and finally to a FM phase with increasing magnetic field [32], in agreement with previous works [9,10,30]. Between the ordered phases and the PM phase, a fluctuation-disordered phase was found where the skyrmion lifetime is finite [24,40,41]. In this study, the critical temperatures between the ordered phases and the fluctuation-disordered phase are independent of the magnetic field [32] which differs from previous experimental work [21,23,24].
Here, we study the temperature dependence of the contributions of skyrmion density (SkD, positive topological charge + Q ) and the antiskyrmion density (ASkD, negative topological charge -Q ) to the net topological charge Q. Q ±  are now real numbers which describe the two topological charge contributions in the whole super cell but cannot identify a particular skyrmion number [42]. We use parallel tempering MC (PTMC) [43,44] to compute Q ± to determine the B-T phase diagram of Pd/Fe/Ir(111). We show that Q ± increases with temperature and the corresponding topological susceptibility c  Q , which measures the correlations and fluctuations, can be used to obtain the critical temperatures in the whole range of magnetic fields. Here, we particularly focus on the effect of the influence of frustration of magnetic exchange interaction on the critical temperature by computing the B-T phase diagram with a full set of exchange coefficients obtained by DFT (J DFT ) which contains magnetic exchange up to the 9th nearest neighbor and an effective exchange parameter (J eff ) [13].

Computational details
To obtain the B-T phase diagram of Pd/Fe/Ir(111), we started from the electronic structure of the system to calculate the interactions between the magnetic moments with DFT calculations obtained in [9,13]. We then use these parameters in our MC simulation in order to calculate the temperature dependence of thermodynamic quantities.

Model
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 7 [45,46]. Our spin models are based on the extended Heisenberg model: 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 2 ) where 6 spin spiral periods are stabilized.

Methods
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 ) of the whole super cell of one configuration. The replicas of adjacent temperatures were swapped 10 4 times (average steps) from which we calculate the average total energy á ñ E and the average magnetization á ñ M . In the following, the braket notation is used for all averaged values. In between these swapping steps, the spin configurations were thermalized with 10 6 Metropolis MC steps. In addition, we calculate the topological charge density q by discretizing equation (1) as explained in [50]. In our atomistic model =   m 1 therefore we calculate q in a unit cell as: where m i is the magnetization of site i, [ ] m m m 1 2 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 and the magnetic susceptibility as: where k 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 which describes the fluctuations and correlations of the SkD and ASkD. Each order parameter is modeled by an arctangent function with a linear noise: where I is the intensity (area underneath the peak), T 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: This fitting procedure leads to a precision of around 10% on T c .

Results
3.1. Stability diagram at low temperature 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 1.6 T s,1 , the system has a spin spiral ground state in agreement with DFT. The spin spiral propagates along the Ḡ -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 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].

Energy and heat capacity
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 á ñ E 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 á ñ E 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 214 K c .

Magnetization
In a second step, we analyze the magnetization á ñ M 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).
The behavior of the magnetization depends on the ground state (see figure 3(a)). Without a magnetic field (B=0 T), the magnetization á ñ M 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)).

Topological charge
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 á ñ Q in figure 4(a). Without magnetic field, á ñ Q 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, á ñ Q is maximum at low temperature and decreases smoothly with temperature. When the magnetic fields stabilize the spin spiral or the FM ground states, á ñ Q 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.
The SkD and ASkD contributions á ñ  Q 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 á ñ Q 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 á ñ Q 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.
To analyze the fluctuations of the topological charge densities SkD and ASkD, we calculate the topological susceptibilities c  Q which are shown in figure 6. We identify the transition temperatures c á ñ  Q ,1 and c á ñ  Q ,2 , as the peak positions of the topological susceptibilities at fixed magnetic fields. As an example, we plot the topological susceptibilities c á ñ  Q for B=1 T in figure 6(a). We find a broad peak with its maximum at about » á ñ . This peak indicates the transition from the intermediate region to the PM phase, as in the simple Heisenberg model in [42]. However, at á ñ  T Q ,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.
3.5. Effect of magnetic exchange frustration on topological charge á ñ Q 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 á ñ Q 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. The

B-T phase diagram
We summarize the information of the transition temperatures T c , á ñ  T Q ,1 and á ñ  T Q ,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 á ñ Q 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 Q ,1 ) the system shows a transition into the intermediate region, with a non-zero topological charge á ñ Q . 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 c ) or the broad peak of the topological susceptibility ( á ñ  T Q ,2 ). 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.

Conclusion
We could show, that at low temperatures, the critical temperature from the long-range ordered phase to an intermediate region decreases in the case of a spin spiral ground state and increases in case of a FM ground state with increasing magnetic field. We showed that the pairwise creation of SkD and ASkD are good quantities to study the critical temperatures and chiral excitations. In the ordered phases, both the SkD and ASkD distribution are created at the same rate. In the intermediate region, the SkD increases faster than the ASkD which results in a net topological charge. At fixed magnetic field, in a certain temperature range, SkD and ASkD is created asymmetrically which yields a net topological charge even at magnetic fields without a SkX ground state. We showed, that the critical temperatures are increased by roughly a factor of 2, when frustration of exchange is taken into account. In addition, we find that the maximum of topological charge is in the SkX lattice phase for the full set of exchange coefficients and in contrast to that in the intermediate region with FM ground state, when only an effective exchange parameter is taken into account. We interpret these differences as the occurrence of different metastable states, depending on the stabilization mechanism. Our simulations provide quantities that could be measured that would distinguish between the two models.