Extended spin model in atomistic simulations of alloys

TL;DR

This paper introduces an extended atomistic spin model for alloys that interpolates between ferromagnetic and paramagnetic states, enabling accurate finite temperature magnetic property simulations with computational efficiency.

Contribution

The paper presents a novel extended spin model based on a parameterized Heisenberg Hamiltonian that captures magnetic behavior across temperature ranges for alloys and complex magnetic materials.

Findings

Good agreement with experimental Curie temperatures.

Effective modeling of magnetic properties in alloys like Permalloy.

Differences in magnetic moments of Fe and Ni in alloys are captured.

Abstract

An extended atomistic spin model allowing for studies of the finite temperature magnetic properties of alloys is proposed. The model is obtained by extending the Heisenberg Hamiltonian via a parameterization from a first principles basis, interpolating from both the low temperature ferromagnetic and the high temperature paramagnetic reference states. This allows us to treat magnetic systems with varying degree of itinerant character within the model. Satisfactory agreement with both previous theoretical studies and experiments are obtained in terms of Curie temperatures and paramagnetic properties. The proposed model is not restricted to elements but is also applied to binary alloys, such as the technologically important material Permalloy, where significant differences in the finite magnetic properties of Fe and Ni magnetic moments are found. The proposed model strives to find the right compromise between accuracy and computational feasibility for accurate modeling, even for complex magnetic alloys and compounds.

Figures

Tables

Table 1. Table 1: Curie temperatures T c in K from calculations and comparison to experiments (Expts.) and previous studies (Theory). N (LSF) denotes calculations without (with) longitudinal fluctuations and rescaled exchange parameters. DLM and FM denotes from which reference state the exchange interactions are calculated.

Tc (K)	DLM				FM				Refs.
	LDA		GGA		LDA		GGA
System	N	LSF	N	LSF	N	LSF	N	LSF	Expts.	Theory
Fe	1132	926	1275	1068	1238	939	1412	1083	1043a	1060b,1065c,900d,1095e
Co	1441	807	1611	979	1475	791	1637	902	1388a	1080b,1280d,1012e
Ni	509	432	580	478	486	431	578	459	633a	510b,615c,412e
Py	645	482	747	572	706	484	761	560	850f	650g
Fe0.3Co0.7	1564	927	1686	1105	1605	902	1733	1080	-	1490d
Fe0.5Co0.5	1650	1005	1862	1215	1634	982	1885	1165	1253-1370h	1600d
Fe0.7Co0.3	1689	1060	1879	1267	1656	1010	1847	1211	-	1490d
a Reference Wohlfarth,1980
b Reference Rosengaard and Johansson,1997
c Reference Ruban et al.,2007
d Reference Ležaić et al.,2007
e Reference Uhl and Kübler,1996
f Reference Adachi et al.,1991
g Reference Kudrnovský et al., 2008 (RPA)
h Reference Kawahara et al.,2003

Table 2. Table 2: Spin stiffness D in (meVÅ 2 ) from calculations and comparison to experiments (Expts.) and previous studies (Theory). DLM and FM denotes from which reference state the exchange interactions are calculated and LDA and GGA denote different exchange-correlation potentials in calculations.

D (meVÅ2)	DLM		FM		Refs.
System	LDA	GGA	LDA	GGA	Expts.	Theory
Fe	320	368	466	573	314a	247b,250c
Co	614	676	675	723	510d	502b,663c
Ni	895	924	707	924	550e	739b,756c
Py	655	693	620	611	390f	-
Fe0.3Co0.7	552	600	611	701	476g	-
Fe0.5Co0.5	536	588	526	646	800g	-
Fe0.7Co0.3	496	536	466	472	470g	-
a Reference Stringfellow,1968
b Reference Rosengaard and Johansson,1997
c Reference Pajda et al.,2001
d Reference Wohlfarth,1980
e Reference Mitchell and Paul,1985
f Reference Adachi et al.,1991
g Reference Liu et al.,1994

Table 3. Table 3: The averaged tilted angle between the nearest neighbours θ N ​ N subscript 𝜃 𝑁 𝑁 \theta_{NN} at T c subscript 𝑇 𝑐 T_{c} and at T = 1.25 ​ T c 𝑇 1.25 subscript 𝑇 𝑐 T=1.25T_{c} and correlation length κ 𝜅 \kappa at T = 1.25 ​ T c 𝑇 1.25 subscript 𝑇 𝑐 T=1.25T_{c} . Our calculation including (without) LSF are denoted LSF (N) and compared to previous theory (Ref.) and experiments (Exp.).

	${\theta }_{NN}$ (∘)		$\kappa$ (Å-1) at $T=1.25T_c$
System	$T_c$	1.25 $T_c$	LSF	N	Ref.	Exp.
Fe	70	81	0.40	0.36	0.38a	0.40 b
Co	69	85	0.38	0.25	0.34a	-
Ni	78	86	0.24	0.28	0.28a	0.24 c
Py	76	85	0.26	0.27	-	-
a Reference Rosengaard and Johansson,1997
b Reference Shirane et al.,1986
c Reference Steinsvoll et al.,1984

Equations

$$H_{tot}=H_0+H_{LSF}(\{|{\mathbf{m}}_i|\})+H_2({\mathbf{m}}_i,{\mathbf{m}}_j),$$

$$\begin{array}{rl}{H}_{tot}=& \sum _i{J}_1(\{|{\mathbf{m}}_i|\})\\ & -\sum _{ij}(1+k)J_{ij}(|{\mathbf{m}}_i|,|{\mathbf{m}}_j|){\mathbf{m}}_i\cdot {\mathbf{m}}_j.\end{array}$$

$$Z=\int g(m)\mathrm{d}m\exp \left(-\frac{H_{tot}}{k_{B}T}\right),$$

Full text

I Introduction

First principles calculations and atomistic modelling play an important role for understanding the properties of magnetic compounds and alloys. In particular, the finite temperature properties determine the usefulness of the material for technological applications. However, a complete theory for finite temperature magnetic properties remains a formidable challenge in condensed matter physics despite the vast advances in the recent decades. Such a theory need to include not only transversal but also longitudinal spin fluctuations (LSF) of the magnetic moments in responsible for the moment formation in the material. Magnetic moments can roughly be classified in two different types, either being well localized, or exhibiting a more itinerant character. Previous studies have typically described localized moment systems by means of the Heisenberg modelHeisenberg (1928) and itinerant magnets by band theoriesBloch (1929); Stoner (1936); Slater (1936a, b); Stoner (1938), respectively. However, several effortsMoriya (1979); Uhl and Kübler (1996); Rosengaard and Johansson (1997); Ruban et al. (2007) have been directed towards a unification of the existing approaches with the goal of being able to describe the two different scenarios within a single spin model.

The first step was taken by MoriyaMoriya (1979) describing the magnetic phase transition based on low-energy excitations. In the frame work of Ginzburg-Landau theory, the Hamiltonian is expanded in the even powers of the local moment magnetization $\mathbf{M}(\mathbf{r})$ where the expansion coefficients are obtained by calculating the constrained total energy $E(M)$Mohn et al. (1990). This theoretical effort was followed by Uhl and KüblerUhl and Kübler (1996), who formulated a spin fluctuation theory in which the exchange parameters of the Hamiltonian were calculated from first-principles of spin-spiral states. Inspired by this work, Rosengaard and JohanssonRosengaard and Johansson (1997) proposed a Hamiltonian where a LSF term was introduced in addition to the Heisenberg term and subsequently calculated thermodynamic properties by means of Monte Carlo simulations. Ruban et al.Ruban et al. (2007), developed an extended scheme where a more advanced treatment of the LSF was introduced. In this latter approach, the disordered local moment was used as the reference state of mapping the total energies and thus the model can be expected to give a more accurate description of the paramagnetic state compared to the FM ground state.

A large number of works have later been focused on the effect of longitudinal fluctuations in the description of the half metallic properties of the Heusler alloysLežaić et al. (2006, 2013); Sandratskii et al. (2007); Sandratskii (2008). Common to these works is the proposal of a model Hamiltonian that considers a separate treatment of the strong-moment components and the induced moment, such that it can effectively take into account the fluctuations on the induced moments by renormalization of the exchange interactions. A similar technique has been applied to systems with multiple sublattices in other studiesPolesya et al. (2010, 2016). A further extension includes calculation of the ordering temperature not only to the ferromagnetic-nonmagnetic transition but also for the Néel temperature of antiferromagnetic systemsPolesya et al. (2016); Khmelevskyi (2016a). Parallel to the studies of equilibrium magnetic properties from Monte Carlo methods, there has recently been attempts to include longitudinal fluctuations also in dynamical simulations in the framework of Langevin spin dynamicsMa and Dudarev (2012, 2014). However, so far it has only been applied for simple elemental metals as Fe and Gd.

Experimentally, the degree of itinerancy is most frequently studied in the paramagnetic susceptibility, either by comparing the effective moment extracted from the Curie-Weiss constant with the true local moment, or the deviation of the Curie-Weiss law. Recently, a great deal of progress has been achieved to study the Stoner excitations through the transversal dynamical susceptibility from first principles within the framework of linear response time dependent density functional theory Buczek et al. (2011); Lounis et al. (2015), or from tight binding theory Costa et al. (2004) . However, due to the complexity of such calculations one faces with difficulties applying the theory to complex magnetic compounds and alloys and to include temperature effects. An issue that often occurs when mapping total energies from first principles to a spin model is that the calculated exchange interactions sometimes render distinctly different results depending on the starting reference state of the mapping. Hence, it would be very attractive to construct an extended spin model where the reference state dependence is eradicated and at the same time include effects of LSF.

In this paper, we propose a general extension of the Heisenberg model that indeed removes the dependence on the reference state and includes not only the low energy transversal excitations but also the longitudinal spin fluctuations. The scheme combines first principles and Monte Carlo (MC) simulations and provides a more rigorous description of the finite temperature magnetic properties compared to models that includes only transversal fluctuations. In particular, it is expected that a wider range of materials can be investigated within our model, including materials with induced moments.

The paper is organized as follows. Our extended spin model and the procedure of calculating the parameters from first principles are presented in Sec. II.1. The modified algorithm in the Monte Carlo for updating the moments in each step is clarified in Sec. II.2. Some practical issues including the approximations made from the mapping of the system on the extended spin model are discussed in Sec. II.3 and III.1. In the numerical results, Sec. III.2 we compare the critical temperature of each of the considered systems and study the influence of imposing the LSF term. The energy and moment distributions as functions of temperature are shown in Sec. III.4 and Sec. III.5. Paramagnetic properties, specifically the magnetic short range order (MSRO) are discussed and compared with other works in Sec. III.6. Sec. III.7 contains technical remarks and comparisons to existing methods found in the literature and we conclude and summarize our work in Sec. IV.

II Formalism

II.1

Formulation of the problem

A magnetic material at finite temperatures displays two different types of excitations, longitudinal and transversal fluctuations, as illustrated in Fig. 1.

The starting point for the simulations is the construction of a model Hamiltonian that includes both the longitudinal and transversal spin fluctuations of the magnetic moments $\mathbf{m}$, such that the total Hamiltonian can be expressed as

[TABLE]

where $\mathbf{m}_{i}$ denotes magnetic moment at site $i$, $H_{0}$ is the nonmagnetic reference energy, $H_{LSF}$ and $H_{2}$ the energies (Hamiltonian) of longitudinal and transversal spin fluctuations, respectively. The nonmagnetic reference energy can be disregarded from the problem since it only sets a reference energy which does not change during the simulations. The proposed form of the Hamiltonian, Eq. (1) takes the form

[TABLE]

In Eq. (2), the first term describes the total energy of the homogeneous paramagnetic state which in the model corresponds to the LSF energy that in turn depends on the parameter $J_{1}$ which depends on the magnitudes of the local moments of each magnetic atom type of the system. The second term is the standard Heisenberg Hamiltonian of transversal fluctuations with the important distinction that the exchange parameters $J_{ij}$ now are explicitly dependent on the size of the magnetic moments and “rescaled” with the parameter $k$, as will be explained below. All the parameters entering the Hamiltonian are in our model parameterized from first principles calculations from which the total energy is mapped to the Hamiltonian. When dealing only with transversal fluctuations, a common weakness is the reference state dependence of the calculated exchange interactionsHeine et al. (1990). This weakness does not come as a surprise since the Heisenberg model assumes rigid moments where the size of the magnetic moment does not change upon rotation, i.e. strictly speaking only valid for localized spins. Unfortunately, most of the real materials possess some degree of itinerancy where the size of the moment changes depending on internal orientation rendering Heisenberg model less applicable. For instance, one severe and common example is Ni where the magnetic moment vanishes for tilting angle between two nearest neighbour moments of around 60 degrees or more Rosengaard and Johansson (1997); Halilov et al. (1998).

Another large important class of materials for which a pure Heisenberg description fails to describe the properties correctly are systems which exhibit induced moments, i.e. moments that only appear in the presence of other localized moments creating a local internal magnetic field that polarizes the atoms. Even for systems where the Heisenberg description works reasonablely well, such as Fe, the calculated exchange interactions do depend on the reference state Heine et al. (1990); Szilva et al. (2013) and hence the calculated critical temperatures could differ to some extent. Szilva et. al.Szilva et al. (2013) proposed a formalism that allows the calculation of the exchange interaction from a general non-collinear magnetic state that better represents the configurations at finite temperatures by adding a biquadratic term. Here, in this study we take an alternative approach to soften the reference state dependence and at the same time include longitudinal spin fluctuations. This is achieved by using the extended spin Hamiltonian, Eq. (2) and employing two reference states in the mapping, as illustrated in Fig. 2. The most natural reference states are the ferromagnetic (FM) low temperature state and the paramagnetic high temperature state. The latter is referred as the disordered local moment (DLM) state. Both states are naturally treated in a straight forward manner in the present computational framework.

The proposed Hamiltonian Eq. (1), is designed in such a way that it is exact for both the FM and DLM states, in the sense the total energies from DFT are recovered, and it interpolates between the two limits at intermediate temperatures. This is most easily realized in the DLM limit where the Heisenberg term is equal to zero such that the only surviving term is the sum of $J_{1}$, consistent with the model proposed by Shallcross et. al Shallcross et al. (2005) and Ruban et. al Ruban et al. (2007). The parameter $k$ can then be viewed as a rescaling parameter that tunes the interpolation between the DLM and FM limits. However, it is worth stressing that the model does not contain any free parameters. Calculating exchange parameters starting from either FM or DLM state will indeed result in different sets of $J_{ij}$. However, since the overall Hamiltonian is fixed, it will result in different values of the rescaling factor $k$ such that the “rescaled” exchange parameters as well as the exchange energy are similar, independently from the starting reference state. The most accurate LSF term in the Hamiltonian describes the energy cost or gain of forming a local magnetic moment with size $|\mathbf{m}|$ embedded in a host with the average moment $\langle m\rangle$. This choice was employed by Shallcross et. al Shallcross et al. (2005) and Ruban et. al Ruban et al. (2007) but has the main disadvantage that it becomes prohibitively complicated for systems with several magnetic components, such as binary and ternary alloys. An approximation to simplify the model is to assume that the embedded moment has the same size as the host, in such a way that the number of indices is reduced to the number of components without the additional index of the average moment. For one atom per cell, this becomes the well known Landau expansion, in which the energy is expanded in even powers of the size of the magnetic moment $|\mathbf{m}|=m$, i.e. $E(m)\approx a_{1}m^{2}+a_{2}m^{4}+...$, where $a_{i}$ are expansion coefficients. By employing constrained density functional theory using the Fixed Spin Moment (FSM) technique, $E(m)$ is easily calculated and analytically fitted to the Landau expansion expression. Even with this simplified LSF term in the Monte Carlo simulations, a number of studiesRosengaard and Johansson (1997); Ležaić et al. (2006); Ma and Dudarev (2012) has been demonstrated with rather satisfactory results for systems with one atom per cell, such as Fe.

Alternative methods to treat longitudinal fluctuations include the spin cluster expansion (SCE) technique Drautz and Fähnle (2004), in which the magnetic energy from DFT is rigorously expanded in series of multispin interactions. For the case of a weak induced moment in presence of a strong moment, specialized schemes have been designed by Lezaic et. al Ležaić et al. (2013) and Polesya et. al Polesya et al. (2010).

However, as far as we are aware, there has not been any systematic attempts to expand the formalism to the general case with several magnetic components, which is the main motivation for the present work. To be more specific, we here expand the formalism to binary alloy systems. This covers many of the technologically important alloys but also the essential class of materials consisting of strong and weak magnetic moments, where the weak moment is induced by the strong moment from an internal magnetic field. The extension to several magnetic components is straight-forward, although the computational effort drastically increases. Nevertheless, by employing modern computational techniques as discussed below it is feasible without too much computational effort. The crucial step in the parametrization of the full LSF energy surface is to constrain the size of the magnetic moment on one component at the time. For simplicity, we restrict the discussion to binary alloys with magnetic components $\mathbf{m}^{\alpha}$ and $\mathbf{m}^{\beta}$, respectively. This could, as an example, represent the Fe and Co magnetic moments in a FeCo alloy. To obtain reasonable numerical accuracy, around 20 different moment sizes of each component are required, such that the total number of configurations is $20^{2}=400$ to fully parametrize the LSF energy surface. However, each configuration is independent of each other, so clever scripting and making use of the harvesting power of modern parallel computers makes the parametrization efficient despite the vast number of total configurations required. The moment sizes are arranged in a $n\times m$ Cartesian grid with indices $(m_{i}^{\alpha},m_{j}^{\beta})$, where $\alpha=\{1,..n\}$ and $\beta=\{1,...,m\}$. n and m are the number of moment sizes performed for $m_{i}$ and $m_{j}$, respectively. Indices $i$ and $j$ denote the site indices. Electronic structure calculations based on density functional theory for each grid point yield the input parameters $J_{1},k$ and $J_{ij}$ in the Hamiltonian, Eq. (2), for that particular grid point.

II.2 Atomistic simulations including longitudinal moment fluctuations

Following the full Hamiltonian parametrization, one can use a statistical treatment, for instance Monte Carlo (MC) simulations to work out the equilibrium properties or atomistic spin dynamics (ASD) for the dynamical properties. For systems with one magnetic atom per cell, MC simulations have been performed in a number of previous studiesRosengaard and Johansson (1997); Ležaić et al. (2006); Ruban et al. (2007); Ležaić et al. (2013) while ASD simulations are much less exploredMa and Dudarev (2012, 2014). In this study, we restrict ourselves for the moment to MC simulations while ASD simulations are left for a future study. Using the Metropolis algorithm, each trial move of every selected atomic moment not only consists of the normal rotation but also a random change of the magnitude of the moment. Since the Hamiltonian is only defined at the grid points, we have employed bilinear interpolation of the intermediate points in the simulations. As mentioned above, we are using a slightly simplified LSF Hamiltonian and it has some practical consequences. When the magnitude of a trial moment is changed in the Metropolis update, i.e. $m_{i}^{\alpha^{\prime}}=m_{i}^{\alpha}+\delta m$, we make the approximation that we use the average magnitude $\langle m^{\beta}\rangle$ of the other component in the binary alloys around the site $i$, such that the parameters $J_{1},k$ and $J_{ij}$ yield the configurations corresponding to $(m_{i}^{\alpha},\langle m^{\beta}\rangle)$ and after the change $(m_{i}^{\alpha^{\prime}},\langle m^{\beta}\rangle)$, respectively. This approximation does not cause any significant issues for the homogeneous systems since the magnitude of magnetic moments in average does not change much within a sphere from the trial atom point of view. Comparing this approximation against doing interpolation of $J_{ij}$ between each pair $m_{i}$ and $m_{j}$, only minor differences were found but it gains a magnitude faster in the computational efforts. The $J_{1}$ term may be a bit more sensitive since it is assumed that moments belonging to the same component $\alpha$ change the moment size instantaneously with the trial moment. Comparison of different strategies are discussed in more detail in Sec. III.7.

II.3 Details of the calculations

Electronic structure calculations based on density functional theory (DFT) were employed for extracting all material dependent parameters in the Hamiltoninan. We used an implementation of multiple-scattering theory based on the Korringa-Kohn-Rostoker (KKR) formalism Ruban et al. (2002). We employed both the local spin density approximation (LDA) and the generalized gradient approximation (GGA) using the PBE parametrizationPerdew et al. (1996) as the exchange-correlation potential. Another method normally referred to as Fixed Spin Moment (FSM) Dederichs et al. (1984); Schwarz and Mohn (1984) where the magnitude of the magnetic moment was varied by employing a small constraining field with strength adjusted in such a way that it reaches the prescribed target. Around 20 different moment sizes were used for each magnetic component. The chemical and magnetic disorder were treated using the coherent potential approximation (CPA). A binary alloy with components A and B can be written as $A_{x}B_{1-x}$, where $x$ is the concentration of A and this configuration is used for the ferromagnetic state (FM). The disordered local moment (DLM) configuration is obtained by forming a four-component alloy, $A^{\uparrow}_{x/2}A^{\downarrow}_{x/2}B^{\uparrow}_{(1-x)/2}B^{\downarrow}_{(1-x)/2}$, with $\uparrow$ ($\downarrow$) denoting moments pointing along (opposite) to the quantization axis. The exchange interactions $J_{ij}$ were calculated using the LKAG formalism Liechtenstein et al. (1984, 1987). The experimental lattice constants were used for all systems neglecting any thermal expansion. Moreover, the electronic entropy effects were neglected to avoid additional temperature dependence of the parameters.

The atomistic Monte Carlo simulations making use of the Metropolis algorithm were performed and implemented in the Uppsala Atomistic Spin Dynamics (UppASD) software package Skubic et al. (2008); Eriksson et al. (2017). In each trial move in the Metropolis update, in order to speed up the calculations, one of the three modes was randomly selected between a pure transversal, a pure longitudinal or a combination of the two. The transversal updates were in turn chosen randomly to either an uniform rotation, a Gaussian distribution rotation or a spin flip for further speed up. Critical temperatures were extracted using the fourth-order size-dependent Binder cumulantBinder (1981) and since we are dealing with random alloys, for each system and temperature, an ensemble averaging over around 10 different disorder configurations were performed. Calculation of spin stiffness followed the procedure outlined in Ref.[Pajda et al., 2001] generalized to random alloys using the same methode as in Ref.[Dürrenfeld et al., 2015] averaged over 1000 disorder configurations.

Evaluation of thermodynamic properties for a model including longitudinal fluctuations involves technically intricate details. The longitudinal part of the partition function over which a functional integration is performed can be written asMurata and Doniach (1972); Sandratskii (2008); Wysocki et al. (2008); Dietermann et al. (2012); Ruban et al. (2013); Khmelevskyi (2016b)

[TABLE]

where $g(m)$ is the phase space measure (PSM)Wysocki et al. (2008). In this study we use $g(m)=1$, that is the so-called uniform PSM or Murata-Doniach metric that corresponds to the decoupled treatment of longitudinal and transversal spin fluctuations. However, as was pointed out in Ref.[Wysocki et al., 2008], magnetic moments are not canonical variables and therefore the PSM is not known. Including PSM with Jacobian factor, $g(m)=m^{2}$, couples the longitudinal and transversal spin fluctuations and it has been argued that this choice could lead to improved resultsSandratskii (2008); Ruban et al. (2013). Although not used for the majority of simulations in this study, we did implement this PSM in our MC program and we will briefly discuss some results and differences to the uniform PSM in Sec. III.7.

III Results

For the results in this section, we consider the elemental materials Fe, Co and Ni, where Fe crystallizes in the bcc structure and Co and Ni in fcc structure. Moreover, we investigate the Permalloy in the fcc structure, which is a binary alloy with composition 81%Ni and 19% Fe. Finally, we investigate a binary alloy of Fe and Co in the bcc structure, with Co concentrations ranging from 30% to 70%.

III.1 The total energy landscape from the first-principles calculations

The energy surfaces relevant to the parameters $J_{1}$ and rescaling factor $k$ are shown in Fig. 3 for Ni. We show that the starting point of using different sets of exchange parameters from DLM and FM state can be eliminated using the rescaled interactions. Figure 3a and 3b indicate that the total energy of FM state has a shallow minimum of the moment size around 0.7 $\mu B$, while there is no such a minimum in the DLM state, indicating that in the DLM state no magnetic moments can be sustained unless additional fluctuations are present in the system. The total exchange energy is not explicitly plotted here, but from Eq. (2) it is apparent that it is the total energy difference between the FM and DLM states. Inspection of individual nearest neighbour exchange interaction corresponding to the total energy minimum for FM, it is found that the DLM reference state gives approximately two times larger $J_{ij}$ than in the FM state. This fact is well reflected on the calculated Curie temperature based on the traditional Heisenberg model. In Fig. 3c and 3d, the rescaling factor $k$ of the exchange interactions are displayed in the case when the “bare” exchange interactions are calculated from the DLM and FM states, respectively. A negative (positive) value of $k$ corresponds to reduction (addition) of the exchange such that it forces the exchange energy being the same regardless of chosen initial reference state. Calculating $J_{ij}$ from the FM state gives a larger surface roughness of $k$ compare to calculation from the DLM state due to larger numerical sensitivity, however since the overall total energy surface is smooth, it is not expected to cause any particular numerical difficulties in the simulations. Although we show only Ni as an example of the the energy parameterization, the procedure is similar for the other systems with the key feature that the overall spin model gives correct total energy in both the FM and DLM limits.

III.2 Critical temperatures and spin stiffness

In Tab. 1, our calculated Curie temperatures, $T_{c}$ together with previously published experimental and theoretical results are collected. In all of our calculations, the rescaled exchange parameters were used and we compare the two cases when LSF is included or not in the Hamiltonian from Eq. (2). A few general trends of the calculated results are immediately noticed. First of all, the inclusion of LSF always lowers the calculated $T_{c}$ values and underestimates them compared to experiments, with the exception of Fe where excellent agreement is found. Secondly, within the same functional (LDA or GGA), there is only a minor difference in $T_{c}$ starting from either DLM or FM configuration, which is one of the design goals of the model. Without LSF and rescaling of the exchange parameters, the well known dependence of chosen reference state is evident. For example, calculated $T_{c}$ of Ni in GGA is around 400 K (900 K) with “bare” exchange parameters calculated from FM (DLM) reference, compared to around 580 K using rescaling for both reference states. For the DLM state, the moments were fixed to the same value as in FM, since without constraints or inclusion of LSF term the DLM moment vanishes for Ni. In general, for all the considered systems, without rescaling the exchange interactions from FM configuration underestimate $T_{c}$ while there is an overestimation using exchange interactions from DLM. Within the same level of approximation and methodology, our calculations corresponds in general rather well with previous published results.

The quantitatively underestimation of $T_{c}$ with LSF included is however slightly disappointing but could be understood in the present model from relatively simple arguments. If the LSF energy would not manifest a minimum, as in Ni for example, or having a shallow minimum like as in Co, then the system may gain energy of by reducing the moment at finite temperature. However, at the same time, the exchange energy is decreasing upon decreasing moment and fluctuations. On the opposite, increasing moment size gains exchange energy but costs LSF energy, so there is a delicate energy balance between the two contributions. Consequently, with LSF included it does provide the system an alternative path of lowering the energy thus making the system more magnetically “soft”. It could well be that the softening of the moment in the present model is slightly overestimated since the LSF energy is only an approximation to some extent and depends on the average magnetic moments of the local environment. Moreover, using the uniform phase space measure as done here causes an oversampling of configurations with smaller moments due to the decoupling of the transversal and longitudinal fluctuations.

FexCo1-x alloys in bcc structure are known for having a maximum average moment around $x\approx 0.7$ according to the Slater-Pauling relationship, which is found both in experiments Adachi et al. (1991) and from calculationsAbrikosov et al. (1996). These alloys have among the largest Curie temperatures of any discovered transition metal alloy systems, even larger than the individual host elements Fe and Co. This is mostly confirmed by our calculations, with a possible exception of Fe0.3Co0.7 in GGA that it is very close to the elemental Fe. On a technical note, the results of Fe, Co and Ni in Table 1 were obtained treating the system as having two components. This makes the evaluation of the exchange and LSF energies slightly more accurate than treating the system as having single component, despite describing the same element. A comparison between the two are discussed in more details in Sec. III.7.

In Tab. 2, results of our calculated values of the spin stiffness together with experimental and previous published values are compared. Spin stiffness is sensitive to the individual exchange interactions just as the sum of all exchange interactions gives the exchange strength which is relevant to the Curie temperature. The difference being that the spin stiffness contains an additional spatial dependence in the sum. The calculated spin stiffness values are in general in reasonable agreement with both experiments and previous calculations for systems where such data exist. Fe is found to have too large spin stiffness when the exchange parameters are calculated from the FM configuration. This is due to the fact that the individual exchange parameters are not only rather different in DLM or FM configuration, but also rather sensitive to the volume and the basis set in this case. Although the exchange interactions are scaled uniformly to yield the correct total energy, that does not universally translate to improving the spin stiffness. The values of Ni are also overestimated compare to experiments, as found in previous calculations, but the most likely explanation is the use of LDA or GGA that fails to fully describe the electronic structure of Ni being a moderately correlated metal.

In the case of random alloys, Py and FexCo1-x, the spin stiffnesses were calculated as an average over 1000 different disorder configurations. Perhaps not so surprising, the spin stiffness for Py is rather similar to Ni while values for the FexCo1-x alloys are similar to those of elemental Fe and Co. Regarding experiments, inelastic neutron scattering experiments is considered the most straight forward way to measure spin waves and thus the spin stiffness. For the elemental materials Fe, Co and Ni, there are plenty of existing experimental studies and they all rather much agree to each other. However, for random alloys such as Py and Fe-Co, we find only very few studies of direct measurements of $D$. For the measured Fe-Co systems, a thin film geometry was used which could affect the results and that could in turn explain the difference from the calculated bulk value.

III.3 Finite temperature transversal and longitudinal fluctuations

To further elaborate on the temperature induced magnetic phase transition, we show in Fig. 4 the total energy decomposition in terms of the exchange and LSF energies, as a function of the temperature. The exchange term gives the energy gain by aligning moments with each other, while the LSF energy is a measure of the formation (reduction) of a local moment. By resolving the exchange and LSF energies, the stability of the local moment at a certain temperature can be estimated. In the case of Fe, Fig. 4a, the LSF energy is more or less constant, while the exchange energy is monotonous increasing with respect to the temperature and vanishes at high temperature where all moments are randomly distributed. The energies confirm that Fe is forming a stable moment in both the FM and DLM states. In contrast to Fe, Ni as shown in Fig. 4b reveals a much more precarious energy balance. Upon increasing the temperature to around $T_{c}$, a part of the exchange energy increase is compensated by decreasing the LSF energy, causing an reduction of the magnetic moment. For Py, shown in Fig. 4c, the interplay of the LSF and exchange energies behaves similar as to elemental Ni where part of the loss of the exchange energy is compensated by the LSF energy. As a result, the moments of Py in average have a tendency to be reduced at finite temperatures.

III.4 Energy distribution

In the previous section we discussed the average decomposed energies as function of temperature. Here, we study in more detail the distribution of energies during the simulations.

In Fig. 5, the composition of site resolved energies of Py from a snapshot during the Monte Carlo simulations is displayed for two temperatures, at low temperature ($10\ K$) and around Curie temperature ($550\ K$, see Tab. 1). If the energies are averaged both over all the atoms and in time, the results in Fig. 4 are recovered. The spread of the magnitude of moments are discussed in more details in Sec. III.5.

At low temperatures, the thermal fluctuations are very small and consequently the spread of energies are small. The total energy is fluctuating around the minimum of the energy surface and the largest difference between Ni and Fe moments is the larger exchange energy associated with Fe. At elevated temperature close to the Curie point the situation is rather different. First of all, due to the large thermal fluctuations the spread of both energies and moment magnitudes are significant. Also the directions of moments have large spread and since the figure shows a particular snapshot, there are certain moments that are anti-aligned to each other causing positive site exchange energies. All these fluctuations makes the system different from pure DLM configurations where the moment distribution are random but with fixed magnitude of moments. The high temperature snapshot configurations could be viewed as many DLM-like configurations superimposed on each other. This has also been explored directly from DFT calculations in other studies. For instance, in Ref.[Dong et al., 2015] DLM configurations combined with fluctuating moments and statistical models were used to calculate high temperature properties such as elastic constants and average magnetic moment in Fe. As a contrast, in the simulations employed here, all these fluctuations are naturally arising as outputs. The site resolved LSF energy shows very different behaviour for Ni and Fe atoms. As in DLM, the positive curvature of Ni means that it prefers not to sustain any magnetic moments on its own, nevertheless any finite Ni moments are stabilized by the entropy of spin fluctuations. On the other hand, the LSF energy of Fe has a minimum around 2.8 $\mu_{B}$ and as a consequence, the Fe moments are magnetically more resilient than Ni, and which in turn are induced from the Fe moments. It is worth noting that the individual output LSF energies on both Fe and Ni in Py also follow the Landau-type expansion as expected from the model.

III.5 Average magnetization, local moment and moment distribution

The extended Heisenberg model where LSF is taken into account allows us to examine the dependency of the change in the local moment size on temperatures by performing the MC simulations. The magnetization and the local moments size with respect to temperature are displayed in Fig. 6(a) for Fe, Co and Ni. We observe a linear reduction on the magnetization at low temperature in all the considered cases which does not follow the Bloch’s $T^{3/2}$ as expected from using classical statistics. However, here we are focusing at the high temperatures regime for temperature around $T_{c}$ and above where the classical spin model is well justified. The average moment size of Fe is reduced only weekly of about $10\%$ up to $T_{c}$. In contrast, the reduction of the average moments of Co and Ni in percentage increase by a factor of two compared to Fe. Our results are consistent with what have been found in Ref.[Rosengaard and Johansson, 1997] where the changes on the Ni size is continuous at Tc, however a different behaviour for Ni was shown in Ref.[Ruban et al., 2007]. Despite of a rather constant average moment size, the variation on the level of individual local moment sizes is rather different. This is shown in the Fig. 6(b), where we have calculated the standard deviation of the local moment size as a function of temperature. Compare to Fe and Co where the standard deviation are relatively small at low temperatures, Ni has a large standard deviation even at low temperatures. This fact indicates that the description of a “rigid moment” in the longitudinal direction is appropriate for Fe, intermediate for Co but not so for Ni, as expected since Ni moments being the most itinerant of the three.

In Fig. 7(a), the component-resolved local moment size and magnetization are calculated for Py. We focus on a few intricate features that could only be displayed in random alloys. First of all, the size of Fe moments in Py is calculated to be larger than it is in the elemental Fe while the moment size of Ni has roughly the same value, as expected from other studiesPan et al. (2016). Moreover, compared to the elemental Fe where moment sizes are rather constant, the local Fe moments in Py are more strongly reduced, in particular above $T_{c}$. Fig. 7(b) the moment size distributions of Fe and Ni in Py are shown for two temperatures, $T=10\ K$ and $T=550\ K$. At $T=10\ K$, the longitudinal fluctuations are small for both Fe and Ni resulting in narrow distribution. At elevated temperatures, the longitudinal fluctuations are more energetically favored. At $T=550\ K$, right below the $T_{c}$ of the system, the moment distribution of Fe shows an asymmetry with a long tail to small moments. The distribution profile of Ni is rather different than Fe. It has a significant shift to small moments indicating that Ni prefer to eventually lose its moment.

III.6 Impact of LSF on the high temperature properties

In Fig. 8 we show the magnetic short-range ordering (MSRO) in terms of the average tilting angle between the nearest neighbour moments as a function of the temperature relatively to the ordering temperature $T_{c}$. It is found that above $T_{c}$ the LSF term has negligible effect on the tilting angle, moreover it is only determined by the underlying structures, which can be seen from the fact that all systems in fcc lattices the temperature dependent $\theta_{NN}$ is essentially the same. Our results confirm the conclusion in Ref. [Wysocki et al., 2008] where it was reported that the MSRO is expected to be weak for spin models of similar type as we employ in this study. In fact, the degree of MSRO is expected to be more or less constant for a large range of itinerancy. To be more specific, we obtain $cos\theta_{NN}=0.2$ at $T=1.1T_{c}$ for Fe which is in an excellent agreement with the value calculated in Ref. [Wysocki et al., 2008]. At temperatures below $T_{c}$, in the case of Fe and Ni, including LSF renders a reduced tiling angle compare to without. Co sticks out from the other systems by showing a rather different behavior where the $\theta_{NN}$ increases gradually at the low temperature regime while rapidly at around $T_{c}$.

Apart from the average tilting angle that reveals only the nearest neighbour ordering, the magnetic order can also be evaluated through the correlation length of the real space correlation function $\langle\mathbf{m}_{0}\cdot\mathbf{m}_{j}\rangle$ fitting to the Ornstein-Zernike expressionOrnstein and Zernike (1914) $\sim exp(-\kappa d)/d$, where $d$ is the intersite distance and $\kappa$ (Å*-1*) the inverse correlation length. A comparison of the correlation length may not be elucidated directly from Fig. 9, however $\kappa$ can be analyzed from the fitting of the Ornstein-Zernike expression. In Tab. 3 the summarized inverse correlation length $\kappa$ obtained for $T/T_{c}=1.25$ are listed. Ni is found having the longest correlation length in comparison with the other systems. Compared with the previous theoretical workRosengaard and Johansson (1997), we obtained a good agreement in general. Despite we have taken into account both contributions from the nearest- and the next nearest-shells, whereas it did not in Ref. [Rosengaard and Johansson, 1997], the $\kappa$ is proved to be insensitive to this aspect. The surprisingly good agreement with the experimental value however, should not be over-credited. Due to small finite size effects present in the simulation data, we included a constant term in the least-square fitting procedure that shifts the entire function. In addition, while the monotonous decease of the Ornstein-Zernike expression is suitable for describing the long-range tail of the correlation function, it is unable to describe the oscillation noises due to the underlying crystal structure. The cumulative effects from those factors give rise to an estimated error bar of about $\pm 10\%$.

III.7 Remarks of treatment of exchange, phase space measure and connection to previous works

In this section, we will discuss in more details about some technical challenges and design choices in our implementation and comparison with previous works. From the outset, our implementation was designed to treat arbitrary number of magnetic components which give us some freedom for single-component systems such as Fe, Co and Ni. These systems could in our implementation be treated as binary alloys which could slightly affect the energetics of the exchange and LSF terms within each trial update in the Monte Carlo simulations. As mentioned in Sec. II.2, the parameters in the simulations depend on the magnitude of the trial moment as well as the average magnitude of the surrounding moments of each of the other components. If there is only a single component, then it reduces to the local moment magnitude and the LSF term is nothing else than the Landau-type expression that was employed, for example, in Ref.[Rosengaard and Johansson, 1997] for Fe,Co and Ni, [Ma and Dudarev, 2012] for Fe and in [Derlet, 2012] for FeRh. Using the rescaled exchange parameters from LDA, we repeat calculations of $T_{c}$ for Fe and Ni but as treated as single component. We obtained 935 K and 308 K, respectively, compared to 926 K and 432 K as from Tab. 1 treated as binary. Fe is basically unaffected by the simplified treatment due to the stable moments in which the LSF energy has a deep minimum, whereas Ni is found to be even softer than using the more advanced binary component treatment. If instead the “bare” exchange interactions from the DLM state without rescaling are employed, the values of $T_{c}$ changes to 1120 K and 617 K for Fe and Ni, respectively. It is worth noting that the total energy in the FM limit is however not correct in this case. The value of Fe then becomes slightly overestimated compared to experiments (1043 K) while for Ni, excellent results that are remarkable close to the experimental value of 633 K is found. This rather good agreement can only be considered as fortuitous where shortcomings of the energetics from the underlying LDA (or GGA) calculations are balancing out each other. The value of Ni is also remarkable close to what was found in Ref.[Ruban et al., 2007] were similar exchange interactions from the DLM state was employed, but in that study a more elaborate LSF energy term was employed which indeed may improve the description of the high temperature properties.

As pointed out in Section. II.3, the choice of PSM does matter in models including longitudinal spin fluctuations. Using the PSM including Jacobian weight instead of uniform put larger statistical weight to states with larger moments and in general that causes larger resistance to decreasing the moment. As a consequence, for all systems, apart from Fe that shows minor differences due to its stable moment, the critical temperature are considerably higher with the Jacobian included. For instance, the $T_{c}$ of Ni (LDA) is increased from 432 K to 613 K and in Py from 482 K to 708 K bringing the values much closer to experiment (633 K and 850 K, respectively). A more detailed comparison of different choices of PSM in the simulations is however left for a future study.

IV Summary

We have constructed a general framework for atomistic simulations of multi-component random alloys including longitudinal fluctuations. The main ingredients consist of mapping total energies from electronic structure calculations for a number of fixed magnetic moments in both the low temperature ferromagnetic limit and the high temperature paramagnetic limit. The model is then “exact” in the both limits and interpolating between the two at intermediate temperatures. Within the model, each magnetic moment magnitude is allowed to change in a fashion that approximately corresponds to Stoner excitations in real materials and the model is constructed in such a way that the strength of the exchange energy is independent on the starting reference state of the mapping, which is typically not the case found in previous studies.

The computational framework has been implemented in Monte Carlo simulations and applied to the elemental transition metals Fe,Co and Ni, together with the binary random alloys Py and Fe-Co alloys. Regarding Curie temperatures, with the exception of Fe, the calculated values in general are slightly underestimated compared to experiments. However, the values are rather insensitive from the starting reference state which is an attractive feature lacking in “normal” simulations for most materials. The simulations do qualitatively describe rather well the high temperature magnetic properties such as the correlation where excellent results compared to experiments are found.

Acknowledgements.

The authors acknowledge valuable discussions with Levente Vitos. The work was financed through the VR (Swedish Research Council), SeRC (Swedish e-Science Research Centre) and GGS (Göran Gustafssons Stiftelser). A.B. acknowledges support from eSSENCE and the CEA-Enhanced Eurotalents programme, co-funded by FP7 Marie Skłodowska-Curie COFUND Programme (Grant Agreement n° 600382). The computations were performed on resources provided by SNIC (Swedish National Infrastructure for computing) at NSC (National Supercomputer Centre) in Linköping, Sweden.