Properties of Liquid Iron along the Melting Line up to the Earth-core Pressures

0.15 MB·English

by  Yu. D. Fomin

Properties of Liquid Iron along the Melting Line up to the Earth-core Pressures
Yu. D. Fomin
Institute for High Pressure Physics, Russian Academy of Sciences, Troitsk 142190, Moscow, Russia
V. N. Ryzhov and V. V. Brazhkin
Institute for High Pressure Physics, Russian Academy of Sciences, Troitsk 142190, Moscow, Russia and
Moscow Institute of Physics and Technology, 141700 Moscow, Russia
(Dated: January 31, 2013)
3 Wereportamoleculardynamicsstudyoftransportcoefficientsandinfinitefrequencyshearmod-
1 ulus of liquid iron at high temperatures and high pressures. We observe a simultaneous rise of
0
both shear viscosity and diffusion coefficient along the melting line and estimate if liquid iron can
2
vitrify underEarth-core conditions. Weshow that in frames of themodel studied in our work iron

n demonstratesamoderateincreaseofviscosityalongthemeltingline. Itisalsodemonstratedthatin
a thelimitofhightemperaturesandhighpressurestheliquidironbehavessimilartothesoftspheres
J system with exponent n≈4.6.

0
3 PACSnumbers: 61.20.Gy, 61.20.Ne,64.60.Kw

]
i I. INTRODUCTION idea is to employ some model of iron in molecular dy-
c
s namics.
-
rl The behavior of iron at the Earth core pressure- Several authors carried out MD simulations of iron at
mt temperatureconditionsisatopicofhotdebatesformany Earth-core conditions with different empirical potentials
decades. Apparently the behavior of iron is of great im- or by means of DFT method. In Refs. [4, 5] the authors
.
t portance for understanding of the phenomena occurring use the same parametrization of embedded atom poten-
a
m in the inner and outer core. However, it appears to be tial (EAM) for iron proposed by Sutton and Chen [6].
verydifficulttofindanyunambiguousinformationabout TheauthorsofRef. [4]andRef. [5]havechosendifferent
-
d ironatsuchextremeconditions. Theproblemisclear: it density-temperaturepointswhichmakesmoredifficultto
n is impossible to carry out direct experiments at so high comparetheirs results. However,the viscositydata from
o temperatures. Therefore one has to use extrapolation both articles look to be consistent with each other. The
c
[ of the lower temperatures results. This leads to a great viscosities at the Earth core conditions obtained in both
variationofthemostprincipledata. Eventhelocationof cited papers are of the orderof magnitude of0.01 Pa·s.
1 themeltinglineofironathighpressuresisnotclear: the In Ref. [5] the viscosities are also compared to the ab-
v
8 results from diamondanvil experiments give the melting initio calculations of Alfe et. al [7] (Table I of Ref. [7]).
5 temperature approximatelytwice smaller then the shock Onecanseefromthiscomparisonthatthedatafromclas-
1 waveexperiments(apossibleexplanationofthisdiscrep- sicalMDcalculationsandfromab-initioMDarecloseto
7 ancy was proposed in Ref. [1]). each other and that there is no systematic deviation of
.
1 Thesituationwithtransportpropertiesofironathigh EAM data from the ab-initio ones.
0 pressureis evenmuchworse. The difference betweenthe
3 viscosities obtained by different methods achieves 1014 A lot of information about iron at high pressure was
1 recently obtained by ab-initio simulations. It includes
[2]. Seccoclassifiedtheironviscosityestimationsinthree
v: groups [3]: the ones from geodesic measurements and the melting line calculations [8, 9], transport properties
i of iron [7] and mixtures of iron with other elements [10].
X seismologicalinvestigations give the Earth-coreviscosity
upto1011Pa·s. TheviscositiesobtainedfromtheEarth However, ab-initio simulations are very computationally
r
a magnetic field are of the order of 2.7·107 Pa·s. Finally expensiveanddonotallowtostudythepropertiesofiron
in a vast region of pressure - temperature points.
the theoretical predictions give the iron viscosities from
2.5·10−3 up to 50 Pa·s [3]. Obviously the discrepancy The goal of the present study is to perform a system-
of 1014 between different methods can not be recognized aticstudyofliquidironpropertiesalongthe meltingline
as acceptable. in a wide rangeof temperatures andpressures. This will
One of the possible sources of errors in the high pres- allowustoseethechangeswhichtakeplaceinamodelof
sure iron viscosity estimation can originate from the ex- realisticliquidunderhugechangeinpressureandtemper-
trapolation of low pressure data far beyond the range of ature. Since we use the same model and all data points
pressures where we have experimental data. In this re- are obtained directly we do not use any extrapolation
spectitisreasonabletofindasystemwhichwecanstudy procedures which allows us to see the general trends in
with reasonable precision in wide range of pressures and the liquid iron behavior along the melting line up to the
temperatures. The most obvious way to implement this very high pressures.
2
II. SYSTEMS AND METHODS All of these distributions written to the file were used to
estimate the viscosity.
All simulations reported in this work were done by
As it was stated above the goal of this article is to
LAMMPS simulation package [13].
studythepropertiesofliquidironalongthemeltingline.
In the literature there is a lot of different data on the
melting line of iron obtained by different groups and us-
III. RESULTS AND DISCUSSION
ing different methods. Several authors reported melting
line of iron from molecular dynamic simulation.
As it was mentioned in the introduction we are inter-
The most extensive simulation of melting line of iron
esting in the behavior of the transportcoefficients of liq-
was performed by Belonosho et. al [11]. Basing on the
uid iron along the melting line. In our previous publica-
EAM potential for iron introduced in this work the au-
tionsweanalyzedthebehaviorofsimpleliquidsalongthe
thors computed the melting line from P = 60 Gpa up
melting curve [14, 15]. Two simple models were studied
to 400 GPa. This corresponds to the temperatures up
- soft spheres (Φ(r) =ε(σ)n with n=12) and Lennard-
to 8000 K which is even above the estimated Earth core Jones (Φ(r) = ε·((σ)12r−(σ)6)) liquids. It was shown
temperature. r r
that in these models both shear viscosity and diffusion
Basing on these well documented data from Ref. [11] coefficient grow up upon increasing pressure along the
wecomputethepropertiesofironalongthe meltingline. melting line. However, even if these model liquids can
For doing this we simulate a system of 3456 iron atoms be used as a rough model to analyze the most general
in the (P,T) points located along the melting line from trends in liquids, the simplicity of these models can re-
Tlow = 2500 K up to Thigh = 8000 K with the tempera- sultinlargequalitativeerrorsathighpressuresandhigh
ture step dT =100 K. Firstly we simulate the system in temperatures comparing to the experimental liquids. In
NPT ensemble to find the equilibrium density of liquid this respect the present work is aimed to understanding
at the given temperature and pressure. At this stage we the high temperature - high pressure behavior of an ex-
carry out 106 MD steps with the step size dt=0.001ps. ampleofrealliquid. Wechooseaparticularcaseofliquid
When the equilibrium density is found we simulate the iron due to its importance in geophysicalinvestigations.
system in NVT ensemble at this density in order to cal- Fig.1(a)showsthe meltinglineofiron[11]andFig.1
culate the equilibrium structure and infinite frequency (b) demonstrates the liquid branchofthe melting line in
shear modulus of the liquid. At this stage the system the density - temperature plane. This line was obtained
is propagated 106 time steps with dt = 0.0002 ps. Fi- by performing NPT simulations at the PT data points
nallyusingtheobtainedfinalstructureatthechosenden- fromRef. [11]. Asitisexpectedthedensityofthe liquid
sity and temperature as initial conditions we carry out quickly rises along the melting curve.
NVE simulationofthe sampletofindthe diffusioncoef- Figs. 2 (a) and (b) represent the diffusion coefficient
ficient. Fordoingthiswesimulatethesystemfor2.6·106 and shear viscosity of liquid iron along the melting line.
steps with dt = 0.0002 ps. The infinite frequency shear As in the case of simple models both the diffusion coeffi-
modulus is evaluated as Ginf = βV < Px2y >, where cient and the viscosity rapidly increase with increasing
β = 1/(kBT) and Pxy is off-diagonal pressure compo- the temperature. The viscosity rise is especially dra-
nent. Thediffusioncoefficientiscomputedfromtheslope matic: at the highest temperature studied (8000 K) it
ofmeansquaredisplacementoftheparticlesviaEinstein is 2.5 times higher then at the lowest one (2500 K). At
relation. the same time the diffusion coefficient increases just 1.5
Oneofthecentralquantitiesofouranalysisistheshear times. Anyway the liquid becomes more viscous and
viscosity of liquid iron along the melting line. In order more diffusive at the same time.
to compute the viscosity we employ the Reverse Non- The question of simultaneous growing up of diffusion
equilibrium MD method also known as Mu¨ller-Plathe coefficient and shear viscosity was raised up in our pre-
method[12]. Inthismethodanartificialmomentumflux vious publication [15]. This question is interesting in
is imposed in the system which results in a linear veloc- conjunctionwithglasstransition. Themostcommoncri-
ity profile of particles. The viscosity coefficient can be teriaofglasstransitionstatesthattheliquidexperiences
calculated from the slope of the velocity profile and the glasstransitionwhenitsviscosityreachessomeveryhigh
momentumtransferredtothesystem[12]. Fortheviscos- value. A typical convention is that the viscosity of glass
ity calculation the system was simulated for 2·106 time transitionis1013Poise. However,itisimplicitlyassumed
steps with dt = 0.0002 ps. The momentum transfer was that the liquid looses its diffusivity under this viscosity
undertaken every 10 steps. The temperature was held grows. However, in case of high temperatures and high
by couplingto Berendsenthermostatwith time constant pressures both viscosity and diffusivity increase, so one
t = 10·dt. The first half of the simulation was used comes to a contradiction to the usual common view on
B
for equilibration while during the second half the veloc- glass transition.
ity distribution was written in the file every 100 steps. InordertosolvethiscontradictionweproposedinRef.
3
(a)
0.65 (a)
8000
0.60
7000
s
0.55
/
6000 2
m
K
0.50
T,5000 -8 0
1
4000 , 0.45
D
3000 0.40
2000 0.35
0 100 200 300 400 2000 3000 4000 5000 6000 7000 8000
P,GPa T,K
(b) (b)
16
8000
14
7000
12
s
6000
*
K a
P 10
T,5000 m
,
8
4000
3000 6
2000 4
9 10 11 12 13 14 2000 3000 4000 5000 6000 7000 8000
3
,g/cm T,K
FIG.1: (Coloronline)(a)IronmeltinglineasobtainedinRef. FIG. 2: (Color online) Diffusion (a) and shear viscosity (b)
[11] in P −T oordinates. (b) Liquid branch of the melting along the melting line. Squares - MD data, continuous line -
line. Squares - MD data, continuous line - soft spheres-like soft spheres approximation (see the text).
approximation (see thetext);
η
[15]touseonemorecriteriumofglasstransition: theliq- τ = , (1)
G
uid undergoes the glass transition if the relaxation time inf
becomes as long as the typical time of experiment. Dif-
whereη istheviscosityoftheliquidandG theinfinite
inf
ferentpublicationsproposetouseor100secondsor1000
frequency shear modulus. In the majority of experimen-
secondsastheglasstransitionrelaxationtime. Following
tal situations it is supposed that the viscosity changes
this definitionone needsto seethe behaviorofthe relax-
much faster then G and the relaxationtime is mainly
inf
ation time in order to understand if the liquid vitrifies
determined by the viscosity behavior. This case the vis-
following the melting line up to extremely high temper-
cositycriteriumofglasstransitionbecomesequivalentto
atures - high pressures limit.
the relaxation time one. However, as it was shown in
The relaxationtime can be computed via Maxwell re- Ref. [15] the infinite frequency shear modulus of liquid
lation [16] can dramatically grow along the melting line.
4
soft spheres model is especially simple since it demon-
strates a set of scaling properties along the melting line
[17–20]. Thescalingrelationsforallquantitiespresented
25 (a)
here are given in Refs. [14, 15]. Here we repeat them for
the sake of completeness:
20
P ∼T1+3/n, (2)
a
P 15
G

,nf10 D ∼T1/2−1/n, (3)
Gi
5
η ∼T1/2+2/n, (4)
0
2000 3000 4000 5000 6000 7000 8000
G ∼T1+3/n. (5)
inf
T,K
In order to see the relations between the simplest
model studied and the current liquid iron system we fit
(b) all the quantities above (D, η, P and G ) to the rela-
inf
2.0 tions of the form X = a·Tα +b, X is the quantity of
interest and α is the correspondent exponent from eqs.
1.8 (2) - (5). In order to get all of the exponents consistent
with the case of soft spheres all quantities were fitted
1.6 simultaneously. The results of such fitting are given in
Figs. 1 - 3. The exponent coefficient n is found to be
s
1.4
p equal to n=4.568. One can see that except the melting
, pressure all of the quantities are well represented by the
1.2
softspheres-likescalingrelations. Inthe caseof pressure
the deviation does not exceed 12% in the whole range
1.0
of temperatures considered in this work. However, the
0.8 slope dP fromthescalingformulaandfromMDdataare
dT
very different which means that the melting line itself is
0.6 poorly represented by the scaling law. At the same time
2000 3000 4000 5000 6000 7000 8000
the transport coefficients and elastic properties (G )
inf
are well described by the soft spheres-like model.
T,K
It is well known that the structure of liquid metals
FIG. 3: (Color online) Infinite frequency shear modulus (a) can be well approximated by simple hard spheres model
andrelaxationtime(b)alongthemeltingline. Squares-MD
[21]. In Ref. [22] experimental measurements of liquid
data, continuous line - soft spheres approximation (see the
ironstructuresfactorswerereportedandthecomparison
text).
of experimental curves with the hard spheres model was
done. As it followsfrom Fig. 3 of Ref. [22] the structure
Fig.3(a)showstheinfinitefrequencyshearmodulusof factorsofironcanbesufficientlywellrepresentedbyhard
iron along the melting line. One can see that G dras- spheres ones which proves that the main contribution
inf
tically increases with increasing the temperature. The intotheliquidstructurecomesfromrepulsivepartofthe
ratio of G at the highest and the lowesttemperatures interaction. Itiswellknownthatthestructureofliquidis
inf
isapproximately10,whilefortheviscosityitis2.5. Asa closely related to its transport properties and therefore
result in spite of the rise of viscosity the relaxation time one can expect that simple purely repulsive models of
still decreases (Fig. 3 (b)). liquid can reproduce the diffusion and shear viscosity of
iron sufficiently well.
From the results presented above one can see that the At the same time melting line is strongly affected by
qualitativebehaviorofliquidironathightemperatures- thepresenceofattractivetermsintheinterparticleinter-
highpressuresisequivalenttothe behaviorofsimple liq- actionpotential[23]whichmeansthatapurelyrepulsive
uids such as soft spheres and Lennard-Jones ones. The model such as soft spheres should fail to reproduce the
5
meltingcurveofasystemwithbothrepulsiveandattrac- tionsconfirmthisspeculationandproposethattheexact
tive interactions. resultsforsoftspheresreportedinourpreviouswork[15]
One of the experimental studies of the liquid iron vis- canbe extrapolatedto the hightemperature-highpres-
cosityathighpressurewasreportedinRef. [24]. Theau- sure limit of liquids in general.
thorsofthispapermeasuredtheviscosityattemperature
Y.F is grateful to A. B. Belonoshko (Theoretical
ashighas 2050K andfound thatthe changein viscosity
Physics KTH, Sweden) for sharing his results for the
comparing to the room temperature is small. Basing on
ironmelting line and V.V. Stegailov(JIHT RAS) for his
this result they concluded that the statement proposed
help with simulations. Y.F. also thanks the JointSuper-
by Poirier that the viscosity of liquid is nearly constant
computing Center of the Russian Academy of Sciences
along the melting line [25] is correct. From our results
andtheRussianScientificCenterKurchatovInstitutefor
wecanconcludethatitisjustpartiallytrue. Theviscos-
computationalfacilities. Theworkwassupportedinpart
ity change along the melting line is not fast: it increases
by the Russian Foundation for Basic Research (Grants
2.5 times on 3.2-fold temperature change. However, we
No 13-02-00579, No 13-02-00913, No 11-02-00341-a and
observe the systematic rise of viscosity along the melt-
No 11-02-00303)and the Ministry of Education and Sci-
ing line so one can not claim that it is constant: if on
ence of Russian Federation, projects 8370 and 8512.
measures the viscosities along the melting line for large
enough temperature interval one will clearly see the rise
of viscosity. However, the temperatures studied in our
workrangefrom2500K up to 8000K whichexceedsthe
range of temperatures reported in most of experimental [1] S.V. Starikov and V. V. Stegailov, Phys. Rev. B, 80,
works. It means that in the range of temperatures ex- 220104 (2009).
plored in experiments the viscosity change can be small [2] D.E. Smylie, V.V. Brazhkin, A. Palmer, Uspekhi Fiz-
enoughto use Poirierstatement with sufficientaccuracy. icgeskikh Naus 179, 91 (2009).
[3] R.A. Secco ”Viscosity of the outer core”, in Mineral
Physics and Crystallography: A Handbook of Physi-
CONCLUSIONS cal Constans (AGU Refer- ence Shelf, Vol. 2, Ed. T J
Ahrens)(Washington,DC:AmericanGeophysicalUnion,
1995) p.218.
The present article represents a molecular dynamics [4] Y.Zhang,G.Guo,G.Nie,Phys.Chem.Minerals27,164
study oftransportcoefficients andglasstransitionofliq- (2000).
uidironinthelimitofhightemperatures-highpressures [5] C. Desgranges and J. Delhommelle, Phys. Rev. B 76,
along the melting line. We show that both shear viscos- 172102 (2007).
[6] A. P. Sutton and S. Chen, Philos. Mag. Lett. 61, 139
ity and diffusion coefficient increase along the melting
(1990).
line. However, due to very rapid increase of the infi- [7] D. Alfe, G. Kresse, and M. J. Gillan, Phys. Rev. B 61,
nite frequency shear modulus the relaxation time drops
132 (2000).
quickly with increasing the pressure along the melting [8] D. Alfe, G.D. Price and M.J. Gillan, Phys. Rev. B 65,
curve which means that liquid iron becomes harder to 165118 (2002).
vitrify at higher temperatures and higher pressures. [9] D. Alfe, ”Iron at Earth’s Core Conditions from First
Principles Calculations”, Book Series: Reviews in Min-
Itisworthtonotethatthemagnitudeofviscositieswe
earology & Geochemistry 107, 337-354 (2010).
obtainareconsistentwithothersimulationsofliquidiron
[10] M. Pozzo, Ch.Davies, D. Gubbins and D. Alfe, Phys.
at high temperatures and high pressures [4, 5, 7]. How- Rev. B 87, 014110 (2013).
ever,it looksthatallsimulationsofironatsuchextreme [11] A. B. Belonoshko, R. Ahuja, and B. Johansson, Phys.
conditions strongly underestimate the shear viscosity. Rev. Lett.84, 3638 (2000).
Surprisingly, the behavior of liquid iron at Earth-core [12] Muller-Plathe, PhysRev E, 59, 4894 (1999).
like temperatures and pressures can be sufficiently well [13] S. Plimpton, Fast Parallel Algorithms for Short-Range
Molecular Dynamics, J. Comp. Phys., 117, 1 (1995);
qualitativelydescribedbysoftspheresmodelwhichisone
http://lammps.sandia.gov/
of the simplest models of liquid. By fitting the MD data
[14] Yu.D.Fomin,V.V.Brazhkin,V.N.RyzhovJEPTLett.,
to the soft spheres scaling relation we find that liquid 95, 320 (2012).
iron is qualitatively similar to the soft spheres with the [15] Yu.D.Fomin,V.V.Brazhkin,V.N.RyzhovPhys.Rev.
exponent n=4.568. E 86, 011503 (2012).
The most important difference between the soft [16] J.-P.Hansen,I.R.McDonald.”Theoryofsimpleliquids”,
spheres and liquid iron represented by the EAM poten- Academic Press (1986).
[17] O. Klein, Medd. Vetenskapsakad. Nobelinst. 5, No. 6
tial Ref. [11] is in the collective nature of the later. It
(1919).
is well knownthat in the limit of high pressuresthe par-
[18] T.H. Berlin and E.W. Montroll, J. Chem. Phys. 20, 75
ticles come very close to each other and the system is
(1952).
dominated by the repulsive excluded volume effects and [19] V. V. Zhakhovsky, Zh. Eksp. Teor. Fiz. 105, 1615
the collective effects can become negligible. Our simula- (1994)[JETP 105, 1615 (1994)].
6
[20] W.G. Hoover, M. Ross and K. W. Jonson, Journal of Phys., J. Chem. Phys.134, 044523 (2011).
Chemical Physics. 52, 4931 (1970). [24] M. D. Rutter, R. A. Secco, H. Liu, T. Uchida, M. L.
[21] N.W.AshcroftandJ.Lekner,Phys.Rev.145,83(1966). Rivers, St. R. Sutton, and Y. Wang, Phys. Rev. B, 66,
[22] G.Shen,V.B. Prakapenka,M.L. RiversandSt.R.Sut- 060102(R) (2002).
ton, Phys.Rev.Lett. 92, 185701 (2004). [25] J.P. Poirier, Geophys. J. 92, 99 (1988).
[23] Yu.D.Fomin,E.N.Tsiok, and V.N.Ryzhov,J. Chem.