text:
# Chiral exponents in 𝑂(𝑁) × 𝑂(𝑚) spin models at 𝑂(1/𝑁²)

Abstract. The critical exponents corresponding to chirality are computed at $`O(1/N^2)`$ in $`d`$-dimensions at the stable chiral fixed point of a scalar field theory with an $`O(N)`$ $`\times `$ $`O(m)`$ symmetry. Padé-Borel estimates for the exponents are given in three dimensions for the Landau-Ginzburg-Wilson model at $`m`$ $`=`$ $`2`$.

LTH 551

Scalar quantum field theories have been widely used to determine information at criticality for a large range of critical phenomena. The primary tool of computation is the renormalization group equation where the field theory is renormalized in perturbation theory to a certain loop order. The resulting renormalization group functions are then evaluated at the appropriate fixed point, such as that of Wilson and Fisher, , and a numerical estimate provided for the corresponding critical exponents, . A recent, comprehensive and detailed review of the use of the renormalization group in this area is given in . For the much studied Heisenberg system, the $`O(N)`$ $`\varphi ^4`$ theory has been renormalized to five loops in conventional dimensionally regularized perturbation theory,  and references therein, and to six loops in  where the spacetime dimension is fixed to be three. The resulting critical exponents after appropriate resummation are in very reasonable agreement with experiment and other methods such as Monte Carlo simulations and the high temperature series expansion, demonstrating the success of the approach. More recently, there has been renewed interest in understanding the phase transitions in spin models which are similar to the usual Heisenberg spin systems. These are spin systems where there is an additional chiral symmetry which gives rise to fixed points in the system over and above the usual Gaussian free field and Heisenberg fixed points. In other words there is a chiral transition in systems with noncollinear order. Early work on this system was carried out by Kawamura in . From the field theoretic point of view the underlying quantum field theory of such a spin system is a $`\varphi ^4`$ model with an $`O(N)`$ $`\times `$ $`O(m)`$ symmetry where to relate to the chiral property one sets $`m`$ $`=`$ $`2`$ thereby reducing it to a Landau-Ginzburg-Wilson theory. However, the more general models have been studied for arbitrary $`m`$ to try and unravel the true fixed point structure of the field theory. This is important as there appears to be conflicting evidence from theoretical calculations, either using Monte Carlo techniques or renormalization group methods, as to the number of stable fixed points and their existence and as to whether they are first or second order. Prior to  the underlying field theory had only been renormalized to two loops in the $`\overline{\text{MS}}`$ scheme. In  this was extended to three loops and the fixed dimension calculation of critical exponents was performed to six loops in . Moreover, in  the chiral exponents associated with the stable chiral fixed point of the $`O(N)`$ $`\times `$ $`O(m)`$ model have been computed to six loops in three dimensions. They have been shown to be in good agreement with other results, . Briefly the method of calculation is to compute the anomalous dimension of a chiral bilinear operator and evaluate it at the stable chiral fixed point when $`m`$ $`=`$ $`2`$.

Aside from explicit perturbation theory there are other field theory methods with which one can estimate critical exponents. One such method is the large $`N`$ expansion where the exponents are determined order by order in powers of $`1/N`$. The critical point large $`N`$ technique developed in  has proved to be an excellent tool to compute critical exponents. For instance, in the $`O(N)`$ $`\varphi ^4`$ theory the wave function exponent or anomalous dimension, $`\eta `$, is known at $`O(1/N^3)`$ in arbitrary spacetime dimension $`d`$, . Through the critical renormalization group equation the coefficients of the $`ϵ`$-expansion of such exponents are in one-to-one agreement with the coefficients of the corresponding perturbative renormalization group functions evaluated at the same fixed point, where $`d`$ $`=`$ $`4`$ $``$ $`2ϵ`$. Given the recent interest in the $`O(N)`$ $`\times `$ $`O(m)`$ field theory it is the purpose of this article to calculate the same chiral exponents of  in the large $`N`$ expansion at $`O(1/N^2)`$ in $`d`$-dimensions using the approach of . Indeed in  the exponents $`\eta `$ and $`\nu `$ have already been computed at $`O(1/N^2)`$ in $`d`$-dimensions at each fixed point. Our calculation involves determining the critical exponent corresponding to the anomalous dimension of the chiral operator and parallels a similar calculation in the $`O(N)`$ $`\varphi ^4`$ theory itself where the critical exponent of the bilinear operator corresponding to crossover behaviour was determined at $`O(1/N^2)`$, .

The (massless) Lagrangian underlying the chiral phase transition in frustrated spin systems with noncollinear order is based on that for the Landau-Ginzburg-Wilson $`O(N)`$ $`\times `$ $`O(2)`$ model. Therefore, we take

$$L=\frac{1}{2}^\mu \varphi ^{a\alpha }_\mu \varphi ^{a\alpha }+\frac{\overline{g}}{4!}\left(\varphi ^{a\alpha }\varphi ^{a\alpha }\right)^2+\frac{g_T}{4!}\left[\left(\varphi ^{a\alpha }\varphi ^{a\beta }\right)^2\left(\varphi ^{a\alpha }\varphi ^{a\alpha }\right)^2\right]$$
(1)
where $`1`$ $``$ $`a`$ $``$ $`N`$, $`1`$ $``$ $`\alpha `$ $``$ $`m`$ and $`\overline{g}`$ and $`g_T`$ are independent coupling constants. Alternatively, it can be reformulated in terms of the auxiliary fields $`\sigma `$ and $`T^{\alpha \beta }`$ by, ,

$$L=\frac{1}{2}^\mu \varphi ^{a\alpha }_\mu \varphi ^{a\alpha }+\frac{1}{2}\sigma \varphi ^{a\alpha }\varphi ^{a\alpha }+\frac{1}{2}T^{\alpha \beta }\varphi ^{a\alpha }\varphi ^{a\beta }\frac{3\sigma ^2}{2g}\frac{3T^{\alpha \beta }T^{\alpha \beta }}{2g_T}$$
(2)
where in our notation $`\overline{g}`$ $`=`$ $`g`$ $`+`$ $`(m1)g_T/m`$ and $`T^{\alpha \beta }`$ is an $`O(m)`$ symmetric and traceless auxiliary field. Following the treatment of  there are in fact four fixed points in (2) which are the usual Gaussian and Heisenberg ones together with the stable chiral and unstable antichiral fixed points. For the purposes of this article we are interested in the stable chiral fixed point. In the large $`N`$ critical point approach of  the propagators of the fields of (2) take the following asymptotic scaling forms in the critical region in momentum space

$$\varphi (k)\frac{A}{(k^2)^{\mu \alpha }},\sigma (k)\frac{B}{(k^2)^{\mu \beta }},T^{\alpha \beta ,\sigma \rho }(k)=X^{\alpha \beta ,\sigma \rho }T(k)$$
(3)
where

$$T(k)\frac{C}{(k^2)^{\mu \gamma }}$$
(4)
and

$$X^{\alpha \beta ,\sigma \rho }=\frac{1}{2}\left(\delta ^{\alpha \sigma }\delta ^{\beta \rho }+\delta ^{\alpha \rho }\delta ^{\beta \sigma }\frac{2}{m}\delta ^{\alpha \beta }\delta ^{\sigma \rho }\right)$$
(5)
with the propagator of the field denoted by the same letter. The quantities $`A`$, $`B`$ and $`C`$ are the momentum independent amplitudes and the exponents are defined by

$$\alpha =\mu 1+\frac{\eta }{2},\beta =2\eta \chi ,\gamma =2\eta \chi _T$$
(6)
where $`d`$ $`=`$ $`2\mu `$ and $`\chi `$ and $`\chi _T`$ are the respective anomalous dimensions of the $`\sigma `$ and $`T^{\alpha \beta }`$ vertices of (2). With these forms, (3) and (4), various critical exponents have already been computed in  in the $`1/N`$ expansion. Within this large $`N`$ formalism the four fixed points can be identified with whether the auxiliary fields propagate or not. For example, at the Gaussian fixed point the field theory is free and therefore neither $`\sigma `$ nor $`T^{\alpha \beta }`$ propagate. Likewise the Heisenberg model is clearly recovered from (2) if $`T^{\alpha \beta }`$ is set to zero. For the other non-trivial fixed points, the unstable antichiral one is characterized by $`\sigma `$ $`=`$ $`0`$ and $`T^{\alpha \beta }`$ $``$ $`0`$ and the stable chiral point is given by $`\sigma `$ $``$ $`0`$ and $`T^{\alpha \beta }`$ $``$ $`0`$. Therefore, we will compute the anomalous dimension of the chiral operator for the full Lagrangian, (2). This bilinear operator, $`C^{ab,\alpha \beta }`$, is defined by, ,

$$C^{ab,\alpha \beta }=\varphi ^{a\alpha }\varphi ^{b\beta }\varphi ^{a\beta }\varphi ^{b\alpha }$$
(7)
and is traceless. Denoting the anomalous dimension exponent of this operator by $`\eta _c`$ then the chiral exponent, $`\varphi _c`$, of  is defined by

$$\varphi _c=(2\eta _c)\nu $$
(8)
where

$$\eta _c=\eta +\eta _𝒪.$$
(9)
The exponent $`\eta _𝒪`$ corresponds to the anomalous dimension of the bare operator $`C^{ab,\alpha \beta }`$ itself.

In the large $`N`$ determination of the critical exponent of $`C^{ab,\alpha \beta }`$ the operator is inserted in a two point function where the lines of the graphs are replaced by the massless propagators (3). As we are interested in $`\eta _𝒪`$ at $`O(1/N^2)`$ there is only one topology at $`O(1/N)`$ and six basic topologies at $`O(1/N^2)`$. By topology we mean the independent type of diagrams built from a $`\sigma \varphi ^2`$ style vertex. Of these $`O(1/N^2)`$ topologies four are $`2`$ loop and two are $`3`$ loop. This mixture in number of loops is a consequence of the way one counts powers of $`1/N`$. In the critical point method each combination of amplitudes $`z`$ $`=`$ $`A^2B`$ and $`y`$ $`=`$ $`A^2C`$ counts $`O(1/N)`$ whilst a closed $`\varphi ^{a\alpha }`$ field loop is $`O(N)`$. Ordinarily this would mean that there are topologies up to five loops. However, due to the traceless nature of $`C^{ab,\alpha \beta }`$, diagrams where the operator is inserted in a closed $`\varphi ^{a\alpha }`$ loop yield a zero group theory factor. In addition to substituting the critical propagators into the contributing Feynman diagrams a large $`N`$ regularization needs to be introduced. This is provided by the formal shift $`\chi `$ $``$ $`\chi `$ $`+`$ $`\mathrm{\Delta }`$, $`\chi _T`$ $``$ $`\chi _T`$ $`+`$ $`\mathrm{\Delta }`$. Therefore, each diagram when computed will involve double and simple poles in the parameter $`\mathrm{\Delta }`$ which plays a similar role to $`ϵ`$ in dimensional regularization. These are renormalized in a fashion akin to conventional perturbation theory. Consequently the scaling behaviour of the remaining finite Green’s function is exploited to deduce $`\eta _𝒪`$ at $`O(1/N^2)`$, . This method which we have summarized has been used to evaluate a variety of operator dimensions in several models. For instance, we have written a computer algebra programme in the symbolic manipulation language Form, , to compute $`\eta _𝒪`$. As a non-trivial check on its correctness it has been used to re-evaluate the exponent $`\chi `$ at $`O(1/N^2)`$ at the stable chiral fixed point. Moreover, we have used the same programme to compute the crossover exponent at $`O(1/N^2)`$ at the Heisenberg fixed point where the resulting exponent agreed with four loop $`\overline{\text{MS}}`$ perturbation theory, .

The result of our computation is

$`\eta _𝒪`$ $`=`$ $`{\displaystyle \frac{\mu \eta _1}{(m+1)(\mu 2)N}}`$
$`+[({\displaystyle \frac{\mu ^2(2\mu 1)}{2(m+1)(\mu 2)^2}}{\displaystyle \frac{\mu (2\mu ^25\mu +4)}{2(\mu 2)^2}})v^{}\mu +{\displaystyle \frac{4}{(\mu 2)}}+{\displaystyle \frac{3}{(\mu 2)^2}}`$
$`+{\displaystyle \frac{2}{(\mu 2)^3}}{\displaystyle \frac{3}{2(\mu 1)}}+{\displaystyle \frac{\mu ^2(2\mu ^39\mu ^2+9\mu 1)}{2(m+1)(\mu 1)(\mu 2)^3}}]{\displaystyle \frac{\eta _1}{(m+1)N^2}}+O({\displaystyle \frac{1}{N^3}})`$
where we use the same notation and definitions of  with

$$v^{}=\psi (2\mu 2)+\psi (2\mu )\psi (\mu 2)\psi (2)$$
(11)
and $`\psi (x)`$ $`=`$ $`d\mathrm{ln}\mathrm{\Gamma }(x)/dx`$ where $`\mathrm{\Gamma }(x)`$ is the Euler gamma function. The anomalous dimension $`\eta `$ is given by

$$\eta =\underset{i=1}{\overset{\mathrm{}}{}}\frac{\eta _i}{N^i}$$
(12)
where, ,

$$\eta _1=\frac{2(m+1)\mathrm{\Gamma }(2\mu 2)}{\mathrm{\Gamma }(\mu +1)\mathrm{\Gamma }(\mu 1)\mathrm{\Gamma }(\mu 2)\mathrm{\Gamma }(2\mu )}.$$
(13)
Using the scaling relations (8) and (9), we find

$`\varphi _c`$ $`=`$ $`{\displaystyle \frac{1}{(\mu 1)}}+{\displaystyle \frac{\mu (3m+2)\eta _1}{2(m+1)(\mu 1)(\mu 2)N}}`$ (14)
$`+[({\displaystyle \frac{3\mu ^2(5\mu ^214\mu +10)}{2(\mu 1)(\mu 2)^3}}+{\displaystyle \frac{3\mu ^2(m+3)(3\mu 4)}{2(m+1)^2(\mu 2)^3}})R_1`$
$`{\displaystyle \frac{\mu ^2(m^2+3m+4)(2\mu 3)^2}{2(m+1)^2(\mu 1)(\mu 2)^3}}\left[R_2+R_3^2\right]`$
$`+(\mu (6\mu 7)(3\mu 5)+{\displaystyle \frac{2(2\mu ^311\mu ^2+12\mu 2)}{(m+1)}}`$
$`+{\displaystyle \frac{\mu (2\mu ^227\mu +29)}{(m+1)^2}}){\displaystyle \frac{\mu v^{}}{4(\mu 1)^2(\mu 2)^2}}`$
$`{\displaystyle \frac{1}{2}}+{\displaystyle \frac{25}{4(\mu 2)}}{\displaystyle \frac{23}{(\mu 2)^2}}{\displaystyle \frac{7}{(\mu 2)^3}}{\displaystyle \frac{8}{(\mu 1)}}`$
$`+{\displaystyle \frac{7}{4(\mu 1)^2}}+{\displaystyle \frac{3}{2(\mu 1)^3}}{\displaystyle \frac{\mu (27\mu ^383\mu ^2+63\mu 6)}{4(m+1)(\mu 1)^2(\mu 2)^3}}`$
$`{\displaystyle \frac{\mu ^2(6\mu ^3+9\mu ^267\mu +49)}{4(m+1)^2(\mu 1)^2(\mu 2)^3}}]{\displaystyle \frac{\eta _1^2}{N^2}}+O({\displaystyle \frac{1}{N^3}})`$
where

$`R_1`$ $`=`$ $`\psi ^{}(\mu 1)\psi ^{}(1)`$
$`R_2`$ $`=`$ $`\psi ^{}(2\mu 3)\psi ^{}(2\mu )\psi ^{}(\mu 1)+\psi ^{}(1)`$
$`R_3`$ $`=`$ $`\psi (2\mu 3)+\psi (2\mu )\psi (\mu 1)\psi (1).`$ (15)
For completeness we note that at the stable chiral fixed point, ,

$`\eta _2`$ $`=`$ $`[({\displaystyle \frac{1}{2}}\mu {\displaystyle \frac{1}{(\mu 2)}}{\displaystyle \frac{\mu (m+3)(2\mu 1)}{2(m+1)^2(\mu 2)}})v^{}`$ (16)
$`{\displaystyle \frac{(\mu 1)}{2\mu }}\mu {\displaystyle \frac{2}{(\mu 2)}}+{\displaystyle \frac{1}{2(\mu 1)}}{\displaystyle \frac{\mu (m+3)(2\mu ^25\mu +1)}{2(m+1)^2(\mu 2)^2}}].`$
The expression for $`\varphi _c`$ at $`O(1/N^2)`$ represents the main result of this article and to $`O(1/N)`$ it is in exact agreement with equation (11) of . It can be expanded in $`d`$ $`=`$ $`4`$ $``$ $`2ϵ`$ dimensions in preparation for comparison with explicit $`\overline{\text{MS}}`$ perturbative calculations when they become available. Therefore, we have

$`\varphi _c`$ $`=`$ $`1+ϵ+ϵ^2+ϵ^3+ϵ^4(3m+2)\left[ϵϵ^32(1\zeta (3))ϵ^4\right]{\displaystyle \frac{1}{N}}`$ (17)
$`+[(3m^28m44)ϵ+{\displaystyle \frac{1}{2}}(24m^225m207)ϵ^2`$
$`\left(12(2m^2+7m+11)\zeta (3)+{\displaystyle \frac{1}{4}}(20m^2+123m+339)\right)ϵ^3`$
$`+(80(m^2+3m+4)\zeta (5)18(2m^2+7m+11)\zeta (4)+4(7m^2+58m+149)\zeta (3)`$
$`{\displaystyle \frac{1}{8}}(218m^2+411m+273))ϵ^4]{\displaystyle \frac{1}{N^2}}+O({\displaystyle \frac{ϵ^5}{N^3}})`$
for arbitrary $`m`$ where the order symbol indicates the terms omitted in the independent expansions in powers of $`ϵ`$ and $`1/N`$ and $`\zeta (n)`$ is the Riemann zeta function. For $`m`$ $`=`$ $`2`$ this reduces to

$`\varphi _c`$ $`=`$ $`1+ϵ+ϵ^2+ϵ^3+ϵ^48\left[ϵϵ^32(1\zeta (3))ϵ^4\right]{\displaystyle \frac{1}{N}}`$ (18)
$`+[48ϵ{\displaystyle \frac{161}{2}}ϵ^2({\displaystyle \frac{665}{4}}+396\zeta (3))ϵ^3`$
$`+(1120\zeta (5)594\zeta (4)+1172\zeta (3){\displaystyle \frac{1976}{8}}ϵ^4)]{\displaystyle \frac{1}{N^2}}+O({\displaystyle \frac{ϵ^5}{N^3}}).`$
The expansion to $`O(1/N)`$ illustrates a general property of the general expression for $`\varphi _c`$ at this order. Comparing (14) with the critical exponent of the corresponding bilinear operator which relates to crossing behaviour at the Heisenberg fixed point, , it is in exact agreement for $`m`$ $`=`$ $`2`$, with the differences in values only appearing at $`O(1/N^2)`$ in arbitrary $`d`$. Thus the $`ϵ`$ expansions for both exponents at $`O(1/N)`$ will be the same.

Next we can evaluate our expressions for $`m`$ $`=`$ $`2`$ in three dimensions to obtain

$`\varphi _c`$ $`=`$ $`2{\displaystyle \frac{32}{\pi ^2N}}{\displaystyle \frac{56[2\pi ^2+7]}{\pi ^4N^2}}+O\left({\displaystyle \frac{1}{N^3}}\right)`$
$`\eta _c`$ $`=`$ $`{\displaystyle \frac{232}{3\pi ^4N^2}}+O\left({\displaystyle \frac{1}{N^3}}\right)`$
$`\eta _𝒪`$ $`=`$ $`{\displaystyle \frac{4}{\pi ^2N}}{\displaystyle \frac{56}{\pi ^4N^2}}+O\left({\displaystyle \frac{1}{N^3}}\right)`$ (19)
where the intermediate exponents are given by

$`\eta `$ $`=`$ $`{\displaystyle \frac{4}{\pi ^2N}}{\displaystyle \frac{64}{3\pi ^4N^2}}+O\left({\displaystyle \frac{1}{N^3}}\right)`$
$`\nu `$ $`=`$ $`1{\displaystyle \frac{16}{\pi ^2N}}{\displaystyle \frac{8[21\pi ^2+88]}{3\pi ^4N^2}}+O\left({\displaystyle \frac{1}{N^3}}\right).`$ (20)
It is interesting to note that there is no $`O(1/N)`$ correction to $`\eta _c`$. There are two related chiral exponents defined by the scaling relations in arbitrary dimensions

$$\beta _c=2\mu \nu \varphi _c,\gamma _c=2\varphi _c2\mu \nu .$$
(21)
Using (19), (20) and (21) we find

$`\beta _c`$ $`=`$ $`1{\displaystyle \frac{16}{\pi ^2N}}{\displaystyle \frac{8[7\pi ^2+39]}{\pi ^4N^2}}+O\left({\displaystyle \frac{1}{N^3}}\right)`$
$`\gamma _c`$ $`=`$ $`1{\displaystyle \frac{16}{\pi ^2N}}{\displaystyle \frac{8[7\pi ^2+10]}{\pi ^4N^2}}+O\left({\displaystyle \frac{1}{N^3}}\right)`$ (22)
in three dimensions. From the results for each of $`\varphi _c`$, $`\beta _c`$ and $`\gamma _c`$ it is easy to observe that the coefficient of the $`O(1/N^2)`$ term is large and therefore each series is diverging. One way to improve the convergence is to apply the Padé-Borel resummation technique. For a series of the form $`_{n=0}^{\mathrm{}}a_nx^n`$ it can be rewritten in terms of the associated Borel function

$$\underset{n=0}{\overset{\mathrm{}}{}}a_nx^n=\frac{1}{x}_0^{\mathrm{}}𝑑te^{t/x}\underset{n=0}{\overset{\mathrm{}}{}}\frac{a_nt^n}{n!}$$
(23)
and then the integrand can be replaced by a Padé approximant if a limited number of coefficients $`\{a_n\}`$ are known. Hence, for $`\varphi _c`$ we have

$$\varphi _c=2N_0^{\mathrm{}}𝑑t\frac{e^{Nt}}{\left[1a_1t+(a_1^2\frac{1}{2}a_2)t^2\right]}$$
(24)
where

$$a_1=\frac{16}{\pi ^2},a_2=\frac{28[2\pi ^2+7]}{\pi ^4}.$$
(25)
We have evaluated this estimate for $`\varphi _c`$ for various values of $`N`$ and recorded them in Table 1 together with the estimates for $`\beta _c`$ and $`\gamma _c`$ using the same method.

One has to be careful in interpreting the results in Table 1. Unlike the widely studied Heisenberg model, it is believed that for low values of $`N`$ the stable chiral fixed point of (1) is not analytically connected to the one used in the large $`N`$ expansion. Indeed not only does the work of  suggest this but recent exact renormalization group analysis corroborates this position, . In other words there is evidence that the large $`N`$ expansion is not predictive for $`N`$ $`=`$ $`2`$ and $`3`$, . Nevertheless if one sets this point of view aside for $`N`$ $`=`$ $`2`$ and $`3`$ then the results of Table 1 follow from (24). On the other hand the values for say $`N`$ $`=`$ $`6`$ and larger would be expected to be more reliable. Though there are no experimental or theoretical results to compare with in these cases.

To conclude with we have provided the $`O(1/N^2)`$ expressions for the chiral exponents of the $`O(N)`$ $`\times `$ $`O(m)`$ $`\varphi ^4`$ theory in arbitrary dimensions at the stable chiral fixed point. At $`O(1/N)`$ the result for $`m`$ $`=`$ $`2`$ is in exact agreement with the expression of the crossover exponent at the Heisenberg fixed point in the same model. However, in three dimensions for the specific case $`m`$ $`=`$ $`2`$, the $`O(1/N^2)`$ corrections are large and therefore the series is not converging quickly. Indeed comparison with the crossover exponent, , shows that the convergence for the $`m`$ $`=`$ $`2`$ chiral exponent is slower and gets slower as $`m`$ increases.


format: text

text:
# First Principles Force Field for Metallic Tantalum

## I Introduction

Despite decades of experimental and theoretical research on the mechanical properties of materials many questions remain open, particularly in the relation between atomistic processes (involving dislocations, grain boundaries, cracks, their mobility and interactions thereof) and the macroscopic behavior (plastic deformation, failure, etc.). Macroscopic plasticity and failure are well characterized experimentally and described using a variety of mesoscale and macroscale models with parameters obtained empirically. These models and their parameters should ultimately be derivable in terms of the fundamental physics of atomic interactions as described by quantum mechanics (QM). Unfortunately despite the enormous progress in ab initio QM, such calculations are too computationally demanding to study directly most processes relevant to plasticity and failure. To do so would require that millions of atoms be described for nanoseconds or longer, calculations that are impossible to consider today. In order to bridge this gap between atomic interactions and the mechanical properties of macroscopic systems we use first principles QM data to derive a FF with which energies and forces can be calculated in a computationally efficient way given only the atomic positions. This allows us to use classical molecular dynamics (MD) to simulate the various atomistic processes governing the mechanical and thermodynamical properties of materials. Since the FF is derived by fitting a wide range of QM data, we expect to obtain an accurate description of the atomic interactions which, with molecular dynamics, can provide insight and constitutive equations to be used in the mesoscale and macroscale models of plasticity and failure. We illustrate this procedure here for Ta. Ta is chosen because it exhibits only a single crystalline phase over the interesting range of temperatures and pressures. This makes the validation of theoretical predictions against experimental data under a wide range of conditions more clear.

Several authors have used ab initio methods to study various mechanical properties of Ta. Söderlind and Moriarty used the full potential linear muffin-tin orbital method within the GGA approximation and with spin orbit interactions to calculate different zero temperature properties of of Ta, including the EOS of different crystalline phases, elastic constants, and shear strength.  Vacancy formation migration energies have also been calculated from first principles.  Recently, Ismail-Beigi and Arias  calculated ab initio dislocation core energy and structures for Ta and Mo.

On the other hand, several FF have been developed for bcc metals. The first principles based multi-ion analytic model generalized pseudopotential (MGPT) force field has been used to calculate several mechanical and thermodynamical properties of various metals,  including Ta. One of the most popular many-body force fields for metals is the Embedded Atom Model (EAM), proposed by Daw and Baskes in 1984.  Most EAM FF have been based on experimental data regarding structures at and near equilibrium \[for bcc metals, see for example Ref. \]. One of the main advantages of EAM force fields is that they are computationally very efficient which allows MD simulations of large systems for long times. The modified embedded atom model (MEAM), which includes angular dependence of the electronic density, was developed by Baskes and coworkers  and applied to a variety of materials including bcc transition metals .

In this paper we propose a strategy to derive accurate many body FF based on QM calculations; these FF can be used to simulate the dynamics of systems with millions of atoms. We use accurate QM (LAPW GGA method) to calculate various mechanical properties of Ta which require a small number of atoms: namely the zero temperature EOS for bcc, fcc and hcp crystalline phases, and elastic constants. We also use a mixed-basis pseudopotential method to calculate vacancy formation energies. We find that we can describe all this first principles data using a classical many body EAM FF (named qEAM FF) with good overall accuracy. We then illustrate the use of the qEAM FF with MD simulations to study various properties as a function of pressure and temperature; such as the melting curve of Ta up to pressures of $`300`$ GPa, and thermal expansivity.

This paper is organized as follows. In Section II we present first principles results for the equation of state for bcc, fcc and hcp phases, elastic constants and volume relaxed vacancy formation energy and enthalpy for a wide pressure range. In Section III we develop the qEAM FF based on ab initio results. In Section IV we calculate the thermal expansivity of Ta using the qEAM FF and ab initio calculations. Section V presents the calculation of the melting curve of Ta for pressures up to $`300`$ GPa using molecular dynamics. Finally, in section VI, conclusions are drawn.

## II Quantum mechanics results

We computed the static equation of state of Ta for different crystalline phases using the linearized augmented plane wave (LAPW) method.  LAPW is an all-electron method, with no essential shape approximations for the charge density or potential, and is easily converged. The $`5p`$, $`4f`$, $`5d`$ and $`6s`$ states were treated as band states, and the deeper states were treated as soft core electrons. Here we used the Perdew, Burke, and Ernzerhof (PBE) implementation of the generalized gradient approximation  for the exchange-correlation potential. A 16x16x16 special k-point mesh was used, giving 140 k-points within the irreducible Brillouin zone of the bcc lattice. Tests demonstrated convergence with this mesh. The convergence parameter RKmax was 9 giving about 1800 planewaves and 200 basis functions per atom at zero pressure.

Total energies were computed for bcc, fcc and hcp phases at 20 volumes. For the fcc and hcp phases, 12x12x12 and 16x16x12 k-point meshes were used for Brillouin zone integrations giving 182 and 180 k-points within the irreducible zone respectively. For the hcp phase, the ideal c/a ratio was used, and at two different volumes also c/a was optimized. We found that the change in the energy due to the this optimization is less than 40 meV/atom around zero pressure and decreases by pressure. We have tabulated the ab initio energy-volume data for the different phases and have made them available as supplementary material, see Ref. . We have fitted our energy-volume data to Rose’s universal equation of state;  the obtained zero pressure volume ($`V_0`$), zero temperature bulk modulus ($`B_T`$), and its derivative with respect to pressure ($`B_T^{^{}}`$) are shown in Table I. We also show in Table I the results obtained by Söderlind and Moriarty using full potential linear muffin-tin orbital method within the GGA approximation and with spin orbit interactions (denoted as FP LMTO GGA SC), and room temperature experimental values by Cynn and Yoo.  Our LAPW calculations of the bcc equation of state agree well with the experimental values and previous theoretical calculations.

Static elastic constants \[cs=(c11-c12)/2 and c44\] were obtained from strain energies by straining the bcc cell with volume conserving tetragonal and orthorhombic. We calculate cs using tetragonal strain of the cubic bcc lattice:

$`𝐚`$ $`=`$ $`a(1+ϵ,0,0),`$ (1)
$`𝐛`$ $`=`$ $`a(0,1+ϵ,0),`$ (2)
$`𝐜`$ $`=`$ $`a(0,0,1/(1+ϵ)^2),`$ (3)
where $`a`$ is the cubic lattice constant of the system, $`𝐚`$, $`𝐛`$, and $`𝐜`$ are the lattice vectors and $`ϵ`$ is the strain. $`c_s`$ is related to the quadratic term of the strain energy:

$$E(ϵ)=E(ϵ=0)+6V(ϵ=0)c_sϵ^2+O(ϵ^3),$$
(4)
where $`E(ϵ=0)`$ is the energy of the unstrained system and $`V(ϵ=0)`$ is its volume. Similarly c44 is obtained from the orthorhombic strain:

$`𝐚`$ $`=`$ $`a(1,ϵ,0),`$ (5)
$`𝐛`$ $`=`$ $`a(ϵ,1,0),`$ (6)
$`𝐜`$ $`=`$ $`a(0,0,1/(1ϵ^2)).`$ (7)
the shear constant $`c_{44}`$ is obtained from:

$$E(ϵ)=E(ϵ=0)+2V(ϵ=0)c_{44}ϵ^2+O(ϵ^3).$$
(8)
The convergence of strain energies with respect to the Brillouin zone integration was carefully checked; we used 16x16x16 k-points meshes in the full Brillouin zone giving 344 and 612 k-points within the irreducible Brillouin zone of tetragonal and orthorhombic lattice respectively. $`c_s`$ and $`B_T=(c_{11}+2c_{12})/3`$ were used to calculate $`c_{11}`$ and $`c_{12}`$; the resulting zero pressure and zero temperature elastic constants are shown in Table I. The elastic constants $`c_s`$ and $`c_{44}`$ as functions of pressure are available as supplementary material, see also Ref. . The zero pressure values and initial slopes are in good agreement with the experimental data of Katahara et al. We find that c44 shows a change in behavior at $``$150 GPa  which is probably due to the electronic transition also evident in the equation of state. The band structure and density of state show a major reconfiguration with pressure. Our results indicate that the elastic constants can be much more sensitive to changes in the occupied states below the Fermi surface than the equation of state, where changes were much more subtle.

The vacancy formation energy was determined  using supercells of 16 and 54 atoms and a mixed-basis pseudopotential (MB-PS) code. Using the GGA approximation,  the zero pressure volume relaxed vacancy formation energy was determined to be 3.25 eV for a 54 atom supercell, and 3.26 eV for a 16 atom supercell, indicating convergence. The effects of structural relaxation on vacancy formation energy are discussed in Ref. . A comparison with recent ab initio and experimental results will be presented in the next section. The ab initio data used to calculate vacancy formation energies is also available as supplementary material.

## III qEAM force field

In order to predict mechanical properties of materials and processes like phase diagrams, dislocations structures and mobility, mechanical failure, etc. it is important to have accurate classical force fields to describe the atomic interactions. With MD simulations it is then possible to study large systems (millions of atoms) for relatively long times (ns).

One of the most popular many-body force fields for metals is the EAM, proposed by Daw and Baskes in 1984.  This approach is computationally efficient and has been used successfully for numerous applications, like calculation of thermodynamic functions, liquid metals, defects, grain boundary structure, fracture, etc, see for example Ref. .

The EAM implementation that we have used is based on the one proposed by Chantasiriwan and Milstein.  This from of EAM was chosen because it can describe third order elastic constants correctly  thus being useful in a wide range of strains. The energy of a given atomic configuration with atom positions $`\{r_i\}`$ is given by:

$$U\{r_i\}=\underset{i}{\overset{N}{}}F(\rho _i)+\underset{i<j}{}\varphi (r_{ij}),$$
(9)
with

$$\rho _i=\underset{ji}{}f(r_{ij})$$
(10)
where $`F(\rho )`$ is the embedding energy, $`\rho _i`$ is the total “electronic density” on the atomic site $`i`$, $`f(r)`$ is the electron density function, $`\varphi (r)`$ is a two-body term and $`r_{ij}=|r_ir_j|`$.

Following Ref.  we took the two-body term as:

$$\varphi (r)=\{\begin{array}{cc}(rr_m)^4_{i=0}^7b_ir^i\hfill & \text{if }r<r_m\hfill \\ 0\hfill & \text{otherwise,}\hfill \end{array}$$
(11)
the factor $`(rr_m)^4`$ ensures that the two body term and its first three derivatives with respect to $`r`$ vanish at the the cut-off radius ($`r_m`$). The optimized form of the two body term shows short range repulsion and longer range attraction.

The “electron density” is:

$$f(r)=\frac{1+a_1\mathrm{cos}(\frac{\alpha r}{V^{1/3}})+a_2\mathrm{sin}(\frac{\alpha r}{V^{1/3}})}{r^\beta },$$
(12)
where $`V`$ is the volume per atom of the system. The parameters of the “electron density” are volume dependent, but structure independent. The importance of the oscillatory behavior of the “electron density function” in embedded atom model-like force fields is related to their ability to correctly account anharmonicities.

Finally the embedding energy as a function of the electronic density is obtained from the reference bcc structure:

$$F(\rho )=U_{Rose}(V)\underset{i<j}{}\varphi (r_{ij}),$$
(13)
where the sum is made for the perfect bcc structure and $`U_{Rose}(V)`$ is Rose’s universal equation of state:

$$U_{Rose}(V)=E_{coh}(1+a^{}+f_3a^3+f_4a^4)e^a^{}$$
(14)
where $`a^{}=(aa_0)/a_0\lambda `$ and $`\lambda =(E_{coh}/9V_0B_T)^{1/2}`$.

We we optimized the parameters in the qEAM FF using the following data:

(a) Zero temperature energy-volume and pressure-volume curves for different crystal structures (bcc, fcc and A15) in a wide pressure range, from $`10`$ GPa to $`500`$ GPa. For the bcc and fcc structures the data from Section II was used, while for A15 we used the results obtained by Söderlind and Moriarty.

(b) Zero temperature, zero pressure elastic coefficients, shown in section II.

(c) Vacancy formation energy at zero pressure.

(d) Energetic of homogeneously sheared bcc crystal, from Ref. .

(c) Unrelaxed (100) surface energy from Ref. .

We fit the parameters entering the qEAM FF energy expression to the training set using an optimization algorithm based on simulated annealing. We define an error function of the form:

$$C=\underset{i}{}C_i\frac{(Q_i^{qEAM}Q_i^{target})^2}{(Q_i^{target})^2},$$
(15)
where the sum runs over the different quantities to be fitted, $`Q_i^{qEAM}`$ is the value given by the qEAM FF for quantity $`i`$, $`Q_i^{target}`$ is the target quantity and $`C_i`$ are weight factors. The EOS information of the different phases account for $`78\%`$ of the total cost; elastic constants and vacancy formation energy account for $`9\%`$ each one; surface energy represents $`3\%`$ of the total cost and energy of sheared bcc crystal represent $`1\%`$.

In Tables II and III we show the optimized qEAM parameters.

A detailed comparison between the qEAM FF and the data it was fitted to is shown in the following subsections, together with other ab initio and experimental results. We also show the calculation of related quantities obtained using the qEAM FF via MD simulations.

### A Equation of state and elastic constants

The most important quantities used to develop the qEAM FF are zero temperature EOS for different crystal phases of Ta in a wide pressure range. We used energy-volume and pressure-volume for bcc Ta from Section II; fcc-bcc energy difference and fcc pressure for different volumes from Section II; and first principles A15-bcc energy difference obtained by Söderlind and Moriarty  using full potential linear muffin-tin orbital method within the GGA approximation with spin orbit interactions. Ta is a bcc metal and no pressure-induced phase transition to other solid structure has been found experimentally or theoretically. Nevertheless, using QM it is possible to calculate the equation of state for different crystalline structures although they may not be thermodynamically stable. Including the EOS of different phases in the FF development leads to an accurate description of the atomic interactions even when the environment of an atom is not that of the stable phase; this could play a key role to correctly describe defects and non-equilibrium processes.

In Figure 1 we show energy \[Figure 1(a)\] and pressure \[Figure 1(b)\] as a function of volume for bcc Ta at $`T=0`$ K. The circles denote LAPW GGA results and the lines the qEAM FF. Figures 2 and 3 show the same results for fcc and A15 Ta respectively. In Figure 3(a) open circles show the A15 energy calculated by Söderlind and Moriarty,  and the filled circles denote the sum of the A15-bcc energy difference from Ref.  and our bcc energy from section II, these are the quantities the FF was actually fitted to. The insets on Figures 2(a) and 3(a) show the fcc-bcc and A15-bcc energy difference as a function of volume obtained with the qEAM FF (lines) and from QM (circles). In Figure 4(a) and 4(b) we show energy-volume and pressure-volume curves for hcp phase respectively; circles denote denote the LAPW GGA results of section II and the lines show the qEAM FF results. Note that although the the hcp data was not included in the training set for the qEAM FF the agreement is very good. Figures 1 to 4 show that the qEAM FF reproduces the zero temperature EOS for the four different phases very well.

We also included in the FF training set the ab initio elastic constants from Section II at zero pressure. Table I shows bcc Ta EOS parameters \[zero pressure volume ($`V_0`$), bulk modulus ($`B_T`$), its derivative with respect to pressure ($`B_T^{}`$) and the elastic constants\] obtained using the qEAM FF together with the QM values from Section II and the ones reported in Ref. . $`V_0`$, $`B_T`$, and $`B_T^{}`$ were obtained fitting Rose’s universal equation of state to the energy-volume data shown in Figure 1. The elastic constants $`c_s`$ and $`c_{44}`$ were calculated with the qEAM FF using the tetragonal and orthorhombic strains shown above \[equations (1) to (4)\].

In Figure 5 we show the elastic constants \[bulk modulus $`B_T`$, $`c_s`$ and $`c_{44}`$\] as a function of pressure obtained with the qEAM FF (filled circles and full lines) and the LAPW results from Section II. While the agreement in $`B_T`$ is excellent and that for $`c_s`$ is good, qEAM FF greatly underestimates $`c_{44}`$ for high pressures; this problem is amplified by the possible electronic phase transition that leads to a change of behavior of $`c_{44}`$ at $``$ 150 PGa (see Section II and Ref. ).

In order to estimate the accuracy of the the numerical calculation of elastic constants we have calculated $`c_{11}`$ independently using a uniaxial strain:

$`𝐚`$ $`=`$ $`a(1+ϵ,0,0),`$ (16)
$`𝐛`$ $`=`$ $`a(0,1,0),`$ (17)
$`𝐜`$ $`=`$ $`a(0,0,1),`$ (18)
the elastic constant $`c_{11}`$ is obtained from:

$$E(ϵ)=E(ϵ=0)+PV(ϵ=0)ϵ+\frac{1}{2}V(ϵ=0)c_{11}ϵ^2+O(ϵ^3),$$
(19)
where $`E(ϵ=0)`$ is the zero strain energy, $`V(ϵ=0)`$ is the volume at zero strain, and P is the zero strain hydrostatic pressure. For zero pressure we calculate $`c_{11}=273.67GPa`$ only 0.4 $`\%`$ higher than the one computed from $`c_s`$ \[equations (1) and (2)\] and $`B_T`$ (calculated form the energy-volume data shown in Figure 1). We find a similar agreement under compression: for pressure P=109.54 GPa (corresponding to a volume per atom $`V(ϵ=0)=13.36\AA ^3`$) we obtain $`c_{11}=819.90`$ GPa only 0.8 $`\%`$ larger than the calculated from $`c_s`$ and $`B_T`$.

We have calculated the T$`=300`$ K EOS using isothermal-isobaric (NPT) MD with a Hoover  thermostat and Rahman-Parrinello barostat.  In Table I we show the T=$`300`$ K zero pressure volume, bulk modulus and its first derivative with respect to volume; we also show recent compressibility data  obtained in a diamond-anvil cell at room temperature and ultrasonic measurements of adiabatic elastic constants.

Zero temperature EOS data obtained with the qEAM FF is available in the supplementary material.

### B Vacancy formation energy

We used experimental values for vacancy formation energy and cohesive energy in the training set used to fit the qEAM FF. The experimental value for the relaxed vacancy formation energy is $`E_{vac}=2.8`$ eV.  From this value we estimated the value of the unrelaxed vacancy formation energy to be $`3`$ eV. This unrelaxed value was used to fit the qEAM FF. The value for the cohesive energy used is $`E_{coh}=8.10`$ eV.

The volume-relaxed vacancy formation energy is defined in computer simulations as:

$$e_{vac}(P)=e_{vac}(N1,P)\frac{N1}{N}e_{xtal}(N,P),$$
(20)
where $`e_{xtal}(N,P)`$ is the energy of the perfect crystal with N atoms at pressure $`P`$ and $`e_{vac}(N1,P)`$ is the energy corresponding to the N-1 atoms system with a vacancy at pressure $`P`$ where the atomic positions are not allowed to relax. The relaxed vacancy formation energy is defined in the same way but with $`e_{vac}(N1,P)`$ being the atom-relaxed energy of the system with a vacancy.

The volume-relaxed vacancy formation energy obtained using the qEAM FF is 3.10 eV, in very good agreement with the target value (3 eV) and only slightly lower than the ab initio value $`3.25`$eV. The relaxed vacancy formation energy obtained with the qEAM FF is $`2.95`$ eV. Previous work by Korhonen, Puska and Nieminen, using a DFT full potential linear muffin-tin orbital method with LDA gives $`3.49`$ eV for the unrelaxed vacancy formation energy;  Satta, Willaime and Gironcoli  (using DFT plane Waves LDA) calculated $`3.51`$ eV and $`2.99`$ eV for unrelaxed and relaxed vacancy formation energies; recently Söderlind, Yang, Moriarty and Wills (also using the full potential linear muffin-tin orbital method with GGA) obtained $`3.74`$ eV for the volume relaxed vacancy formation energy and $`3.2`$ for the relaxed vacancy formation energy. Table IV summarizes our vacancy formation energies together with previous theoretical results and experimental data.

In Figure 6(a) we show the volume-relaxed vacancy formation energy ($`e_{vac}`$) as a function of pressure (P); the thick solid line shows qEAM results and circles denote QM results of Section II. We used a simulation cell containing 54 atoms for the perfect crystal case with periodic boundary conditions. Vacancy formation enthalpy is defined in the same way as $`e_{vac}(P)`$ \[eq. (20)\] replacing energy with enthalpy. In Figure 6(b) we plot the vacancy formation enthalpy $`h_{vac}`$ with respect to pressure. The difference between the vacancy formation enthalpy as obtained by qEAM FF and QM, is due to the difference in the vacancy formation volume obtained from the two methods. The vacancy formation volume \[$`v_{vac}(P)`$\] is, again, defined in the same way as $`e_{vac}(P)`$ \[eq. (20)\] replacing energy with volume. In Figure 7 we plot the vacancy formation volume ($`\mathrm{\Omega }_{vac}^{}{}_{}{}^{f}`$), with respect to pressure. While the pressure dependence of $`v_{vac}(P)`$ calculated using the qEAM FF is very similar to the MB-PS calculation, our Force Field overestimates the vacancy formation volume resulting in higher vacancy formation enthalpy for compressed states.

In order to calculate the vacancy formation enthalpy at finite temperatures as a function of pressure we performed NPT MD simulations using a cell containing $`N=1458`$ atoms with periodic boundary conditions. We performed simulations at 9 different volumes (19, 18.36, 18, 17, 16, 15, 14, 13, 11 $`\AA ^3`$); for each volume we started with the relaxed structure a $`T=0`$ K and heat the system in $`100`$ K steps; for each temperature we performed $`25`$ ps MD simulations and used the last $`20`$ ps for measurements. In Figure 8 we show vacancy formation enthalpy as a function of pressure for $`T=0K`$ (full atomic relaxation), $`T=300`$ K, the volume-relaxed enthalpy is also shown for comparison; the zero pressure values are also shown in Table IV. The fundamental data used to calculate vacancy energy, enthalpy and volume can be found in the supplementary material.

A very important quantity, which determines vacancy mobility, is the vacancy migration energy barrier ($`E_{vac}^{mig}`$). We calculate $`E_{vac}^{mig}`$ using the qEAM FF by marching a nearest neighbor atom towards the position of the vacancy in short steps. At each step the position of the marching atom is fixed, as well as that of a distant, reference, atom, and the positions of all the other atoms are relaxed at constant pressure. In this way we obtain the optimum migration path and energy as a function of displacement. In Figure 9 we show the energy as a function the position of the marching atom for zero pressure and zero temperature; we obtain a vacancy migration energy $`E_{vac}^{mig}1.093`$ eV. The activation energy for self diffusion is defined as $`Q=E_{vac}^f+E_{vac}^{mig}`$, using the qEAM FF we obtain $`Q=4.028`$ eV in good agreement with the experimental value $`3.8\pm 0.3`$ eV,  and ab initio calculations  which give $`3.82`$ eV, see Table IV. At T=300 K we obtain a vacancy migration energy of 1.1 $`\pm `$ 0.5 eV, see Table IV, very similar to the zero temperature value.

### C Energetics of homogeneously sheared bcc crystal

Zero temperature, first principles energetics of a homogeneously sheared bcc crystal in the observed twinning mode was also included in the qEAM training set.

The ideal shear strength is defined to be the stress separating elastic and plastic deformation when a homogeneous shear is applied to a perfect crystal. It gives an upper bound for the shear strength of the material. The shear transformation is in the direction of the observed twining mode and deforms the crystal into itself.  For bcc crystal we use the following transformation of the cell vectors:

$`𝐚={\displaystyle \frac{1}{2}}[\overline{1}11]+{\displaystyle \frac{s}{\sqrt{18}}}[\overline{1}\overline{1}1],`$ (21)
$`𝐛={\displaystyle \frac{1}{2}}[1\overline{1}1]+{\displaystyle \frac{s}{\sqrt{18}}}[\overline{1}\overline{1}1],`$ (22)
$`𝐜={\displaystyle \frac{1}{2}}[11\overline{1}],`$ (23)
when the shear variable $`s`$ is equal to the twinning shear $`s=s_{tw}=2^{1/2}`$ the lattice vectors \[$`a=1/3[\overline{2}12]`$, $`b=1/3[1\overline{2}2]`$ and $`c=\frac{1}{2}[11\overline{1}]`$\] form a bcc structure, twin of the initial one.

In this way one can calculate the energy along the shear path,

$$W(s)=e(V,s)e(V,s=0),$$
(24)
where $`e(V,s)`$ is the energy per atom of the deformed system and $`e(V,s=0)`$ is the perfect crystal energy. The energy barrier associated with this transformation is $`W_{max}=W(s=0.5)`$. The corresponding stress is defined as:

$$\tau (s)=\frac{1}{V}\frac{dW(s)}{ds}.$$
(25)
The ideal shear strength ($`\tau _{max}`$) is defined as the maximum stress along the path. In Figure 10 we show energy and stress as a function of shear using the qEAM FF for zero pressure volume, see also Ref. .

Söderlind and Moriarty  calculated $`W(s)`$ and $`\tau (s)`$ for Ta at different volumes, from first principles. In developing the qEAM FF we used $`W_{max}`$ for V=$`17.6186\AA ^3`$ and V=$`10.909\AA ^3`$ as part of the training set. In Table V we show a comparison between the first principles results  and the ones obtained using the qEAM FF. We can see that the qEAM results are is in very good agreement with the ab initio calculations.

### D Surface energy

The unrelaxed (100) surface energy using the qEAM FF is $`1.971`$ J/m2, lower than the first principles values of $`3.097`$ J/m2 of Ref.  and 3.23 J/m of Ref. . Low surface energy is a common problem in EAM-like force fields, see Ref. . The zero temperature experimental estimate of the surface energy (averaged over different surfaces) is 2.902 J/m2.

## IV Thermal expansion

Thermal expansivity is an important materials property that can be calculated directly from MD simulations. In order to calculate the lattice constant as a function of temperature for zero pressure we performed NPT MD simulations with a computational cell containing 1024 atoms increasing the temperature by 100 K every 25 ps at zero applied pressure. For each temperature the first 5 ps were taken as thermalization and the remaining 20 ps were used for measurements.  ab initio MD simulations are very time consuming and only feasible for small systems and short times; thus in order to calculate the thermal expansion from first principles we take the following approach.  The Helmholtz free energy can be written as:

$$F(V,T)=E_0(V)+F_{el}(V,T)+F_{vib},$$
(26)
where $`E_0(V)`$ is the zero temperature energy (Section II), $`F_{el}(V,T)`$ is the electronic contribution and $`F_{vib}`$ is the vibrational part of the free energy. The electronic contribution is obtained using quantum statistical mechanics. The vibrational part is obtained using the particle in a cell method,  further detail of our calculations can be found in Ref. . This calculations were performed using the mixed basis pseudopotential method and a cell containing 54 atoms.

In Figure 11(a) we show the thermal expansivity ($`\alpha (T)=1/VV/T`$) as a function of temperature obtained form our MD simulations (circles), first principles results (dashed-dotted line), as well as the experimental values (full and dashed lines).  Figure 11(b) shows the linear thermal expansion $`(aa_0)/a_0`$ (where $`a_0`$ is the T$`=300`$ K lattice constant) as a function of temperature obtained using the qEAM FF (circles) as well as the experimental values.  The high temperature experimental results (shown as dashed lines) represent provisional experimental data. Both MD and the particle in cell methods are based on classical mechanics so none of the them are expected to capture the low temperature behavior where the differences between quantum and classical statistical mechanics are important. The force field results agree with experimental data well; for example the qEAM FF underestimates the change in lattice parameter from $`T=300`$ K to $`T=2000`$ K by less than $`0.2\%`$. This is an important result since the thermal expansion in related to anharmonicities of the internal energy. Our ab initio data also agree well with the experimental results; we slightly overestimate thermal expansivity for temperatures lower than T=2000 K.

## V Melting curve of Ta

We studied melting of Ta using one phase and two-phase MD simulations with the qEAM force field.

In Figure 12 we show the results of zero pressure heating a solid sample until it melts and then cooling the melt. The system contains 1024 atoms with periodic boundary conditions with a heating and cooling rate of 100 K per 25 ps. In Figure 12 we plot enthalpy \[Figure 12(a)\] and volume \[Figure 12(b)\] as a function of temperature. The solid superheats on heating and the liquid undercools as expected for a small, periodic system with no defects to act as nucleation centers and high heating and cooling rates.

In order to overcome some of these problems we calculate the melting temperature using the “two phase technique”. We place pre-equilibrated liquid and pre-equilibrated solid samples in a single computational cell. Once this initial configuration is built we perform TPN MD (using a Hoover thermostat and Rahman-Parinello barostat) simulation and observe which phase grows. If the simulation temperature T is lower than the melting temperature at the given pressure \[$`T_{melt}(P)`$\] the liquid will start crystallizing, on the other hand the solid will melt if $`T>T_{melt}(P)`$.

Given that the liquid-solid phase transition is first order, i.e. the enthalpy of the solid and liquid in equilibrium differ by a finite amount, it is very easy to know whether the system is crystallizing or melting by analyzing the time evolution of the total potential energy during the MD run. In Figure 13 we show potential energy as a function of time for zero pressure two-phase simulations at different temperatures. We can see that at T=3100 K (dotted line in Figure 13) the potential energy decreases with time; this means that the solid phase is growing and $`T_{melt}>3100`$ K. On the other hand, for $`T=3200`$ K the energy grows, the system is, then, melting and $`T_{melt}<3200`$ K. At T = 3150 K the energy is rather constant and both phases are close to equilibrium. This value for the zero pressure melting is in very good agreement with experimental results which range from $`3213`$ K to $`3273`$ K.  This is a very important validation of the FF, taking into account that only zero temperature data was used in its development.

The slope of the melting curve is given by the Clausius Clapeyron equation:

$$\left[\frac{dP}{dT}\right]_{melt}=\frac{1}{T}\frac{\mathrm{\Delta }H}{\mathrm{\Delta }V},$$
(27)
where $`\mathrm{\Delta }H`$ and $`\mathrm{\Delta }V`$ are the enthalpy and volume difference between the liquid and solid in equilibrium respectively. From our MD runs, see Figure 12, the slope of the melting curve at zero pressure is $`dT/dP=92.8`$ K/GPa, larger than the the experimental value $`dT/dP=60\pm 10`$ K/GPa.

Using the two phase simulation procedure we calculated the melting temperature for various pressures up to $`300`$ GPa. In Figure 14 we show the melting curve for Ta, see also Ref. . To the best of our knowledge this is the first calculation of melting temperature in Ta for such a wide pressure range. Zero pressure experimental values  are also shown in Figure 14 as empty diamonds. High pressure melting of Ta has been studied experimentally via shock compression;  melting is identified as a change in the velocity of the rarefication wave (from the longitudinal to the bulk sound velocity). The transition was found to occur in the pressure range from $``$250 GPa to 295 GPa. The calculation of the temperature in shock experiments (i.e. along the Hugoniot) is difficult; the electronic contribution to the specific heat has a very strong effect on the melting temperature.  Simple models for the electronic behavior lead to very different temperatures: using the free electron gas model the melting temperature is $``$10000 K while considering band electrons give $`T_{melt}`$ 7500 K, see Ref.. These experimental values are shown in Figure 14 are empty circles. Using the accurate ab initio thermal equation of state obtained using the methods described above and the Rankine-Hugoniot equation Cohen and Gülseren calculated pressure-volume and temperature-pressure curves for Ta under shock conditions. This calculation leads to a melting temperature (using the experimental melting pressure Pmelt=375 GPa) of 8150 K (see square in Figure 14). Our MD results are in good agreement with the high pressure calculation of the melting temperature based on shock experiments and ab initio calculations.

## VI Conclusions

We have presented here a general strategy to derive accurate classical force fields based on ab initio QM mechanical calculations for metallic systems. Force fields allow calculations on systems containing millions of atoms, providing a means to study from an atomistic point of view a variety of processes relevant to the mechanical and thermodynamical properties of metals, such as phase transitions, dislocations dynamics and interactions, failure, crack propagation, etc.

We showed that the qEAM FF describes with good accuracy the EOS of Ta for different crystal phases \[bcc (coordination number 8), fcc (coordination number 12), and A15 (mixed coordination numbers)\], and also elastic constants, vacancy formation energy, and energetics of the deformed bcc lattice in the twinning direction. A critical point in our approach is that in developing the force field we use not only the thermodynamically stable phase but also high energy phases and large strains; a force field that can correctly describe such structures should be appropriate for simulations of defects and non-equilibrium processes.

The large amount of QM data used to derive our force field gives an important measure of its quality.

We used the qEAM FF with MD to calculate the melting curve for Ta in a wide pressure range. The zero pressure melting temperature obtained from our simulations $`T_{melt}=3150\pm 50`$ K is in very good agreement with the experimental result of $`3290\pm 100`$ K. This is an important validation of our approach given the fact that the qEAM FF was derived using zero temperature data.

First-priciples-based force fields represent an important step in ab initio multiscale modeling of materials. We have used the qEAM FF to study spall failure,  crack propagation, and dislocation properties such as core structure and energy, Peierls stress, and kink formation energies  with systems containing as many as 50,000 atoms. We used such calculations with a micromechanical model of plasticity developed by Stainer, Cuitiño and Ortiz  to develop a multiscale model of single crystal plasticity in Ta. We are currently performing shock simulations in systems containing more than a million atoms using a parallel MD code; for a system containing 1,098,500 atoms 1 ps of MD simulations (1000 steps) takes approximately one hour on an SGI Origin 2000 computer using 128 250-MHz R10000 processors.

###### Acknowledgements.

This research was funded by a grant from DOE-ASCI-ASAP (DOE W-7405-ENG-48). The facilities of the MSC are also supported by grants from NSF (MRI CHE 99), ARO (MURI), ARO (DURIP), NASA, Kellogg, Dow Chemical, Seiko Epson, Avery Dennison, Chevron Corp., Asahi Chemical, 3M, GM, and Beckman Institute. Thanks to D. Singh and H. Krakauer for use of their LAPW code. We also thank Per Söderlind and John Moriarty for providing us with LMTO data. Computations were partially performed on the Cray SV1 at the Geophysical Laboratory, supported by NSF grant EAR-9975753 and the W. M. Keck Foundation.

FIGURE CAPTIONS

Figure 1. Zero temperature EOS for bcc Ta, LAPW GGA and qEAM FF results. Energy \[Figure 1(a)\] and pressure \[Figure 1(b)\] as a function of volume. Circles denote LAPW GGA results and lines show qEAM FF results.

Figure 2. Zero temperature EOS for fcc Ta, LAPW GGA and qEAM FF results. Energy \[Figure 2(a)\] and pressure \[Figure 2(b)\] as a function of volume. Circles denote LAPW GGA results and lines show qEAM FF results. The inset of Figure 2(a) shows fcc-bcc energy difference as a function of volume (circles joined by dots denote LAPW GGA results and line qEAM FF).

Figure 3. Zero temperature EOS for A15 Ta, ab initio and qEAM FF results. Energy \[Figure 3(a)\] and pressure \[Figure 3(b)\] as a function of volume. Open circles with dotted line denote full potential muffin-tin orbital calculations results Ref. ; filled circles show the sum of A15-bcc energy difference from  and our calculation of bcc energy using LAPW GGA method (section II). The full line shows qEAM FF results. The inset of Figure 3(a) shows A15-bcc energy difference as a function of volume (circles with dotted line denote QM results  and the line qEAM FF).

Figure 4. Zero temperature EOS for hcp Ta, LAPW GGA and qEAM FF results. Energy \[Figure 4(a)\] and pressure \[Figure 4(b)\] as a function of volume. Circles denote LAPW GGA results, section II, and lines show qEAM FF results. The inset of Figure 4(a) shows hcp-bcc energy difference as a function of volume (circles with dotted line denote LAPW GGA results and line qEAM FF).

Figure 5. Zero temperature elastic constants for Ta, LAPW GGA and qEAM FF results. Circles show bulk modulus \[$`(c_{11}+2c_{12})/3`$\]; diamonds show $`c_{44}`$ and squares represent $`c_s=(c_{11}c_{12})/2`$. qEAM FF results are shown with filled symbols and full lines and ab initio LAPW results with open symbols and dashed lines.

Figure 6. Volume relaxed vacancy formation energy (a) and enthalpy (b) as function of pressure. Dashed lines represent ab initio MB-PS results (section II) and full lines show qEAM results.

Figure 7. Vacancy formation volume as function of pressure. Dashed lines represent ab initio MB-PS results (section II) and full lines show qEAM results.

Figure 8. Vacancy formation enthalpy as function of pressure using the qEAM FF. Solid line shows the volume relaxed result; the dotted line the T=0 K fully relaxed results; dashed line is T=$`300`$ K result.

Figure 9. Vacancy migration energy using qEAM FF. Energy as a function of position of the marching atom at T=0 K and zero pressure. The vacancy migration energy is 1.093 eV.

Figure 10. Ideal shear strength of Ta using qEAM FF at zero temperature and volume V=18.36 $`\AA ^3`$. We show energy $`W(s)`$ \[Figure 10 (a)\] and stress $`\tau (s)`$ \[Figure 10 (b)\] as a function of shear.

Figure 11. Thermal expansion in Ta. (a) thermal expansivity as a function of temperature; circles represent qEAM FF results, the dashed-dotted lines shows mixed basis pseudopotantial calculations using the particle in a cell method, and the line denotes experimental results from Ref. . (b) Linear thermal expansion of Ta. $`(aa_0)/a_0`$ as a function of temperature; circles represents qEAM FF results and line denote experimental results from Ref. . The high temperature experimental data (dashed lines) are provisional values.

Figure 12. Melting of tantalum using the qEAM FF. Enthalpy (a) and volume (b) as a function of temperature for zero pressure; heating of bcc Ta (lower branches) and cooling of liquid Ta (higher branches). Heating and cooling rates are 100 K per 25 ps.

Figure 13. Ta melting using the qEAM FF. Two phase simulations. Time evolution of the potential energy in TPN MD at zero pressure for different temperatures. For T=3100 K (dotted line) the potential energy decreases with time, i.e. the system is crystallizing; for T=3200 K (dashed line) the potential energy grows denoting melting. For $`T=3150T_{melt}(P=0)`$ the potential energy remains constant; in this case the solid and liquid phases coexist in equilibrium.

Figure 14. Ta Melting qEAM FF. Melting curve for Ta up to P$`=300`$ GPa obtained using the two-phase simulation technique. We also show the experimental zero pressure melting temperature  (circles) and the results of shock melting.