Finite temperature elastic constants of paramagnetic materials within the disordered local moment picture from ab initio molecular dynamics calculations

# Finite temperature elastic constants of paramagnetic materials within the disordered local moment picture from ab initio molecular dynamics calculations

## Abstract

We present a theoretical scheme to calculate the elastic constants of magnetic materials in the high-temperature paramagnetic state. Our approach is based on a combination of disordered local moments picture and ab initio molecular dynamics (DLM-MD). Moreover, we investigate a possibility to enhance the efficiency of the simulations of elastic properties using recently introduced method: symmetry imposed force constant temperature dependent effective potential (SIFC-TDEP). We have chosen cubic paramagnetic CrN as a model system. This is done due to its technological importance and its demonstrated strong coupling between magnetic and lattice degrees of freedom. We have studied the temperature dependent single-crystal and polycrystalline elastic constants of paramagentic CrN up to 1200 K. The obtained results at T= 300 K agree well with the experimental values of polycrystalline elastic constants as well as Poisson ratio at room temperature. We observe that the Young’s modulus is strongly dependent on temperature, decreasing by 14% from T=300 K to 1200 K. In addition we have studied the elastic anisotropy of CrN as a function of temperature and we observe that CrN becomes substantially more isotropic as the temperature increases. We demonstrate that the use of Birch law may lead to substantial errors for calculations of temperature induced changes of elastic moduli. The proposed methodology can be used for accurate predictions of mechanical properties of magnetic materials at temperatures above their magnetic order-disorder phase transition.

###### pacs:
71.15.Pd, 65.40.-b, 62.20.de

## I Introduction

Elastic properties, an important part of the mechanical response of a material, are among the major properties to be studied in theorteical simulations. For instance, Barannikova et al.Barannikova et al. (2012) have recently shown that there is a significant correlation between the elastic and plastic processes that are involved simultaneously in deforming alloys. It is known that for magnetic materials the existance of local magnetic moments above the magnetic transition temperature, in the paramgentic state noticeably affects the elastic propertiesRivadulla et al. (2009); Alling et al. (2010a). Thus, a possibility to predict elastic moduli of magnetic materials in their high temperature paramagnetic state as a function of temperature is highly requested.
The main approach to incorporate temperature in theoretical studies of elastic properties of magnetic materials have been through approximations made in ab initio schemes by excluding the implicit effect of lattice vibrationsLeonov et al. (2012), or by including the thermal expansion effects, often using the experimental dataZhang et al. (2011); Ruban and Razumovskiy (2012). While lattice expansion is believed to be the most important contribution to the temperature dependence of elastic moduli, it would be worthwhile to develop methods which enable us to directly investigate the full effect of the temperature. Unfortunately, a consistent description of a paramagnetic state of a magnetic material is a highly non-trivial theoretical taskAbrikosov et al. (2016). In particular, we are not aware of any theoretical calculation of elastic moduli, where lattice vibrations and magnetic disorder are included on the same footing. However, this could be achieved by the disordered local moments molecular dynamics (DLM-MD)Steneteg et al. (2012); Abrikosov et al. (2016).
In this article, we have used DLM-MD method to study the temperature-dependent elastic moduli of paramagnetic B1 CrN chosen as our model system to demonstrate the functionality of our method. CrN, is a very interesting system considering both its industrial applications in hard coatings and its physical properties. Many experimental and theoretical groups have studied CrN because of its wide range of applicationsCorliss et al. (1960); Gall et al. (2002); Mayrhofer et al. (2008); Rivadulla et al. (2009); Wang et al. (2012); Alling et al. (2010b); Steneteg et al. (2012); Zhou et al. (2014). Apart from its valuable applications, CrN exhibits some fascinating fundamental physical properties. Just below room temperature, , CrN undergoes a phase transition from an orthorhombic antiferromagnetic (AFM) phase to a cubic B1 paramagnetic (PM) phaseCorliss et al. (1960); Rivadulla et al. (2009); Wang et al. (2012). The phase transition in CrN is due to magnetic entropy and the structural change is related to magnetic stress. The volume of the unit cell reduces by when CrN transforms from cubic to orthorhombicCorliss et al. (1960); Browne et al. (1970). Shulumba et al have used a combination of DLM-MD and TDEP to study vibrational free energy and phase stability of CrNShulumba et al. (2014). In their study the magnetic and vibrational degrees of freedom are coupled through the forces from DLM-MD. They found the transition temperature of CrN to be around 380 K which is a closer value to the experimental value of 286 KRivadulla et al. (2009) than the result obtained neglecting the effect of lattice vibrations.
Among all the research that has been done on CrN, only a few has studied its elastic propertiesAlling et al. (2010a, b); Steneteg et al. (2012); Wang et al. (2012); Zhou et al. (2014); Wang et al. (2016). Importantly, theoretical studies of the elastic properties of CrN, up to this date, have been carried out only at Zhou et al. (2014). An important explanation for this is that until recently there were no appropriate computational tools to treat simultaneously effects of lattice vibrations and magnetic disorder. Thus, application of the proposed technique to study temperature dependent elastic properties of CrN is justified. In addition, we investigate the possibility to enhance the efficiency of the simulations of elastic properties and present the results obtained from the recently developed method, symmetry imposed force constants temperature dependent effective potential (SIFC-TDEP)Shulumba et al. (2015), combined with the DLM picture, by relating the phonon properties of the magnetic disorder system and elasticity. Moreover, we compare both theoretical schemes with the results from more conventional calculations based on the Birch law in which the effect of temperature on the elastic moduli is introduced through the thermal expansion and discuss the accuracy of all the schemes considered in this study.

## Ii Methodology

### ii.1 Elastic Properties

To simulate the paramagnetic phase we have used the disordered local moments molecular dynamics (DLM-MD)Steneteg et al. (2012). In this method, the local moments are spatially disordered and the magnetic state of the system is modified periodically and rearranged randomly with a specific time step, spin flip time, during the course of MD simulation. Using the experimental volumes at specific temperatures, we then apply five different deformations to the lattice and perform separate DLM-MD calculations for each deformation. More details on the DLM-MD method are given in Sec. II.2.
The stress-strain relation in the Voigt notation is defined asGrimvall (1999)

 σi=∑jCijϵj (1)

where is the stress tensor and s are the elements of the elastic tensor space. In a cubic system, due to the symmetry, the elastic tensor will only have three non-vanishing, unidentical elements , and . To obtain the elastic constants, we have used the following deformation matrixSteneteg et al. (2013)

 ϵη=⎛⎜⎝1+ηη/20η/210001⎞⎟⎠ (2)

Inserting this matrix into eq. 1, the elastic constants are derived as

 dσ1(T)dη=C11(T) (3)
 dσ2(T)dη=C12(T) (4)
 dσ6(T)dη=C44(T) (5)

First we calculate the stress, , for a set of deformations, , with deviating just a few percent, in our calculations 1%, from zero. Then, we obtain the numerical derivative of and extract the elastic constants. For each temperature and the volume at that temperature, we calculate stresses, , for a set of molecular dynamics time steps, . Thus, for each , we will have number of stresses. The derivative is numerically calculated by fitting a line to these points using the least square method.
In principle, single crystal samples are often not available, thus the measurements of individual elastic constants, , are rare. In many cases, polycrystalline materials are studied experimentally for which one may determine the polycrystalline bulk modulus , Young’s modulus and shear modulus . Using DLM-MD theory, we can calculate the elastic properties of a single crystal but by using Voigt and Reuss approaches, we can obtain expressions for bulk and shear moduli in polycrystals. For a cubic system these properties are derived as

 BV=BR=B (6)
 B=C11+2C123 (7)
 GV=C11−C12+3C445 (8)
 GR=5(C11−C12)C443(C11−C12)+4C44 (9)

The Young modulus and the Poisson ratio can also be calculated according to the following relations.

 EV,R=9BGV,R3B+GV,R (10)
 νV,R=3B−2GV,R2(3B+GV,R) (11)

It is also useful to define the elastic anisotropy for polycrystalline materials.

 AV,R=GV−GRGV+GR (12)

The Voingt and the Reuss averaging of elastic constants, Eq. 6-12, will give us the upper and the lower bound of elastic constants, respectively. the Voigt-Reuss-Hill approach (the Hill approximation) combines these two limits by averaging over the Voigt and the Reuss elastic constants, assuming that this average gives a good approximation for the actual macroscopic elastic constantsHirsekorn (1990); Toonder et al. (2000).

 EH=EV+ER2 (13)
 GH=GV+GR2 (14)

and

 νH=EH2GH−1 (15)

### ii.2 Details of DLM-MD Simulations

The DLM-MD method was suggested by Steneteg et al. to simulate the paramagnetic state of magnetic materials at finite temperatures including lattice vibrationsSteneteg et al. (2012). In this method, the disordered local moments (DLM) picture of paramagnetism is implemented in the framework of the ab initio molecular dynamics (MD). The DLM approach is introduced by HubbardHubbard (1981, 1979a, 1979b) and HasegawaHasegawa (1979, 1980) and later on applied by Gyorffy et al. within the coherent potential approximation (CPA) electronic structure frameworkGyorffy et al. (1985). Beside the LDA+U approximation for the exchange correlation energy and the supercell description of the disordered magnetic state, additional approximations in DLM-MD approach need to be noted. First, we consider the local moments to have collinear orientations. Indeed it is well justified to ignore the noncollinear orientations of local moments in paramagnetic materials in a description of thermodynamic potentials, well above the magnetic transition temperatureGyorffy et al. (1985). The fast longitudinal fluctuations of magnetic moments dictating their magnitude are treated as fast electronic degrees of freedom governed by the self-consistent solution of the electronic structure problem at each MD step.
Following the same arguments as in Ref. dccclxxx(28), we assume that the magnetic configuration of the system, even though ergodic, does not cover its phase space uniformly in time but rather gets stuck for a specified time (we denote as spin-flip time, ) near the points with a specific configuration of moments at every site pointing in a random direction and then moves instantly to another disordered but different state of the phase space. In fact, Gyorffy et al. supposed that the changes in the orientational configuration of the moments characterize the motion of temporarily broken ergodicity to a large degree.
It should be noted that on one hand, in the paramagnetic state at high temperatures the magnetic excitations exhibit different characteristics from those of low temperature. As an important example, the relevant time scale of its dynamics can be better estimated by the spin decoherence time, , of individual or pair of moments, rather than the inverse spin-wave frequency related to the collective motion of many spins. This means that the relevant magnetic dynamics are much faster in magnetically disordered materials (of the order of - s) in comparison to that of in magnetically ordered systems (of the order of s).
A further complication is that the paramagnetic state is a high temperature state of a magnetic material for which considering the atomic motions and displcements from ideal lattice sites in the simulations is essentialAlling et al. (2016). The time scale for collective atomic motions constituting lattice vibrations is estimated by the inverse of the Debye frequency (). However, in molecular dynamics simulations, we need to calculate the displacements and forces acting on individual atoms on a much shorter time scale of the order of s for a proper description of Born-Oppenheimer (B-O) dynamics.
Figure. 1 illustrates schematically the basic idea of the DLM-MD technique and demonstrates all the different timescales of the method including the B-O MD time step (), the spin-flip time () as well as a time indicator for the “phonon timescales” corresponding to the period time of the highest frequency phonons. The values we have adopted in our simulations, as shown in the figure, comply with the Born-Oppenheimer approximation indicating that electrons exhibit the fastest timescale as compared to other degrees of freedom. Then comes the timescale of the individual atomic motions ( fs). The change in the temporarily broken ergodicity magnetic state of the disordered local magnetic moments comes next ( fs) and only then comes the timescale of the full collective atomic vibrational modes ( fs).
Strictly speaking, the appropriate value of the spin flip time should either be taken from experiments or calculated from real spin dynamics simulations. However, in Ref. dccclxxx(8) it was shown that should be chosen short enough to assure an adiabatic approximation. In particular, several tests with different spin flip times carried out in Ref. dccclxxx(8) demonstrate that fs gives reliable description for CrN.

We run DLM-MD simulations with total number of 5000 steps, resulting in a length of 5 ps for one simulation. The schematic representation of this algorithm is shown is Fig. 1.
Using DLM-MD, the elastic constants for PM phase is calculated at five different temperatures 300, 600, 800, 1000, 1200. For each temperature, five different values of distortions, have been used for the deformation matrix . We can then extract the values, using Eq. 3-5 and calculate the elastic constants at each temperature.
We note that for DLM-MD calculations of constants, we do not need to use the projected cubic elastic constants as described in Sec. II.3, as we sample the phase space of possible magnetic configurations and average out non-cubic magnetic symmetry on the fly.

### ii.3 Finite Temperature Elastic Constants From Static Calculations Including Thermal Expansion Effects

In order to simulate the paramagnetic state of CrN in a static lattice approximation at T=0 K, we have used the DLM picture combined with magnetic special quasirandom structure (SQS) approachAlling et al. (2010b). In our static DLM-SQS calculations, we have employed a cubic supercell with 108 Cr and 108 N atoms in which the spin-up and spin-down Cr moments are mimicing a random alloy distributionZunger et al. (1990).
To obtain the elastic properties, we apply a set of different distortions to this supercell. As explained in Sec. II.1, after performing the first-principles calculations for each of these supercells, we calculate the stress. Thereafter, we obtain the derivative of the stress numerically. This derivative will provide us with different values.
However, one should bear in mind that due to magnetic disorder the cubic symmetry of the supercell is broken and all the elements of the elastic matrix, , , , , , , , and are not identical and need to be calculated. Then we use the projection technique to determine average cubic elastic constantsTasnádi et al. (2012) of the simulated cubic crystal via

 ¯C11=C11+C22+C333 (16)
 ¯C12=C12+C13+C233 (17)
 ¯C44=C44+C55+C663 (18)

Thermal expansion is one important manifestation of anharmonicity. Statically, one can include the temperature effect on the elastic properties through the thermal expansionShang et al. (2010); Wang et al. (2010). This approach is based on Birch’s lawBirch (1961) that includes temperature effects in elastic properties implicitly via thermal expansionLedbetter (2006). Thus within this approach

 Cij(T)≈C0ij(V(T)) (19)

where , are the elastic constants calculated at zero temperature but at volume corresponding to simulation temperature . Note that, for these set of calculations, we should use the average elastic constants as stated in Eq. 16-18 since we are considering a magnetically disordered state. In this work, we have used the experimental thermal expansion obtained from Ref. dccclxxx(12).

### ii.4 Finite Temperature Elastic Constants From SIFC-TDEP

We concentrate on mechanical properties calculated for single crystal materials from first principles simulations. The elastic constants can be defined either by Hooke’s law, described above or in long wave limit of phonons. In this section we give a brief description of the calculation of the elastic properties through the phonons and introducing the methods that include explicit temperature dependence of phonon properties and therefore elastic constants.

 ^H=U0+∑iαp2iα2mi+12!∑ij∑αβΦijαβuiαujβ, (20)

where and are indices for the Cartesian coordinates and and are the Cartesian components of the displacement of the atom and , and is the potential energy of the static lattice; and are the momentum and mass of atom . The are the interatomic force constants, which express a relation between the force in direction acting on the atom , when atom is displaced in direction .

In long wave limit , the atoms move very slowly with vanishing frequencies and these sound waves have frequencies that are determined by macroscopic elastic constants. We can write the relation between elastic and force constantsLeibfried and Ludwig (1961); Fultz (2010); Maradudin and Fein (1962)

 ~Cijkl=−12V∑nΦij0nrknrln, (21)

where is a full elastic constants matrix, that could be later reduced due to the symmetry and transformed to the Voigt notation. Here the interatomic force constants (IFC) between atoms and in the unit cell ; and are the positions of atoms and in the unit cell n taken in respect to the reference unit cell 0. To calculate interatomic force constants we use the temperature dependent effective potential (TDEP) method Hellman et al. (2013, 2011). The advantage of this method is that the temperature dependence of the IFC is included explicitly. We determine the force constants by minimising the difference in forces from the model Hamiltonian (Eq. (20)) and real system simulated with the DLM-MD. The DLM-MD provides us with a set of displacements and forces at each step, so we minimize the following quantity

 min¯¯ΦΔF=1NtNt∑t=1∣∣F% DLM-MDt−FHt∣∣2==1NtNt∑t=1∣∣F% DLM-MDt−¯¯Φ(UDLM-MDt)∣∣2, (22)

The solution for the force constants is obtained by linear square method that minimises . The symmetry constraints on the force constants matrices are applied Hellman et al. (2013, 2011) and this procedure reduces the computational cost. Using TDEP one can get temperaturevolume dependent IFC. TDEP works for the ordered structures. In disordered systems the symmetry of the crystal is broken. Therefore, we use a generalization of the TDEP method towards disordered systems, the so-called SIFC-TDEPShulumba et al. (2015). Here it is applied to calculate vibrational properties at finite temperatures for magnetically disordered systems and to evaluate their elastic properties. Using the SIFC-TDEP method we extract the effective interatomic force constants () by treating Cr atoms as symmetry equivalent and imposing the full symmetry of the underlying crystal lattice on the IFCs. Obtained effective IFCs are then used to calculate elastic constants (Eq. (21)).

It is well known that there is an issue of using Eq. (21) to calculate accurately absolute values of the finite temperature elastic constants, because the result depends on the finite size of the simulation cells. On the other hand, it has been demonstrated that using the SIFC-TDEP and the original TDEP we are able to accurately determine the temperature dependence of phonon frequencies of CrNShulumba et al. (2014) and other materialsLi et al. (2014); Romero et al. (2015).

Therefore, we employ the finite temperature scaling of elastic constants to obtain the final expression, given by

 Cij(V,T)=Cstatij(V,T0)Cphij(V,T)Cphij(V,T0) (23)

with K in our calculations. Note that one can avoid the scaling of the elastic constants by performing calculations on a larger cell. Thus elastic constants from SIFC-TDEP can be calculated including both volume and temperature dependence of IFC, extracted from SIFC-TDEP.

SIFC-TDEP scheme allows one to calculate the finite temperature elastic constants at a fraction (1/5) of computational cost of DLM-MD, in which the supercell has to be distorted using deformation matrix, Eq. 2. Shulumba et al.Shulumba et al. (2015) showed that 80% of the temperature effect on the elastic constant of TiN could be captured with this method, at 20% of the computational cost. In this work, we extend the approach towards the magnetically disordered systems.

### ii.5 Computational Details

All the static DLM and DLM-MD calculations are carried out within the projected augmented wave method (PAW)Blöchl (1994) as implemented in Vienna Ab-initio Simulation Package (VASP)Kresse and Hafner (1993); Kresse and Furthmüller (1996a, b); Kresse and Joubert (1999). For the electronic exchange-correlation effects we have used a combination of local density approximation with a Hubbard Coloumb term (LDA+U)Anisimov et al. (1991). The effective Hubbard term value, , is chosen to be 3 eV for Cr orbitals which is shown to be the optimal value obtained from a thorough theoretical comparison of the structural and electronic properties of CrN with experimental measurementsAlling et al. (2010b). The energy cutoff is set to 500 eV.
We have used a supercell consisting of repetitions of the conventional cubic cell including 8 atoms, giving in total 108 Cr and 108 N atoms. In the PM phase, the spin-up and spin-down magnetic moments are randomly distributed on Cr atoms. The Brillouin zone is sampled using a Monkhorst-Pack schemePack and Monkhorst (1977) with a -mesh of . In order to maintain the desired temperature in our MD calculations, we have used the canonical ensemble (NVT). We have used the Nose thermostatNosé (1991) with the default mass value as it is implemented in VASP in our simulations.
The thermal expansion is included in our calculations using the experimental lattice constants as a function of temperatureWang et al. (2012).

For SIFC-TDEP calculations we run DLM-MD with the same setting as described above with a difference in supercell size. We used a supercell contained 32 chromium and 32 nitrogen atoms arranged in 222 conventional unit cells. We ran the simulations on a grid of 6 temperatures and six volumes the NVT ensemble. From these simulations we extract IFC as a function of volume and temperature. The effective IFC were found to be smooth and easily interpolated across the whole temperature/volume interval. We calculate the theoretical thermal expansion, which is in good agreement with the experimental valuesWang et al. (2012). Elastic constants from SIFC-TDEP are evaluated for the theoretical thermal expansion.
When dealing with MD simulations, we should make sure that the obtained results are well converged, i.e. that the statistical errors are small. The output of an MD run is reported in terms of of the time average, in our case of the elastic constants. Since the simulation times are of finite size, an statistical imprecision of this average values is expected. Our simulation time is 5 ps and in order to estimate the uncertainty of our MD results, we have used a t distribution with 95% confidence interval taking the correlated nature of each MD time step into account according to the method suggested by Allen and TildesleyAllen and Tildesley (1991). In our case we find that the factor of uncorrelated time steps correspond to between 1 to 200 configurations, which add completely new information to the average values, depending on the temperature. The 95% confidence interval for the mean values of elastic constants and derived properties are given as error bars calculated in this way.

## Iii Results

### iii.1 Single crystal Elastic Constants

Figure 2, shows the obtained temperature dependent elastic constants of the paramagnetic B1 CrN from different methods. We can see that the low-temperature limit of the elastic constants calculated at finite temperatures via DLM-MD method are in good agreement with the corresponding elastic constants obtained from static zero-temperature calculations. The zero-Kelvin values also agree well with earlier theoretical calculations (see, for instance, Ref.Zhou et al. (2014)).

In Tab. 1 we report the calculated elastic properties of nonmagnetic (NM), antiferromagnetic (AFM) and paramagnetic (PM) B1-CrN at zero kelvin. The values of the elastic constants from our AFM and PM calculations show that AFM calculations in general tend to overestimate the elastic constants. Relatively small difference between AFM and PM calculations can be attributed to relatively weak magnetic exchange interactions in B1 CrNLindmaa et al. (2013). One should expect substantially larger effect in systems with higher Neél or Curie temperature. Simultaneously, from the values listed in Tab. 1 we can see that the nonmagnetic calculations give a negative value for . This implies that NM B1-CrN is mechanically unstable as the values does not satisfy the Born-Huang stability criterionBorn and Huang (1988) for which the shear elastic coefficient should be positive. Most other elastic moduli come out very wrong as well in non-magnetic calculations. This indicates the importance of magnetic effects which need to be considered while simulating the PM state of a magnetic material.
Tab. 1 also summarizes the elastic properties of PM B1 CrN at T=300 K. We observe that our calculations overestimate by 51 GPa, 9%, as compared to the experimental valueAlmer et al. (2003); Chen et al. (2004) which is reasonable given LDA+U normal uncertainty. The difference between the theoretical and experimental values for is 53 GPa, 38%, and 75 GPa, 74% for . However, we do not trust the experimental values of and and emphasize that our results are in good agreement with other first-principles calculationsZhou et al. (2014). Moreover, our average Young’s modulus value at 300 K, E=433 GPa is in good agreement with the reported experimental value of 400 GPaHolleck (1986). Considering this excellent agreement between theory and experiment for E, additional independent measurements of the values at least at room temperature are desired.
Note that the numerical accuracy of elastic constants calculated with DLM-MD method is quite high. The statistical errors for and are at most within 0.1% of the mean values which is very small. Both and decrease almost linearly with increasing temperature. This indicates the normal temperature dependence behavior originating from anharmonicityGrimvall (1999). For , on the other hand we do not see any specific trend and it appears to be nearly temperature independent within the error bars. SIFC-TDEP values with all the elastic constants decrease monotonously.
For many systems in the absence of phase transitions for intermediate temperatures, the temperature dependent elastic constants can be fitted to the empirical relationGrimvall (1999)

 Cij(T)=Cij(0)(1−b(T−TN)) (24)

where is a constant. As we get close to the melting temperature, high-order anharmonic effects result in a strong nonlinear temperature dependenceGrimvall (1999). As for the case of CrN, the , we have , in which T is the temperature range in which we do our simulations and K is the melting temperature of CrN. This argument further justifies the reliability of our method because our simulations result in a nearly linear (within numerical accuracy) temperature dependence of elastic constants.
The finite temperature values from DLM-MD are in good agreement with the data obtained from SIFC-TDEP.

The single crystal elastic constants of PM CrN at 1200 K is given in Table. 1. As can be seen from Fig. 2, DLM-MD calculations at T=300 K give 591 GPa, 102 GPa and 141 GPa for , and , respectively. At T=1200 K, Table. 1, SIFC-TDEP gives a larger value of GPa which differs by 23 GPa, 6% from GPa from DLM-MD. The GPa is smaller by 15 GPa, 12% in SIFC-TDEP as compared to GPa in DLM-MD. The is 132 GPa and 135 GPa from DLM-MD and SIFC-TDEP, respectively. Clearly there is an increasing difference between DLM-MD and SIFC-TDEP values as the temperature increases. On the other hand, SIFC-TDEP clearly shows much higher numerical stability and much smoother temperature dependence as compared to the DLM-MD calculations, which shows the usefulness of this numerically efficient technique for calculations of the temperature elastic constants in a not too broad temperature interval.
The red triangles in Fig. 2, show the results from DLM static calculations including thermal expansion. The method gives close values of in comparison to DLM-MD and SIFC-TDEP , but the data for other elastic constants, and differ stronger. This divergence demonstrates that the Birch law may be violated in real systems at finite temperatures. This implies that incorporating the temperature effect through thermal expansion may not be an accurate way to obtain the temperature dependence of elastic properties of magnetic materials in their paramagnetic state. In order to get a good picture for finite temperature elastic properties in PM CrN, we need to treat lattice vibrations, magnetic configuration and the effect of thermal expansion on the same footing as it is done in our DLM-MD simulations.

### iii.2 Polycrystalline Elastic Constants

Using single crystal elastic constants, obtained from our methods, we can calculate the temperature dependent polycrystalline elastic constants and the Poisson ratio for PM B1 CrN. The results are displayed in Fig. 3. Similar to what we see in Fig. 2 for , B, E and G moduli show nearly linear temperature dependence following eq. 24. Poisson ratio shows a little bit of variation as the temperature increases. As stated for the single crystal elastic constants, we see a good agreement between the room temperature values obtained from our DLM-MD calculations and from conventional zero kelvin static calculations. The polycrystalline elastic constants values obtained from DLM-MD and SIFC-TDEP are close to each other suggesting that both methods give similar results.
At zero kelvin, the bulk modulus (B) from AFM calculations, as listed in Tab. 1, is found to be 292 GPa which is significantly larger than the PM value of 273 GPa. There is even a larger difference between the values of bulk moduli obtained from PM and NM calculations. The nonmagnetic calculations give the bulk modulus of 387 GPa. However, it is shown in other studies that simulating the paramagnetic state as nonmagnetic could lead to erroneous conclusions. Rivadulla et al. carried out calculations on CrN considering the PM state as nonmagnetic. Their study gave a bulk modulus value of 340 GPa for the cubic CrN showing a drastic reduction of about 25% as compared to that of the orthorhombic phase (255 GPa)Rivadulla et al. (2009). Alling et al.Alling et al. (2010a) calculated the bulk modulus of cubic CrN using DLM to simulate the paramagnetic state and they found the bulk modulus to be 252 GPa which is very close to the bulk modulus of the orthorhombic phase. These results highlight the importance of taking into account the finite local moments during the simulation of the paramagnetic state.
There are several experimental studies on the Young’s modulus of CrN measured for thin films Holleck (1986); Attar and Johannesson (1995) with preferred orientationsAttar and Johannesson (1995); Cunha et al. (1999); Chen et al. (2004). The reported experimental values for CrN range from 324 GPa to 461 GPa. Our obtained values from room temperature calculations are summarized in Tab. 1. We observe that our results from both DLM-MD and SIFC-TDEP, are well within the range of reported experimental values. Our average static Poisson ratio value, , is in good agreement with the values derived from experimental data, 0.28Fabis (1990); Cunha and Andritschky (1999); Qiu et al. (2013) and 0.24Chen et al. (2004). From Fig. 3, we can observe that all polycrystalline elastic constants have a fairly strong temperature dependence, decreasing by almost 14% in DLM-MD and 8% in SIFC-TDEP between room temperature and 1200 .

### iii.3 Elastic Anisotropy

It is possible to determine the measure of the elastic anisotropy experimentally by the strain ratio. The strain ratio can also be recalculated by the ratio between the Young’s moduli in different directionsTasnadi et al. (2010). We can calculate the from the single-crystal elastic constantsGrimvall (1999). Our DLM-MD simulations at room temperature give GPa, GPa and GPa. Even though, there is a wide discrepancy in the experimentally measured directional Young’s modulus valuesChen et al. (2004), our data are fairly consistent with the experimental values 290 GPaCunha et al. (1999), 520 GPaAttar and Johannesson (1995) and 30020 GPaChen et al. (2004) for , and directions, respectively.
In order to quantify the temperature dependence of elastic anisotropy, we calculated the anisotropy according to Voigt-Reuss-Hill definition, Eq. 12, as well as according to ZenerGrimvall (1999), Fig. 4.

 AZ=2C44C11−C12 (25)

If the material is isotropic, the former is equal to 0 as the temperature increases and the latter is equal to 1. In Fig. 4 we see that CrN becomes more isotropic at higher temperatures.

## Iv Summary and Conclusion

We have used first-principles simulations based on ab initio molecular dynamics (AIMD) in combination with disordered local moments (DLM) method to study the finite temperature elastic properties of magnetic material CrN, in its high-T paramagnetic state. Though these simulations are computationally expensive and are substantially more time consuming as compared to conventional static calculations, we see that performing MD to obtain finite temperature elastic constants is needed to get a proper description of the elastic constants at elevated temperatures. Moreover, we have used the recently developed method, SIFC-TDEP, to study temperature dependent elastic constants of CrN. This method also uses DLM-MD but allows one to calculate elastic constants without simulations at distorted lattices. Thus it has higher computational efficiency. In general, we see that both DLM-MD and SIFC-TDEP give results that are in good agreement with each other and with available experiment. On the other hand, the use of Birch law may give larger errors for calculated elastic constants.
We study the temperature dependent elastic properties of prototypical paramagnetic transition metal nitride, CrN, between room temperature and 1200 which corresponds to operation temperature of cutting tools. We have calculated the single crystal elastic constants of PM cubic CrN, , and , as well as its polycrystalline elastic constants, B, G and E being the bulk, shear and Young’s moduli, respectively and also the Poisson ratio, . We observe that the elastic constants decrease nearly linearly with increasing temperature which is the predicted temperature dependent behavior, caused by anharmonicity. We see that polycrystalline elastic constants decrease by between room temperature and 1200 . Studying the elastic anisotropy, demonstrates that the material becomes substantially more isotropic at elevated temperatures. Therefore, the effect of temperature on elastic properties is strong and should be included in the studies of materials functioning at high-T environments. The proposed technique allows for a reliable inclusion of finite temperature effects in ab initio simulations of elastic properties of magnetic materials.

## Acknowledgments

E.M. would also like to asserts her appreciations to Dr. Ferenc Tasnadi for the very useful discussions. Financial support from the Swedish Research Council (VR) through Grant No. 621-2011-4426 and the Swedish Foundation for Strategic Research (SSF) program SRL Grant No. 10-0026 is Acknowledged. I.A.A. is grateful for the Grant from the Ministry of Education and Science of the Russian Federation (Grant No. 14.Y26.31.0005) and Tomsk State University Academic D. I. Mendeleev Fund Program. B. A. acknowledges financial support by the Swedish Research Council (VR) through the young researcher grant No. 621-2011-4417 and the international career grant No. 330-2014-6336 and Marie Sklodowska Curie Actions, Cofund, Project INCA 600398. Moreover, we acknowledge support from the Swedish Government Strategic Research Area in Materials Science on Functional Materials at Linkoöping University (Faculty Grant SFOMatLiU No 2009 00971). The authors would like to acknowledge the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Center (NSC) for providing the computational resources.

### References

1. S. A. Barannikova, A. V. Ponomareva, L. B. Zuev, Y. K. Vekilov,  and I. A. Abrikosov, Solid State Commun. 152, 784 (2012).
2. F. Rivadulla, M. Bañobre-López, C. X. Quintela, A. Piñeiro, V. Pardo, D. Baldomir, M. A. López-Quintela, J. Rivas, C. A. Ramos, H. Salva, J.-S. Zhou,  and J. B. Goodenough, Nat. Mater. 8, 947 (2009).
3. B. Alling, T. Marten,  and I. A. Abrikosov, Nat. Mater. 9, 283 (2010a).
4. I. Leonov, A. I. Poteryaev, V. I. Anisimov,  and D. Vollhardt, Phys. Rev. B 85, 020401 (2012)arXiv:arXiv:1110.0439v1 .
5. H. Zhang, B. Johansson,  and L. Vitos, Phys. Rev. B 84, 140411 (2011).
6. A. V. Ruban and V. I. Razumovskiy, Phys. Rev. B 85, 174407 (2012).
7. I. A. Abrikosov, A. V. Ponomareva, P. Steneteg, S. A. Barannikova,  and B. Alling, Curr. Opin. Solid State Mater. Sci. 20, 85 (2016).
8. P. Steneteg, B. Alling,  and I. A. Abrikosov, Phys. Rev. B 85, 144404 (2012).
9. L. M. Corliss, N. Elliott,  and J. M. Hastings, Phys. Rev. 117, 929 (1960).
10. D. Gall, C.-S. Shin, R. T. Haasch, I. Petrov,  and J. E. Greene, J. Appl. Phys. 91, 5882 (2002).
11. P. H. Mayrhofer, H. Willmann, L. Hultman,  and C. Mitterer, J. Phys. D. Appl. Phys. 41, 155316 (2008).
12. S. Wang, X. Yu, J. Zhang, M. Chen, J. Zhu, L. Wang, D. He, Z. Lin, R. Zhang, K. Leinenweber,  and Y. Zhao, Phys. Rev. B 86, 064111 (2012).
13. B. Alling, T. Marten,  and I. A. Abrikosov, Phys. Rev. B 82, 184430 (2010b).
14. L. Zhou, F. Körmann, D. Holec, M. Bartosik, B. Grabowski, J. Neugebauer,  and P. H. Mayrhofer, Phys. Rev. B 90, 184102 (2014).
15. J. D. Browne, P. R. Liddell, R. Street,  and T. Mills, Phys. Status Solidi 1, 715 (1970).
16. N. Shulumba, B. Alling, O. Hellman, E. Mozafari, P. Steneteg, M. Odén,  and I. A. Abrikosov, Phys. Rev. B 89, 174108 (2014).
17. S. Wang, X. Yu, J. Zhang, L. Wang, K. Leinenweber, D. He,  and Y. Zhao, Cryst. Growth Des. 16, 351 (2016).
18. N. Shulumba, O. Hellman, L. Rogström, Z. Raza, F. Tasnádi, I. A. Abrikosov,  and M. Odén, Appl. Phys. Lett. 107, 231901 (2015).
19. G. Grimvall, Thermophysical properties of materials (Elsevier Science B.V., The Netherlands, 1999).
20. P. Steneteg, O. Hellman, O. Y. Vekilova, N. Shulumba, F. Tasnádi,  and I. A. Abrikosov, Phys. Rev. B 87, 094114 (2013).
21. S. Hirsekorn, Textures Microstruct. 12, 1 (1990).
22. J. M. J. D. Toonder, J. a. W. V. Dommelen,  and F. P. T. Baaijens, Model. Simul. Mater. Sci. Eng. 7, 909 (2000).
23. J. Hubbard, Phys. Rev. B 23, 5974 (1981).
24. J. Hubbard, Phys. Rev. B 20, 4584 (1979a).
25. J. Hubbard, Phys. Rev. B 19, 2626 (1979b).
26. H. Hasegawa, J. Phys. Soc. Japan 46, 1504 (1979).
27. H. Hasegawa, J. Phys. Soc. Japan 49, 178 (1980).
28. B. L. Gyorffy, A. J. Pindor, J. Staunton, G. M. Stocks,  and H. Winter, J. Phys. F Met. Phys. 15, 1337 (1985).
29. B. Alling, F. Körmann, B. Grabowski, A. Glensk, I. A. Abrikosov,  and J. Neugebauer, Press  (2016).
30. A. Zunger, S.-H. Wei, L. G. Ferreira,  and J. E. Bernard, Phys. Rev. Lett. 65, 353 (1990).
31. F. Tasnádi, M. Odén,  and I. A. Abrikosov, Phys. Rev. B 85, 144112 (2012).
32. S.-L. Shang, H. Zhang, Y. Wang,  and Z.-K. Liu, J. Phys. Condens. Matter 22, 375403 (2010).
33. Y. Wang, J. J. Wang, H. Zhang, V. R. Manga, S. L. Shang, L.-Q. Chen,  and Z.-K. Liu, J. Phys. Condens. Matter 22, 225404 (2010).
34. F. Birch, J. Geophys. Res. 66, 2199 (1961).
35. H. Ledbetter, Mater. Sci. Eng. A 442, 31 (2006).
36. G. Leibfried and W. Ludwig, Solid State Phys., Vol. 12 (Academic Press, New York, 1961) pp. 275–444.
37. B. Fultz, Prog. Mater. Sci. 55, 247 (2010).
38. A. A. Maradudin and A. E. Fein, Phys. Rev. 128, 2589 (1962).
39. O. Hellman, P. Steneteg, I. A. Abrikosov,  and S. I. Simak, Phys. Rev. B 87, 104111 (2013).
40. O. Hellman, I. A. Abrikosov,  and S. I. Simak, Phys. Rev. B 84, 180301 (2011).
41. C. W. Li, O. Hellman, J. Ma, A. F. May, H. B. Cao, X. Chen, A. D. Christianson, G. Ehlers, D. J. Singh, B. C. Sales,  and O. Delaire, Phys. Rev. Lett. 112, 175501 (2014).
42. A. H. Romero, E. K. U. Gross, M. J. Verstraete,  and O. Hellman, Phys. Rev. B 91, 214310 (2015).
43. P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
44. G. Kresse and J. Hafner, Phys. Rev. B 48, 13115 (1993).
45. G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996a).
46. G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996b).
47. G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
48. V. I. Anisimov, J. Zaanen,  and O. K. Andersen, Phys. Rev. B 44, 943 (1991).
49. J. D. Pack and H. J. Monkhorst, Phys. Rev. B 16, 1748 (1977).
50. M. P. Allen and D. J. Tildesley, Computer simulation of liquids (Oxford University press, New York, 1991).
51. J. Almer, U. Lienert, R. L. Peng, C. Schlauer,  and M. Odén, J. Appl. Phys. 94, 697 (2003).
52. H. Y. Chen, C. J. Tsai,  and F. H. Lu, Surf. Coatings Technol. 184, 69 (2004).
53. A. Lindmaa, R. Lizárraga, E. Holmström, I. A. Abrikosov,  and B. Alling, Phys. Rev. B 88, 054414 (2013).
54. M. Born and K. Huang, “Dynamical theory of crystal lattices,”  (1988).
55. F. Attar and T. Johannesson, Thin Solid Films 258, 205 (1995).
56. L. Cunha, M. Andritschky, K. Pischow,  and Z. Wang, Thin Solid Films 355, 465 (1999).
57. L. Cunha and M. Andritschky, Surf. Coatings Technol. 111, 158 (1999).
58. Y. Qiu, S. Zhang, B. Li, Y. Wang, J. W. Lee, F. Li,  and D. Zhao, Surf. Coatings Technol. 231, 357 (2013).
59. F. Tasnadi, I. A. Abrikosov, L. Rogström, J. Almer, M. P. Johansson,  and M. Oden, Appl. Phys. Lett. 97, 231902 (2010).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters