Didier Fraix-Burnet
email: fraix@obs.ujf-grenoble.fr

MULTIVARIATE EVOLUTIONARY ANALYSES IN ASTROPHYSICS

come   L'archive ouverte pluridisciplinaire

Introduction

We are now able to study galaxies in great detail, identifying individual stars, gas and dust clouds, as well as different stellar populations. Imagery brings very fine structural details, and spectroscopy provides the kinematical, physical and chemical conditions of the observed entities at different locations within a galaxy. For more distant objects, information is scarcer, but deep systematic sky surveys gather spectra for millions of galaxies at various redshifts. The amount of data on galaxies, their number, their diversity, their complexity and that of their evolution, suggest that they should be envisaged as a population or an ensemble of populations. This implies the use of the appropriate statistical tools.

Like paleontologists, we observe objects from the distant past (galaxies at high redshift), and like evolutionary biologists, we want to understand their relationships with nearby galaxies, like our own Milky Way. Consequently, a "galactogenesis" can be advantageously approached by considering phylogenetics methods.

2 Why be multivariate?

The description of a given galaxy requires many observables, most of them derived from a spectrum. Usual classifications, often inspired by the Hubble tuning fork, use only a very few properties. Even if bivariate plots or correlations are clear and useful, they are incomplete. Worse, they are merely the projection onto a 2-D diagram of a multivariate parameter space. This projection is generally expected to increase the dispersion of the plot. Anyhow, it is difficult to represent many data with only bivariate plots, and any classification necessarily requires an arbitrary binning of one or several parameters.

Multivariate analyses are still not much used in astrophysics. One basic tool, the Principal Component Analysis, is relatively well-known (e.g. Cabanac et al., 2002;[START_REF] Recio-Blanco | Multivariate analysis of globular cluster horizontal branch morphology: searching for the second parameter[END_REF], but this is not a clustering tool in itself. A very few attempts to apply multivariate clustering methods have been made very recently (Chattopadhyay andChattopadhyay, 2006, 2007;Chattopadhyay et al., 2007[START_REF] Chattopadhyay | Globular clusters in the milky way and dwarf galaxies -a distribution-free statistical comparison[END_REF]Chattopadhyay et al., , 2009a,b;,b;[START_REF] Fraix-Burnet | The environment of formation as a second parameter for globular cluster classification[END_REF][START_REF] Fraix-Burnet | Structures in the fundamental plane of early-type galaxies[END_REF]. Sophisticated statistical tools are used in some areas of astrophysics and are developing steadily, but multivariate analysis and clustering techniques have not much penetrated the community. It is true that the interpretation of the results are not always easy.

3 Why be evolutive?

Evolution, an unavoidable fact, is also not correctly taken into account in most classification methods. By mixing together objects at different stages of evolution, most of the physical significance and usefulness of a classification is lost. In practice, the evolution of galaxies is often limited to the evolution of the properties of the entire population as a function of redshift [START_REF] Bell | Galaxy Assembly. Planets to Cosmology: Essential Science in Hubble's Final Years[END_REF]. Since environment (the expanding Universe) and galaxy properties are so much intricate, this kind of study is relevant to a first approximation. However, recent observations have revealed that galaxies of all kinds do not evolve perfectly in parallel, as illustrated for instance by the so-called downsizing effect which shows that large galaxies formed their stars earlier than small ones (e.g. [START_REF] Neistein | Natural downsizing in hierarchical galaxy formation[END_REF]. New observational instruments now bring multivariate information at different stages of evolution, and in various evolutive environments. In this multivariate context, we believe that the notion of "evolution", easy to understand for a single parameter, is advantageously replaced by "diversification".

The transformation of galaxies is a complex process (Fraix-Burnet et al., 2006b,c) that cannot be disentangled with only a very few observables. For instance, the elliptical shape of galaxy can be obtained through the monolithic collapse of a big cloud of gas, or by big mergers. To find which process has shaped a given galaxy, many observables are required. Only a multivariate and evolutive analysis can distinguish different histories.

Classification, complexity, evolution

Multivariate clustering methods compare objects with a given measure and then gather them according to a proximity criterion. There are two main classes. Firstly, distance analyses are based on the overall similarity derived from the values of the parameters describing the objects. The choice of the most adequate distance measure for the data under study is not unique and remains difficult to justify a priori. The way objects are subsequently grouped together is also not uniquely defined. Secondly, methods based on characters (a trait, a descriptor, an observable, or a property, that can be given at least two states characterizing the evolutionary stages of the object for that character) compare objects in their evolutionary relationships [START_REF] Wiley | The Compleat Cladist: A Primer of Phylogenetic Procedures[END_REF]. Here, the "distance" is an evolutionary cost simply measured by the number of changes of the parameter values (or character states). Groupings are then made on the basis of shared or inherited characteristics, and are most conveniently represented on an evolutionary tree.

Character-based methods like cladistics are better suited to the study of complex objects in evolution, even though the relative evolutionary costs of the different characters is not easy to assess. Distance-based methods are generally faster and often produce comparable results, but the overall similarity is not always adequate to compare evolving objects. In any case, one has to choose a multivariate method, and the results are generally somewhat different depending on this choice [START_REF] Buchanan | Phenetics, cladistics, and the search for the alaskan ancestors of the paleoindians: a reassessment of relationships among the clovis, nenana, and denali archaeological complexes[END_REF]. However, the main goal is to reveal a hidden structure in the data sample, and the relevance of the method is mainly provided by the interpretation and usefulness of the result.

We must note that taking all available parameters blindly can kill the multivariate and evolutive analysis. One dangerous component is a hidden correlation, such as a size effect, that creates a redundancy. A less known caveat is due to spurious correlations, due to independant variables that vary as function of a non-necessarily obvious parameter. This is especially the case with the time or the stage of evolution. Two quantities can be totally unrelated but if they vary both with time in a more or less monotonic way, then they appear to be correlated. For instance, all photometric quantities for galaxies are affected by the stellar evolution. In such a case, a cladistic analysis yields a regular tree showing the regular stellar evolution (e.g. [START_REF] Fraix-Burnet | Determining the evolutionary history of galaxies by astrocladistics: some results on close galaxies[END_REF].

Multivariate evolutionary classification in astrophysics has been pioneered by the author (Fraix-Burnet et al., 2006a[START_REF] Fraix-Burnet | The environment of formation as a second parameter for globular cluster classification[END_REF][START_REF] Fraix-Burnet | Structures in the fundamental plane of early-type galaxies[END_REF][START_REF] Fraix-Burnet | Evolutionary Biology Concept, Modeling, and Application[END_REF]. Called astrocladistics, it is based on cladistics that is heavily developed in evolutionary biology. Astrocladistics has been first applied to galaxies (Fraix-Burnet et al., 2006b[START_REF] Fraix-Burnet | Structures in the fundamental plane of early-type galaxies[END_REF] because they can be shown to follow a transmission with modification process when they are transformed through assembling, internal evolution, interaction, merger or stripping. For each transformation event, stars, gas and dust are transmitted to the new object with some modification of their properties. Cladistics has also been applied to globular clusters [START_REF] Fraix-Burnet | The environment of formation as a second parameter for globular cluster classification[END_REF], where interactions and mergers are probably rare. These are thus simpler stellar systems, even though we have firm evidence that internal evolution can create another generation of stars and that globular clusters can lose mass. Basically, the properties of a globular cluster strongly depend on the environment in which it formed (chemical composition and dynamics), and also on the internal evolution which includes at least the aging of its stellar populations. Since galaxies and globular clusters form in a very evolving environment (Universe, dark matter haloes, galaxy clusters, chemical and dynamical environment), the basic properties of different objects are related to each other by some evolutionary pattern.

A more pertinent physical interpretation

An obvious difficulty for a physicist in general is to intepret the results of multivariate analyses using his models that mostly result from a set of equations and are more conveniently presented by curves on bivariate plots. Interestingly enough, these models are multivariate, especially in astrophysics, and the resolution of the set of equations yields a "population" of possible results often called a grid of models. As a result, some parameters are set to sensible values, and the corresponding models are then compared to some observables. These observables can also have been truncated by setting some other observables in order to simplify the information.

It appears that we must here compare two populations, one of real objects and one of models, in a multivariate space. We show here two examples of multivariate (and evolutive) analyses of astrophysical objects showing that such approaches are both more direct, objective and physically pertinent.

Figure 1 shows the cladogram obtained for globular clusters of our Galaxy and the projection of the partitioning on pair plots for the four parameters used for the analysis: logTe, that measures the temperature of stars that are at a specific point in their evolution, Fe/H, MV that is the total visible intensity (magnitude) and roughly indicates the mass of the globular cluster, and Age that can be measured quite precisely because all stars of a given globular cluster are formed nearly at the same time. However Age is not an intrinsic property discriminating evolutionary groups since it evolves in the same way for all. But we gave it a half weight to arrange the objects within each group [START_REF] Fraix-Burnet | The environment of formation as a second parameter for globular cluster classification[END_REF]. Three groups are identified. The first one (in blue) has on average the lower ratio Fe/H that measures the proportion of heavy atomic elements that are processed within stars. This group is consequently considered as more primitive. It is obivous that this partitioning would be impossible to obtain with only bivariate plots.

Looking at other parameters (such as orbital elements, kinematics, more refined chemical abundances...) revealed clear characteristics that allowed us to infer that each group formed during a particular stage of the assembly history of our Galaxy. The blue group is the older one. It formed during the dissipationless collapse of the protogalaxy. They are located mainly in the outer halo. The red group belongs to the inner halo and the corre- sponding clusters formed at a later stage during the dissipational phase of Galactic collapse, which continued in the halo after the formation of the thick disc and its globular clusters. These clusters were very massive before "star evaporation" took place. The latter group (green) formed during an intermediate and relatively short period and comprises clusters of the disk of our Galaxy (all details in [START_REF] Fraix-Burnet | The environment of formation as a second parameter for globular cluster classification[END_REF].

Another example is given with the fundamental plane of early-type galaxies which is a lonog-known correlation between the central velocity dispersion σ, the surface brightness µ e and the effective radius r e . In addition, the metallicity, as measured with the M g 2 index, plays a role and seems to be correlated with log σ. We performed a K-means cluster analysis and a cladistic analysis in parallel [START_REF] Fraix-Burnet | Structures in the fundamental plane of early-type galaxies[END_REF]. The partionings are remarkably in agreement. We believe the reason is due to the careful choice of the parameters. For cladistics, they must be informative with respect to diversification, and should not be redundant or incompatible. This requirement is logically pertinent also for any cluster analysis.

Cladistics provides in addition the evolutionary relationships between the groups. On Figure 2, the tree is projected onto the log σ vs µ e diagram on which the fundamental plane is seen essentialy face-on. Since galaxies are more complicated than globular clusters, the interpretation of the results and all relations between all possible parameters and within each group takes great advantage of numerical simulations. Here again, we are able to derive the probable history of each group of galaxies as well as their relative level of diversification, giving possible sequences of past transforming events such as mergers, accretions or sweeping (for details, see [START_REF] Fraix-Burnet | Structures in the fundamental plane of early-type galaxies[END_REF].

A quite interesting finding is that most known correlations are different or even absent when we consider groups individually (Figure 3). This proves that they have different evolution histories. Another noticeable fact is that the well-known correlation between Mg2 and log σ is indeed spurious, or historical. It is due to the fact that each parameter changes with the level of diversification as clearly shown by the placement of the groups (see Figure 3).

Conclusion

Undoubtly, the study of galaxies now requires multivariate statistical treatments. Evolution must also be taken into account and the concept of populations seems appropriate and points to the use of methodologies developed elsewhere. Complexity, evolution and classification suggest similar studies as in phylogenetics. Astrocladistics has opened the pathway.

Figure 1 :

 1 Figure 1: Top: cladogram obtained for the globular clusters of our Galaxy. Bottom: projection of the partitioning on pair plots with the four parameters used for the cladistic analysis. From Fraix-Burnet et al. (2009).

Figure 2 :

 2 Figure2: The fundamental plane of 699 galaxies, showing the partitioning and projection of the tree obtained by a cladistic analysis[START_REF] Fraix-Burnet | Structures in the fundamental plane of early-type galaxies[END_REF].

Figure 3 :

 3 Figure 3: Cluster and cladistic analysis of the fundamental plane of earlytype galaxies: bivariate plots showing how correlations differ for each group and for the whole sample. Note in particular the M g 2 vs log σ plot revealing a spurious correlation between these two parameters (Fraix-Burnet et al., 2010).


Ján Šimkanin

Convection in rotating non-uniformly stratified spherical fluid shells in dependence on Ekman and Prandtl numbers

Keywords: rotating convection, non-uniform stratification, penetrative convection, small-scale structures

We present an investigation of rotating convection in non-uniformly stratified spherical shells in dependence on Prandtl, Ekman and Rayleigh numbers. Attention is focused on the case in which the thickness of both sublayers (stable and unstable) is identical, which was not investigated before. Influence of stratification is more evident at small values of Rayleigh numbers (large values of Ekman numbers), at large Rayleigh numbers is almost negligible. For small Rayleigh numbers in the case of uniform stratification a large-scale columnar convection is developed in the whole volume, in the case of non-uniform stratification it is suppressed to the unstably stratified region but slightly penetrates to the stably stratified one, except the case of large Prandtl numbers when the convection significantly penetrate to the stably stratified region. The spiralling nature of convection is a dominant feature at moderate Prandtl numbers, the strong spiralling is typical at small Prandtl numbers but it disappears at large Prandtl numbers. For large Rayleigh numbers (small values of Ekm developed. The multilayer convection (``teleconvection") is not developed in our case of non-uniform stratification because of the significant amount of stable stratification. Our case is a typical extreme case. It is dissimilar to both cases if the stably stratified sublayer is thinner than the unstably stratified one and if the unstably stratified sublayer is thinner than the stably stratified one.

Introduction

The Earth's and planetary fluid interiors are characterized by convective motions. Convection (magnetoconvection) constitutes the driving mechanism of hydromagnetic processes leading to magnetic field generation [START_REF] Roberts | Geodynamo theory and simulations[END_REF]. Buoyancy, the fundamental source of convection or magnetoconvection [START_REF] Jones | Convection-driven geodynamo models[END_REF], results from the complicated processes going on in the Earth's and planetary fluid interiors, e.g., a chemical homogenisation, gravitational differentiation, solidification processes acting on the inner core boundary (e.g., the convection in the mushy layer due to the mentioned solidification processes, see [START_REF] Guba | Nonlinear oscillatory convection in mushy layers[END_REF]), etc. Consequently, the outer liquid Earth's core and the liquid interiors of Giant planets are nonuniformly stratified (this means density stratification, see [START_REF] Fearn | Compositional convection and stratification of the Earth's core[END_REF][START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF]).

In general, it is assumed that the upper part of the outer liquid Earth's core (close to the core-mantle boundary 1 ) is stably stratified (subadiabatic radial temperature gradient) and the lower part (towards the inner core boundary 2 ) unstably stratified (superadiabatic radial temperature gradient). The idea of "a stably stratified ocean" at the top of the outer Earth's core was firstly introduced by [START_REF] Braginsky | Magnetohydrodynamics of the Earth's Core[END_REF]. Models of a non-uniformly stratified fluid shell are an acceptable simplification of the real Earth-like conditions. In the Earth's core the stably stratified sublayer is probably very thin (the outer Earth's core is almost unstably stratified, see [START_REF] Braginsky | Magnetohydrodynamics of the Earth's Core[END_REF][START_REF] Fearn | Compositional convection and stratification of the Earth's core[END_REF][START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Šimkanin | Problem of the rotating magnetoconvection in variously stratified fluid layer revisited[END_REF][START_REF] Šimkanin | On hydromagnetic instabilities and the mean electromotive force in a non-uniformly stratified Earth's core affected by viscosity[END_REF][START_REF] Stanley | Effects of an outer thin stably stratified layer on planetary dynamos[END_REF]). Such a stratification is probably typical for the terrestrial planets [START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Stanley | Numerical dynamo models of Uranus' and Neptune's magnetic fields[END_REF][START_REF] Stanley | Effects of an outer thin stably stratified layer on planetary dynamos[END_REF]. However, in the other planets the ratio of the thickness of the appropriate sublayers (e.g., of the stably stratified to unstably stratified sublayers) and their geometric configuration vary. This is noticeable especially with the Giant planets [START_REF] Stanley | Convective-region geometry as the cause of Uranus' and Neptune's unusual magnetic fields[END_REF][START_REF] Stanley | Numerical dynamo models of Uranus' and Neptune's magnetic fields[END_REF][START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Stanley | Effects of an outer thin stably stratified layer on planetary dynamos[END_REF].

Non-uniform stratification can be simulated thermodynamically also in the Boussinesq models by means of internal heat sources [START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Šimkanin | Problem of the rotating magnetoconvection in variously stratified fluid layer revisited[END_REF][START_REF] Šimkanin | On hydromagnetic instabilities and the mean electromotive force in a non-uniformly stratified Earth's core affected by viscosity[END_REF]. If the stably stratified sublayer is very thin (for a geometric configuration stable/unstable), then behaviour is close to the case of uniform stratification when the whole layer is unstably stratified. Likewise, if the unstably stratified sublayer is very thin, then behaviour is close to the case of uniform stratification when the whole layer is stably stratified. Thus, the effects of non-uniform stratification are noticeable if the thicknesses of the stably and unstably stratified sublayers are comparable [START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Šimkanin | Problem of the rotating magnetoconvection in variously stratified fluid layer revisited[END_REF][START_REF] Šimkanin | On hydromagnetic instabilities and the mean electromotive force in a non-uniformly stratified Earth's core affected by viscosity[END_REF][START_REF] Šimkanin | Control volume method for hydromagnetic dynamos in non-uniformly stratified spherical shells[END_REF].

The study of hydromagnetic dynamo action in Mercury provides another good example of the influence of a stably stratified sublayer [START_REF] Christensen | A deep dynamo generating Mercury's magnetic field[END_REF]. Mercury is characterized by a weak magnetic field. A possible explanation could be given by a hydromagnetic dynamo working in the similar geometric configuration as in our study (stable/unstable) but in this case a larger fraction of spherical shell is stably stratified [START_REF] Christensen | A deep dynamo generating Mercury's magnetic field[END_REF]. In such a case the (magneto)convection and the dynamo action are strongly suppressed in the upper stably stratified sublayer, i.e. magneto(convection) and dynamo run in the small unstably stratified sublayer (close to ICB). Such weak dynamo action and skin-effect (the magnetic field generated in the unstably stratified sublayer permeates through the stably stratified 1 hereinafter referred to as CMB 2 hereinafter referred to as ICB

1

*Manuscript

A c c e p t e d M a n u s c r i p t sublayer where it is damped due to skin-effect) lead to the weak magnetic field observed on the surface of Mercury [START_REF] Christensen | A deep dynamo generating Mercury's magnetic field[END_REF].

For a different geometric configuration the influence of a non-uniform stratification is fundamental [START_REF] Stanley | Convective-region geometry as the cause of Uranus' and Neptune's unusual magnetic fields[END_REF][START_REF] Stanley | Numerical dynamo models of Uranus' and Neptune's magnetic fields[END_REF]. They assumed reverse stratification, i.e. the stably stratified sublayer is surrounded by the unstably stratified one. This configuration leads to non-dipolar and non-axisymmetric magnetic fields, which are typical e.g., for Uranus and Neptune.

The case, in which the thickness of both sublayers is the same, was not investigated in previous studies (e.g., [START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Šimkanin | Problem of the rotating magnetoconvection in variously stratified fluid layer revisited[END_REF][START_REF] Šimkanin | On hydromagnetic instabilities and the mean electromotive force in a non-uniformly stratified Earth's core affected by viscosity[END_REF][START_REF] Šimkanin | Control volume method for hydromagnetic dynamos in non-uniformly stratified spherical shells[END_REF]). Consequently, we located the change of the sign to the middle of the convective shell. A convection in rotating non-uniformly stratified spherical fluid shells was studied in dependence on Prandtl, Ekman and Rayleigh numbers. The model and governing equations are given in Section 2. The numerical results are presented in Section 3. Finally, Section 4 provides the conclusions.

Governing equations and model

Convection (with the velocity V) of incompressible fluid in the Boussinesq approximation in a non-uniformly stratified spherical shell (r i < r < r o ) rotating with angular velocity Ω is described by the system of dimensionless equations:

E ∂V ∂t + (V • ∇) V -∇ 2 V + 21 z × V + ∇P = R a r r o T, (1)
∂T ∂t + (V • ∇) T = 1 P r ∇ 2 T + H, (2)
∇ • V = 0. ( 3
)
The radius of the outer sphere L, is the typical length scale, which makes the dimensionless radius r o = 1; the inner core radius r i is, similar to that of the Earth, equal to 0.35. (r, θ, ϕ) is the spherical system of coordinates, 1 z is the unit vector. The typical time, t, is measured in the unit of L 2 /ν, typical velocity, V, in ν/L, temperature, T, in ∆T, and pressure, P , ρν 2 /L 2 . The dimensionless parameters appearing in Eqs (1)-( 3) are the Prandtl number, P r = ν/κ, the Ekman number, E = ν/ΩL 2 and the modified Rayleigh number R a = αg 0 ∆T L/Ων, where κ is thermal diffusivity, ν is kinematic viscosity, ρ is density, α is the coefficient of thermal expansion, ∆T is the drop of temperature through the shell and g 0 is the gravity acceleration at r = r o . Eqs (1)-( 3) are closed by the non-penetrating and no-slip boundary conditions for the velocity field at the rigid surfaces and fixed temperature boundary conditions.

The last term in Eq. ( 2), H, constitutes the internal heat sources (it is the non-dimensional volumetric heat source strength), which enable simulating thermodynamically the various stratification of the spherical shells also in the Boussinesq models. The outer sphere was assumed to be stratified nonuniformly (it is divided into stably and unstably stratified sub-shells) with constant temperature T i = 1 and T o = 0 at the inner and outer boundaries of the Page 4 of 18 A c c e p t e d M a n u s c r i p t shell, respectively. The width of both sub-shells was considered the same, i.e. ∂T /∂r changes its sign in the middle of the convective shell,

r m = (r i + r o )/2. It is accomplished if H = -23P -1
r . Let us describe why we used this value. The non-uniform stratification is considered due to internal heat sources and in the basic state Eq. (2) gets the form

1 P r ∇ 2 T + H = 0. ( 4
)
If we solve Eq. ( 4) considering T i = 1, T o = 0 and ∂T /∂r = 0 at r = r m , we obtain H = -23P -1 r .

Numerical results

Eqs (1)-(3) were solved using the MAG dynamo code [START_REF] Olson | Highly supercritical thermal convection in a rotating spherical shell: centrifugal vs. radial gravity[END_REF][START_REF] Olson | Magnetoconvection in a rotating spherical shell: structure of flow in the outer core[END_REF][START_REF] Olson | Magnetoconvection and thermal coupling of the Earth's core and mantle[END_REF][START_REF] Christensen | Numerical modelling of the geodynamo: a systematic parameter study[END_REF][START_REF] Christensen | A numerical dynamo benchmark[END_REF][START_REF] Christensen | Scaling properties of convection drivne dynamos in rotating spherical shells and application to planetary magnetic fields[END_REF]. It is a serial version of Gary Glatzmaier's rotating spherical convection/magnetoconvection/dynamo code, modified by Ulrich Christensen and Peter Olson. The code solves the non-dimensional Boussinesq equations for time-dependent thermal convection in a rotating spherical shell filled with an electrically conducting fluid. MAG uses toroidal-poloidal decomposition for velocity and magnetic field with explicit time steps. Linear terms are evaluated spectrally (spherical harmonics plus Chebyshev polynomials in radius) and nonlinear terms are evaluated on a spherical grid (for more details, see [START_REF] Olson | Highly supercritical thermal convection in a rotating spherical shell: centrifugal vs. radial gravity[END_REF][START_REF] Olson | Magnetoconvection in a rotating spherical shell: structure of flow in the outer core[END_REF][START_REF] Olson | Magnetoconvection and thermal coupling of the Earth's core and mantle[END_REF][START_REF] Christensen | Numerical modelling of the geodynamo: a systematic parameter study[END_REF][START_REF] Christensen | A numerical dynamo benchmark[END_REF][START_REF] Christensen | Scaling properties of convection drivne dynamos in rotating spherical shells and application to planetary magnetic fields[END_REF]). The computations were performed on the Sun Grid Engine (LUNA) at the Institute of Physics, Academy of Sciences of CR, Prague; the Nemo cluster (SGI) and PC clusters at the Institute of Geophysics, Academy of Sciences of CR, Prague. Dependence of solutions on various values of the Prandtl number, P r , the Ekman number, E, and the Rayleigh number, R a , was investigated. Computations started from zero initial velocity and were performed for P r = 10 -1 , 1, 10; E = 10 -3 , 10 -5 , and R a = 250, 600 (for E = 10 -3 ) and R a = 5000, 15000 (for E = 10 -5 ). The Prandtl number, the ratio of the thermal diffusion time to the viscous one, is the most significant dynamic factor in our study. According to the value of P r , it is possible to split the solutions to three qualitatively different groups, i.e. convection at low Prandtl numbers, P r = 1 and large Prandtl numbers.

Prandtl number P r = 1

Let us start with the most used value of P r . In the case of P r = 1 the characteristic thermal diffusion time is equal to the characteristic viscous diffusion time (τ κ = τ ν ), i.e. thermal and viscous diffusion processes equally affect the dynamics of convection. Dependence of the mean kinetic energy, E k , on R a , E and the type of stratification (uniform and non-uniform) for P r = 1 is given in Fig. 1. In the case of uniform stratification the whole shell is unstably stratified.

E k increases for given E with increasing of R a and also increases with a decreasing of E (the last one is given by the increasing of the critical Rayleigh Page 5 of 18 number, R c , with decreasing of E). The dependence of E k on the type of stratification is weak, almost negligible (see Fig. 1). The typical space distributions of the radial velocity field component, V r are presented in Fig. 2, Z-vorticity in Fig. 3 (equatorial sections) and the azimuthal velocity field component, V ϕ in Fig. 4 (axi-symmetrical meridional sections). All the figures are snapshots done at time t = 3.5 (3.5 time units). Fig. 2 (and the corresponding Figs. ( 6) and ( 10)) were done using the value of r = 0.69 which lies in the stable region but close to the interface r = r m (r m = 0.675). Consequently, we used the value of r practically equal to r m . Dependence on the type of stratification is now more evident than of E k (see Figs. [START_REF] Busse | Convection in rotating spherical fluid shells and its dynamo states[END_REF][START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF][START_REF] Christensen | Numerical modelling of the geodynamo: a systematic parameter study[END_REF]). It is more visible for small values of R a (large values of E), for large values of R a (small values of E) is almost negligible. For small values of R a in the case of uniform stratification, a large-scale columnar convection is developed in the whole volume, while in the case of non-uniform stratification it is suppressed to the unstably stratified region but slightly penetrates to the Page 6 of 18 stably stratified one (see Figs. [START_REF] Busse | Convection in rotating spherical fluid shells and its dynamo states[END_REF][START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF]). This behaviour is known from previous analyses [START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Šimkanin | Problem of the rotating magnetoconvection in variously stratified fluid layer revisited[END_REF][START_REF] Šimkanin | On hydromagnetic instabilities and the mean electromotive force in a non-uniformly stratified Earth's core affected by viscosity[END_REF][START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF]. Increasing R a it is possible to observe a transition from the large-scale convection to the small-scale convection (Figs. [START_REF] Busse | Convection in rotating spherical fluid shells and its dynamo states[END_REF][START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF][START_REF] Christensen | Numerical modelling of the geodynamo: a systematic parameter study[END_REF] from left to right). The vorticity is the much larger in the cases of high Rayleigh numbers (see Fig. 3). For E = 10 -5 and R a = 5000, 15000 only a small-scale convection is developed.

The spiralling cross section of the columns as a dominant feature at moderate Prandtl numbers, causes the Reynolds stress whose action generates a differential rotation, i.e. the spiralling cross section of the columns causes a differential rotation [START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF][START_REF] Zhang | Spiralling columnar convection in rapidly rotating spherical fluid shells[END_REF]. This feature, observable for small values of R a , is in our case for large values of R a unseen. Some irregularities in the spiralling structure of columnar convection could be probably caused by the differential rotation, so-called "active" and "quiet" zones [START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF]. Let us briefly comment our lack of columnar structures at lower values of E, when for similar values of the Rayleigh numbers, a convection at E = 10 -5 would be much more columnar than convection at E = 10 -3 . Our values of R a are very large in the case of E = 10 -5 , so large that the effect of rotation is scarcely visible (buoyancy forces are much stronger than Coriolis forces). This is apparent in Fig. meridional cross-section of the radial velocity (not axi-symmetric just any section) would make it clearer whether the convection has a columnar structure or not.

The multilayer convection ("teleconvection") is not developed in our case of non-uniform stratification because of the significant amount of stable stratification. An excitation of toroidal motions in the outermost part of the shell ("teleconvection") is a consequence of the radial stratification, without which they cannot take place, and is possible if the unstable stratification dominates over the stable one [START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF].

Low Prandtl numbers

In the case of low Prandtl numbers (P r < 1) the characteristic thermal diffusion time is smaller than the characteristic viscous diffusion time (τ κ < τ ν ), i.e. thermal diffusion processes dominate over viscous ones. Inertial modes become important (see left side of Eq. ( 1), where the inertial term is not negligible in comparison with the viscous one). The main difference between convection at P = 10 -1 (we put P = 10 -1 in the case of low Prandtl numbers) and P = 1 occurs at the transition between regular and irregular patterns [START_REF] Busse | Convection in rotating spherical fluid shells and its dynamo states[END_REF]. Dependence of the mean kinetic energy, E k , on R a , E and the type of stratification (uniform Page 8 of 18 ), R a = 5000, 15000 (both for E = 10 -5 ) and P r = 1 (from left to right). Red (blue) colours indicate positive (negative) values. Minimal and maximal values for the first line are (-18, 18), (-47, 47), (-89, 89), (-264, 264), and for the second line (-27, 27), (-61, 61), (-96, 96), (-304, 304). Snapshots at t = 3.5. and non-uniform) for P r = 10 -1 is given in Fig. 5.

E k again increases for given E with increasing of R a and also increases with a decreasing of E. The dependence of E k on the type of stratification is again weak, almost negligible, but for E = 10 -5 and R a = 15000 is strong (see Fig. 5).

The typical space distributions of the radial velocity field component, V r are presented in Fig. 6, Z-vorticity in Fig. 7 (equatorial sections) and the azimuthal velocity field component, V ϕ in Fig. 8 (axi-symmetrical meridional sections). All the figures are snapshots done at time t = 3.5. Dependence on the type of stratification is now again more evident than of E k (see Figs. [START_REF] Christensen | A deep dynamo generating Mercury's magnetic field[END_REF][START_REF] Christensen | Scaling properties of convection drivne dynamos in rotating spherical shells and application to planetary magnetic fields[END_REF][START_REF] Fearn | Compositional convection and stratification of the Earth's core[END_REF]). For R a = 250 (left columns) and R a = 15000 (right columns) it is strong (except for the radial component of velocity, V r , displayed in Fig. 6 where it is weak). In the case of uniform stratification for R a = 250 a large-scale columnar convection with the strong spiralling nature is developed in the whole Page 9 of 18 volume, in the case of non-uniform stratification it is again suppressed to the unstably stratified region (shifted out from CMB) but slightly penetrates to the stably stratified one (see Figs. [START_REF] Christensen | A deep dynamo generating Mercury's magnetic field[END_REF][START_REF] Christensen | Scaling properties of convection drivne dynamos in rotating spherical shells and application to planetary magnetic fields[END_REF]). The spiralling character of a columnar convection causes a differential rotation [START_REF] Busse | Convection in rotating spherical fluid shells and its dynamo states[END_REF][START_REF] Zhang | Spiralling columnar convection in rapidly rotating spherical fluid shells[END_REF]. For R a = 15000 it is possible to observe a small-scale convection without a strong spiralling because this value of R a is far more strongly supercritical, and buoyancy forces are much stronger than Coriolis forces.

Large Prandtl numbers

In the case of large Prandtl numbers (P r > 1) the characteristic thermal diffusion time is greater than the characteristic viscous diffusion time (τ κ > τ ν ), i.e. viscous diffusion processes dominate over thermal ones. The typical difference between convection at P r = 1 and P r = 10 (we put P = 10 in the case of large Prandtl numbers) is that the convection columns retain their alignment Page 10 of 18 with the axis of rotation with increasing P r , but the spiralling nature of their radial orientation disappears [START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF]. As P r increases the influence of the differential rotation, which dominates the evolution of the convection columns for P r of the order unity and bellow diminishes rapidly [START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF]. As in previous cases E k increases for given E with increasing of R a but its increasing with a decreasing of E is not so strong as for P r = 1 and P r < 1. The dependence of E k on the type of stratification is strong for E = 10 -3 but for E = 10 -5 it is weak, almost negligible (see Fig. 9).

The typical space distributions of the radial velocity field component, V r are presented in Fig. 10,Z-vorticity in Fig. 11 (equatorial sections) and the azimuthal velocity field component, V ϕ in Fig. 12 (axi-symmetrical meridional sections). All the figures are snapshots done at time t = 3.5. Dependence on the type of stratification is now again more evident than of E k (see Figs. [START_REF] Heimpel | Simulation of equatorial and high-latitude jets on Jupiter in a deep convection model[END_REF][START_REF] Jones | Convection-driven geodynamo models[END_REF][START_REF] Olson | Highly supercritical thermal convection in a rotating spherical shell: centrifugal vs. radial gravity[END_REF]). For R a = 250 (left columns) it is strong but for R a = 15000 (right columns) it is weak (except for the azimuthal component of velocity, V ϕ , displayed in Fig. 12 where it is strong). For R a = 250 a convection is developed in the whole volume in the both cases of uniform and non-uniform stratifications. It is interesting that in the case of non-uniform stratification the convection significantly penetrates to the stably stratified region (it is not suppressed to the unstably Page 11 of 18 stratified one as it is typical for P r = 1 and P r < 1, see Figs. [START_REF] Jones | Convection-driven geodynamo models[END_REF], ( 3), ( 7)).

For R a = 15000 it is possible to observe a small-scale-convection (similarly as in previous cases). The azimuthal flow is generally weaker at high values of P r than low ones (see Figs. ( 4), ( 8), ( 12)).

Conclusions

A density stratification plays a significant role in the astrophysical and geophysical processes. It is an important dynamical factor in the processes acting in planetary atmospheres (e.g., [START_REF] Heimpel | Simulation of equatorial and high-latitude jets on Jupiter in a deep convection model[END_REF]), planetary interiors (e.g., [START_REF] Braginsky | Magnetohydrodynamics of the Earth's Core[END_REF][START_REF] Stanley | Convective-region geometry as the cause of Uranus' and Neptune's unusual magnetic fields[END_REF][START_REF] Stanley | Numerical dynamo models of Uranus' and Neptune's magnetic fields[END_REF][START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Stanley | Effects of an outer thin stably stratified layer on planetary dynamos[END_REF][START_REF] Christensen | A deep dynamo generating Mercury's magnetic field[END_REF][START_REF] Šimkanin | Problem of the rotating magnetoconvection in variously stratified fluid layer revisited[END_REF][START_REF] Šimkanin | On hydromagnetic instabilities and the mean electromotive force in a non-uniformly stratified Earth's core affected by viscosity[END_REF]), etc.

A systematical parameter study of rotating convection in uniformly and non-uniformly stratified spherical shells in dependence on Prandtl, Ekman and Rayleigh numbers was done e.g., in [START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Busse | Convection in rotating spherical fluid shells and its dynamo states[END_REF][START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF]. Our attention is focused on the case when the thickness of both sublayers (stable and unstable) is the same. Such a case was not investigated in previous studies (e.g., [START_REF] Zhang | On coupling between the Poincare equation and the heatequation[END_REF][START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Šimkanin | Problem of the rotating magnetoconvection in variously stratified fluid layer revisited[END_REF][START_REF] Šimkanin | On hydromagnetic instabilities and the mean electromotive force in a non-uniformly stratified Earth's core affected by viscosity[END_REF][START_REF] Šimkanin | Control volume method for hydromagnetic dynamos in non-uniformly stratified spherical shells[END_REF]). It is true that neither planet nor satellite with this geometric configuration has been observed. Nevertheless, it is interesting to investigate such a case in order to answer a question whether its behaviour is closer to the case of a stable stratification or an unstable one. Thus, we investigate our model in dependence Page 12 of 18 on Prandtl, Ekman and Rayleigh numbers, and the type of stratification (in order to compare an influence of a non-uniform stratification, we investigated also the case of an unstably stratified shell).

The Prandtl number is the most significant dynamic factor in our study. Consequently, we split the solutions to three qualitatively different groups, i.e. convection at low, moderate and large Prandtl numbers. The mean kinetic energy, E k , increases for given E with increasing of R a and also increases with a decreasing of E (the last one is given by the increasing of the critical Rayleigh number, R c , with decreasing of E). The dependence of E k on the type of stratification is weak, almost negligible except the cases of E = 10 -5 , R a = 15000 for P r = 10 -1 and E = 10 -3 for P r = 10 when it is strong (see Figs.

(1), ( 5), ( 9)). Dependence on the type of stratification is more evident on the typical space distributions of the velocity field and vorticity, mainly for small values of R a (large values of E). For large values of R a (small values of E) is almost negligible. For small values of R a in the case of uniform stratification a Page 13 of 18 large-scale columnar convection is developed in the whole volume, in the case of non-uniform stratification it is suppressed to the unstably stratified region but slightly penetrates to the stably stratified one, except the case of P r = 10 when the convection significantly penetrates to the stably stratified region (it is not suppressed to the unstably stratified one). The vorticity is the much larger in the cases of high Rayleigh numbers for all P r (see Figs.

(3), ( 7), ( 11)). The azimuthal flow is generally weaker at high values of P r than low ones (see Figs. ( 4), ( 8), ( 12)). The spiralling cross section of the columns which causes a differential rotation, is a dominant feature at P r = 1 and becomes stronger at P r = 10 -1 , but it disappears at P r = 10. For large values of R a (small values of E) only a small-scale convection is developed (it is possible to observe smallscale structures). However, it is necessary to remark again that our values of R a for E = 10 -5 are far more strongly supercritical, and buoyancy forces are much stronger than Coriolis forces (the effect of rotation is scarcely visible). This is apparent in Figs. structure, but instead resemble non-rotating convection, e.g., no spiralling even at low P r (which is characterized by strong spiralling) in the rightmost figures of Figs.

(3) and [START_REF] Christensen | Scaling properties of convection drivne dynamos in rotating spherical shells and application to planetary magnetic fields[END_REF]. Consequently, the plots of a meridional cross-section of the radial velocity would make it clearer whether the convection has a columnar structure or not. The multilayer convection ("teleconvection") is not developed in our case of non-uniform stratification because of the significant amount of stable stratification.

Our results are in agreement with previous analyses done in the uniformly (e.g., [START_REF] Busse | Convection in rotating spherical fluid shells and its dynamo states[END_REF]) or non-uniformly stratified spherical shells (e.g., [START_REF] Zhang | On coupling between the Poincare equation and the heatequation[END_REF][START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Šimkanin | Problem of the rotating magnetoconvection in variously stratified fluid layer revisited[END_REF][START_REF] Šimkanin | On hydromagnetic instabilities and the mean electromotive force in a non-uniformly stratified Earth's core affected by viscosity[END_REF][START_REF] Šimkanin | Control volume method for hydromagnetic dynamos in non-uniformly stratified spherical shells[END_REF][START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF]). The convective motions were qualitatively described in [START_REF] Zhang | On coupling between the Poincare equation and the heatequation[END_REF][START_REF] Zhang | Teleconvection: remotely driven thermal convection in rotating stratified spherical layers[END_REF][START_REF] Busse | 3.4. Convection in rotating spherical fluid shells[END_REF], etc. We showed that the case characterized by the identical thicknesses of both sublayers is a typical extreme case. It is dissimilar to the case if the stably stratified sublayer is thinner than the unstably stratified one (although the convection penetrates to the stably stratified sublayer, it is more suppressed to the unstably stratified sublayer than in the case if r m is closer CMB) and also dissimilar to the case when the unstably stratified sublayer is thinner than the stably stratified one (although the convection is significantly suppressed to the unstably stratified sublayer, still penetrates to the stably stratified sublayer more than in the case Page 15 of 18

Figure 1 :

 1 Figure 1: Dependence of the mean kinetic energy, E k , on the Rayleigh number, R a , the Ekman number, E, and the type of stratification (uniform -UNI and non-uniform -NON-UNI) for P r = 1.

Figure 2 :

 2 Figure 2: Space distributions of radial velocity field components V r at r = 0.69 for the cases of uniform (top) and non-uniform (bottom) stratification, and for R a = 250, 600 (both for E = 10 -3 ), R a = 5000, 15000 (both for E = 10 -5 ) and P r = 1 (from left to right). Red (blue) colours indicate positive (negative) values. Minimal and maximal values for the first line are (-48, 48), (-136, 136), (-261, 261), (-791, 791), and for the second line (-78, 78), (-173, 173), (-296, 296), (-905, 905). Snapshots at t = 3.5.

Figure 3 :

 3 Figure 3: Equatorial sections of Z-vorticity for the cases of uniform (top) and non-uniform (bottom) stratification, and for R a = 250, 600 (both for E = 10 -3 ), R a = 5000, 15000 (both for E = 10 -5 ) and P r = 1 (from left to right). Red (blue) colours indicate positive (negative) values. Minimal and maximal values for the first line are (-108, 108), (-115, 115), (-349, 349), (-864, 864), and for the second line (-113, 113), (-123, 123), (-401, 401), (-957, 957). Snapshots at t = 3.5.

A

  c c e p t e d M a n u s c r i p t AZIMUTHAL VELOCITY AZIMUTHAL VELOCITY AZIMUTHAL VELOCITY AZIMUTHAL VELOCITY AZIMUTHAL VELOCITY AZIMUTHAL VELOCITY AZIMUTHAL VELOCITY AZIMUTHAL VELOCITY

Figure 4 :

 4 Figure4: Axi-symmetrical meridional sections of azimuthal velocity field components V ϕ for the cases of uniform (top) and non-uniform (bottom) stratification, and for R a = 250, 600 (both for E = 10 -3 ), R a = 5000, 15000 (both for E = 10 -5 ) and P r = 1 (from left to right). Red (blue) colours indicate positive (negative) values. Minimal and maximal values for the first line are (-18, 18), (-47, 47), (-89, 89),(-264, 264), and for the second line (-27, 27), (-61, 61), (-96, 96),(-304, 304). Snapshots at t = 3.5.

Figure 5 :

 5 Figure 5: Dependence of the mean kinetic energy, E k , on the Rayleigh number, R a , the Ekman number, E, and the type of stratification (uniform -UNI and non-uniform -NON-UNI) for P r = 10 -1 .

A c c eFigure 6 :

 6 Figure 6: Space distributions of radial velocity field components V r at r = 0.69 for the cases of uniform (top) and non-uniform (bottom) stratification, and for R a = 250 (for E = 10 -3 ), R a = 15000 (for E = 10 -5 ) and P r = 10 -1 (from left to right). Red (blue) colours indicate positive (negative) values. Minimal and maximal values for the first line are (-807, 807), (-2621, 2621), and for the second line (-779, 779), (-1440, 1440). Snapshots at t = 3.5.

Figure 7 :

 7 Figure 7: Equatorial sections of Z-vorticity for the cases of uniform (top) and non-uniform (bottom) stratification, and for R a = 250 (for E = 10 -3 ), R a = 15000 (for E = 10 -5 ) and P r = 10 -1 (from left to right). Red (blue) colours indicate positive (negative) values. Minimal and maximal values for the first line are (-1203, 1203), (-4095, 4095), and for the second line (-1128, 1128), (-2114, 2114). Snapshots at t = 3.5.

Figure 8 :

 8 Figure 8: Axi-symmetrical meridional sections of azimuthal velocity field components V ϕ for the cases of uniform (top) and non-uniform (bottom) stratification, and for R a = 250 (for E = 10 -3 ), R a = 15000 (for E = 10 -5 ) and P r = 10 -1 (from left to right). Red (blue) colours indicate positive (negative) values. Minimal and maximal values for the first line are (-271, 271), (-876, 876), and for the second line (-262, 262), (-485, 485). Snapshots at t = 3.5.

Figure 9 :

 9 Figure 9: Dependence of the mean kinetic energy, E k , on the Rayleigh number, R a , the Ekman number, E, and the type of stratification (uniform -UNI and non-uniform -NON-UNI) for P r = 10.

( 2 )Figure 10 :

 210 Figure 10: Space distributions of radial velocity field components V r at r = 0.69 for the cases of uniform (top) and non-uniform (bottom) stratification, and for R a = 250 (for E = 10 -3 ), R a = 15000 (for E = 10 -5 ) and P r = 10 (from left to right). Red (blue) colours indicate positive (negative) values. Minimal and maximal values for the first line are (-5.5, 5.5),(-38, 38), and for the second line (-14, 14),(-42, 42). Snapshots at t = 3.5.

Figure 11 :Figure 12 : 5 .

 11125 Figure 11: Equatorial sections of Z-vorticity for the cases of uniform (top) and non-uniform (bottom) stratification, and for R a = 250 (for E = 10 -3 ), R a = 15000 (for E = 10 -5 ) and P r = 10 (from left to right). Red (blue) colours indicate positive (negative) values. Minimal and maximal values for the first line are (-8.3, 8.3), (-54, 54), and for the second line (-21, 21), (-65, 65). Snapshots at t = 3.5.

Acknowledgements

This research has been supported by the Grant Agency of the Academy of Sciences of the Czech Republic (Grant No. IAA300120704). We would like to thank G. Glatzmaier, U. Christensen, P. Olson, and CIG for MAG dynamo code, and the Institute of Physics of Academy of Sciences of CR, Prague for CPU time on the Sun Grid Engine (the Luna project). Our thanks includes also two unknown referees, whose comments greatly improved this paper.


Alexander S Côté

Lidunka Vočadlo

David Dobson

Dario Alfè

John P Brodholt

Ab initio lattice dynamics calculations on the combined effect of temperature and silicon on the stability of different iron phases in the Earth's inner core

come

Ab initio lattice dynamics calculations on the combined

effect of temperature and silicon on the stability of different iron phases in the Earth's inner core

Constraining the crystalline phase in the inner core has a long history (e.g. Anderson 1985[START_REF] Boehler | The phase diagram of iron to 430 kbar[END_REF][START_REF] Saxena | Experimental evidence for a new iron phase and implications for the Earth's core[END_REF], 2000[START_REF] Stixrude | Constraints on the crystalline structure of the inner core: mechanical instability of BCC iron at high pressure[END_REF][START_REF] Andrault | The orthorhombic structure of iron: An in situ study at high temperature and high pressure[END_REF], Alfè et al. 2000a,b, Lin et al. 2002[START_REF] Vočadlo | First principles calculations on crystalline and liquid iron at Earth's core conditions[END_REF][START_REF] Vočadlo | The structure of iron under the conditions of the Earth's inner core[END_REF], 2003[START_REF] Belonoshko | Stability of the body-centredcubic phase of iron in the Earth's inner core[END_REF], 2008[START_REF] Dubrovinsky | Body-centred cubic iron-nickel alloy in Earth's core[END_REF], Côté et al. 2008a,b). Under ambient conditions, the iron crystal has a body-centred-cubic (bcc) structure. Increasing temperature at ambient pressure causes a transformation to the face-centred-cubic (fcc) phase above ~1150 K. At pressures above ~11 GPa, (at room temperature), it adopts a hexagonal-close-packed (hcp) structure which persists to Earth's inner core pressures. Pure Fe in the bcc phase was previously excluded as a possible crystal structure in the core, because it is mechanically unstable at high pressures [START_REF] Stixrude | Iron at high pressure: Linearizedaugmented-plane-wave computations in the generalized-gradient approximation[END_REF][START_REF] Söderlind | First-principles theory of iron up to earth-core pressures: Structural, vibrational, and elastic properties[END_REF]. However, the effect of temperature has proven to be very important, and the stable phase at core pressures and temperatures remains ambiguous. High temperature has been shown to decrease the free energy difference between bcc and hcp phases (Vočadlo et al., 2003), and molecular dynamics calculations have found bcc to be mechanically stable at core conditions [START_REF] Belonoshko | Stability of the body-centredcubic phase of iron in the Earth's inner core[END_REF], Vočadlo et al., 2003). Recently the face-centred cubic (fcc) phase was proposed as a third possible structure in the inner core, with a free energy between those of bcc and hcp [START_REF] Vočadlo | Ab initio calculations on the free energy and high P-T elasticity of face-centred-cubic iron[END_REF]. Since the core is likely to be an iron alloy, the effect of alloying elements on the crystal structure has to be considered. In a previous ab initio study we showed that at core pressure and zero temperature, light elements (Si, S, C, O) tend to stabilise the bcc phase with respect to hcp, but hcp still remains the most stable phase (Côté, et al. 2008a,b). Here we present the results of ab initio lattice dynamics (LD) calculations, to evaluate the effect of temperature on the crystal structure of iron and Fe-Si alloys at core pressures.

A c c e p t e d M a n u s c r i p t 3

Computational Methods

The Density Functional Theory (DFT) calculations in this work were carried out using the VASP code [START_REF] Kresse | Efficient iterative schemes for ab initio totalenergy calculations using a plane-wave basis set[END_REF], which incorporates the projected augmented wave (PAW) method [START_REF] Blöchl | Projector augmented-wave method[END_REF]. The exchange-correlation energy was represented by the generalised gradient approximation (GGA) [START_REF] Wang | Correlation hole of the spin-polarized electron gas, with exact small-wave-vector and high-density scaling[END_REF]. VASP calculates the ground state (T = 0) for each set of ionic positions and the electronic free energy is taken as the quantity to be minimised. Care was taken when choosing the k-point sampling grid and cutoff energies in order to have an energy convergence of no more than 1 meV/atom. Specifically, we used a 400eV cutoff energy throughout the calculations. For pure iron bcc and fcc in the 2atom cells, 12x12x12 (equivalent to 56 k-points in the irreducible Brillouin zone (IBZ)) and 12x12x8 (84 k-points in the IBZ) k-point grids were used respectively.

For hcp a 4-atom cell was used, with a 15x9x9 grid (200 k-points in the IBZ). Larger primitive cells (8-16 atoms) were used for the alloyed structures, and correspondingly smaller k-point sampling grids.

We calculated the phonon frequencies for different structures and examined their vibrational stability. The phonon calculations were carried out using the code PHON [START_REF] Alfè | PHON: A program to calculate phonons using the small displacement method[END_REF], which calculates the phonon vibrational frequencies ω ks along a predefined path in the 1 st Brillouin zone, by using the small displacement method to construct force constant matrices. The force-constant matrix Φ isα, jtβ is the central quantity in the calculation of the phonon frequencies, since the (squares of the) frequencies at wavevector k are the eigenvalues of the dynamical matrix D sα,tβ , defined as:

A c c e p t e d M a n u s c r i p t

              i s 0 i t 0 j jt , is t s t , s τ τ i M M 1 ) ( R R k Φ k D exp (1)
where R i 0 is a vector of the lattice connecting different primitive cells and τ s is the position of the atom s in the primitive cell. If we have the complete force-constant matrix, then Dsα, tβ and hence the frequencies ω ks can be obtained at any . In principle, the elements of Φ isα, jtβ are non-zero for arbitrarily large separations

| R j 0 + τ t -R i 0 -τ s |
, but in practice they decay rapidly with separation, so PHON neglects all the elements beyond a certain cut-off distance resulting in a force constant matrix of finite extent.

Our calculations were extended to nonzero temperatures by using the quasi-harmonic approximation, where the volume and temperature dependence of the Helmholtz free energy F(V,T) was calculated as a sum of the contributions due to static compression, F 0 (V), thermal excitation of electrons, F el (V,T) and thermal excitation of phonons, F vib (V,T):

F(V,T) = F 0 (V) + F el (V,T) + F vib (V,T) (2)
The last term in (2) is defined as:

                 k k , 2 sinh 2 ln ) , ( s B s B vib T k h T k T V F (3)
where k B is the Boltzman constant, and ω sk is the frequency of the phonon mode for wave vector k and volume V. The contribution to the free energy from electronic excitations at different electronic temperatures F el (V,T) was included by using a Fermi smearing function when calculating the ground-state energy in VASP.

The Gibbs free energy was then obtained by

G(P,T) = F(V,T)+ PV (4).
In order to account for the contribution of thermal pressure at different temperatures and obtain the last term of (4) accurately, we performed free energy calculations at 6 different volumes, corresponding to pressures in the range 200-400 GPa. We then fitted the results to an E(V) equation of state, from where we were able to find the correct volume for any pressure at various temperatures. The total free energies for core pressures and temperatures could thus be evaluated.

It should be noted that the calculations are quasi-harmonic and, therefore, neglect anharmonicity. The melting temperature of Fe under core pressures is about 6000 K so anharmonicity may start to become important after about 4000 K or so. However, the total anharmonic contribution to the free-energy of hcp Fe at 6000 K is only of the order of 60 meV/atom [START_REF] Alfè | Thermodynamics of hexagonal close packed iron under Earth's core conditions[END_REF]. Since fcc is a close-packed structure like hcp, it should be affected by anharmonicity in a similar way. It is thus unlikely for the free-energy differences to be larger than 10 meV/atom or so, even at very high temperatures. BCC, however, is not close packed, and so the anharmonic contribution to the free-energy may be higher.

Results.

A c c e p t e d M a n u s c r i p t

i. Vibrational stability

To get accurate results in the harmonic approximation, very accurate forces are needed. We therefore used a very low tolerance criterion for the total energy (10 -7 eV), and an added grid in the Fourier transform mesh which serves to reduce the noise in the forces, which were set to converge below 10 -5 eV/Å. The atoms were displaced a small distance (0.01Å) from their equilibrium positions, in order for the forces to remain within the harmonic approximation, and the forces were calculated. For all cases, the calculations were performed on supercells large enough to avoid selfinteraction due to the periodic boundary conditions, typically 3 times the size of the primitive unit cell in each direction. The resulting supercells contained up to 288 atoms.

Firstly we calculated phonon frequencies for pure Fe at different pressures and electronic temperatures. The calculated dispersion relations for bcc Fe for certain special directions of the 1 st Brillouin zone are displayed in figure 1. [insert figure 1 here]

As reported previously (Vočadlo et al., 2003), at low temperatures the bcc phase (figure 1a) becomes vibrationally unstable at high pressures above ~200 GPa.

Furthermore, as can be seen from the figure (1b), at pressures below 250 GPa, high electronic temperatures contribute to destabilising the structure. This is manifested in a soft mode in the [110] (M-Γ) direction, which starts from zero pressure, and is due to bcc transforming to the fcc phase. Approaching core pressures, high temperatures do act towards eliminating some imaginary frequencies in certain directions and stabilising the structure, but overall the bcc phase remains dynamically unstable. High pressures in combination with electronic excitations at high temperatures also act to A c c e p t e d M a n u s c r i p t 7 suppress magnetic moments, and this was confirmed by spin polarised calculations in this work, as well as in previous studies (Vočadlo et al., 2003).

The hcp and fcc phases of pure iron are both dynamically stable at core pressures, and the phonon dispersion relations at T el = 5500K and 330 GPa are presented in figure 2.

[insert figure 2 here]

We then substitutionally inserted small concentrations of silicon at high pressure in all three phases. As can be seen in figure 3, the imaginary frequencies in the bcc phase disappear with increasing Si concentration, making the Fe-Si alloy vibrationally stable at 10.4 wt. % (12.5 at.%) and above. This is consistent with the experimental results

for FeSi, where a bcc (B2) structure is stable at core conditions [START_REF] Dobson | A new high-pressure phase of FeSi[END_REF] [insert figure 3 here]. In contrast, we find that increasing Si concentration de-stabilises the hcp structure; at 25 at.% Si, hcp is completely unstable. It should be noted that in the alloyed structure, magnetic moments disappear with high pressure even at low temperatures.

ii. Thermodynamical stability

To compare the thermodynamical stability of bcc, hcp and fcc we used a 10.4 wt.% We then focused on hcp and fcc -thermodynamically the most stable structures -at three different Si concentrations, 3.2 wt.%, 6.7 wt% and 10.4 wt.%. We repeated the phonon calculation for all concentrations and derived the free energies. We can see from figure 6 that as the Si concentration increases, fcc becomes the most stable phase at concentrations greater than ~7 wt.% (~13 at.%), at inner core pressures and temperatures. [insert figure 6 here] Figure 6 represents the pseudounivariant at 5500K, where one phase directly transforms into the other, but in reality the two phases would mix to form a solid solution. Thus, the free energy of the solid solution was derived by adding the configurational term k B T(c si ln(c si ) + c Fe ln(c Fe )) to the energies of the individual phases, where c represents the atomic % concentration. The solid solution energies were then fitted to a polynomial for different concentrations, and by taking the common tangent of the two curves, the coexisting compositions of the fcc and hcp phases were calculated. Using these, we derived an approximate phase diagram for the Fe-Si system, shown in figure 7, for 330 GPa and 360 GPa. Also plotted in figure 7 are the Fe-Si solid solution melting temperatures and compositions, predicted by Alfè et al. (2002). Note that the narrow width of our predicted two phase region is completely consistent with the narrow solid solution melting curve inferred from Alfè al (2008), and data from that study are added in figure 7 with open circles. It should be pointed out, however, that subsequent experimental studies were not able to observe fcc [START_REF] Lin | Phase relations of Fe-Si alloy in the Earth's core[END_REF][START_REF] Kuwayama | Phase relations of iron-silicon alloys at high pressure and high temperature[END_REF]). [insert figure 7 here] 4. Discussion.

Earlier ab initio finite temperature molecular dynamics calculations using the technique of thermodynamic integration to calculate free energies (Vočadlo et al., 2003) showed that the thermodynamically stable phase of pure iron at core conditions is the hcp phase. This study confirms that result. Earlier ab initio studies also showed that light elements act to decrease the free-energy difference bcc and hcp, but these were only performed at 0 K, and hcp still remained the stable phase (Côté et al. 2008a,b). It was, however, hypothesised that the combination of high temperature and light elements may be sufficient to make either bcc or fcc stable over the hcp phase.

This is what we have tested here. We find that at a concentration of ~7 wt.% silicon, the combination of high temperatures (over 4000 K) and light element enrichment act to make the fcc phase thermodynamically stable at core conditions. Hcp remains the stable phase at lower pressures and temperatures, as seen in a recent experimental study by [START_REF] Lin | Phase relations of Fe-Si alloy in the Earth's core[END_REF] where they used a concentration of ~8wt.% silicon and reached pressures of 240 GPa and temperatures of 3000 K. We also find that the bcc structure becomes vibrationally stable as more light elements are added, so depending on the concentration of the alloying element and the effect of anhamonicity, it cannot be discarded as an inner core phase.

A c c e p t e d M a n u s c r i p t 10

The light element concentration is generally thought to be about 4 wt.% (7.7 at.%) for the inner core and about 10 wt.% (18.1 at.%) for the outer core [START_REF] Lin | Ironsilicon alloy in the Earth's core?[END_REF][START_REF] Badro | Effect of light elements on the sound velocities in solid iron: Implications for the composition of Earth's core[END_REF][START_REF] Chen | Thermal expansion of iron-rich alloys and implications for the Earth's core[END_REF]. First of all, figure 7 shows that even if Si was the only light element in the inner core, a composition of ~8 at.% places the inner core firmly in the hcp-field. The fcc phase does not become stable until the concentration of Si reaches 11 at.%. Secondly, the higher concentration of light element in the outer core could mean that the core is crystallising fcc and hcp at the eutectic. However, this would predict an inner composition that would be too high. Taken at face value, figure 7 predicts that the inner core must be in the hcp phase and that fcc is not stable.

However, this phase diagram is for two components only. The other light elements are also known to stabilise bcc, so they may also stabilise fcc. If so, the width of the twophase region should widen, and it only needs to widen by ~3 at% for the inner core to enter the two-phase field. This is supported by the results of Lao and Cohen (presented at APS, 2009) who also showed that nickel (Ni) also strongly stabilises the fcc phase, and pushes the core into the two phase region. We suggest, therefore, that the inner core exists in the two-phase region, with hcp and fcc coexisting throughout the core. The phase diagram also shows that the proportion of fcc to hcp changes slightly with temperature; the deeper inner core may, therefore, contain a greater proportion of fcc that the shallower inner core.

Finally, the S-wave velocity of the inner core is known to be less than that of all possible candidate phases, including the effect of light elements; this motivated [START_REF] Vočadlo | Ab initio calculations of the elasticity of iron and iron alloys at inner core conditions: Evidence for a partially molten inner core?[END_REF] to suggest that the inner core contains a small fraction of partial melt.

However, the coexistence of fcc and hcp may allow the elastic velocities to be reduced via the mechanism proposed by [START_REF] Li | Effect of phase transitions on compressional-wave velocities in the Earth's mantle[END_REF] for the Earth's mantle. Li and Weidner showed that P-wave velocities of a material will be

  iron -silicon alloy (3 Si atoms and 13 Fe atoms), a concentration at which all three structures are vibrationally stable. [insert figure 4 here]. As there is more than one possible Si defect, a number of trials were performed with all the different configurations of the defect atoms in all phases, and the most favourable configurations at high pressure were chosen. They are shown in figure 4. The Gibbs free energy vs. temperature at a pressure of 330 GPa is plotted in figure 5, relative to the hcp phase [insert figure 5 here].

  that the free-energy difference between hcp and the other two phases diminishes significantly with increasing temperature. The free energy difference between hcp and bcc with 10.4 wt% Si reduces to only 36 meV/atom at 330 GPa and 5500 K, although the free energies are still in favour of hcp. More interestingly, fcc becomes the most stable structure at temperatures above 4000K.

  's results. The phase diagram is also consistent with experiments by Asanuma et.

  the pressure increase during the passage of the wave causes a small amount of the material to transform from one phase to another. When the two phases are coexisting, this can be a very efficient process, as they showed for the olivine to wadsleyite transition. Although Li and Weidner demonstrated this for compressional waves, shear waves should have an analogous, and potentially larger effect. This mechanism may well occur between fcc and hcp during the passage of a seismic wave in the inner core, effectively reducing the velocity from its elastic, unrelaxed values, and explaining the low observed S-wave velocities in the inner core.

Figure 1 :

 1 Figure captions

Figure 2 :

 2 Figure 2: Phonon dispersion of (a) hcp and (b) fcc phases of pure iron at Earth's inner

Figure 3 :

 3 Figure 3: Phonon dispersion relations for the bcc phase at 330 GPa for (a) pure iron,

Figure 4 :

 4 Figure 4: The most stable configurations used for (a) bcc, viewed along <101>, (b)

Figure 5 :

 5 Figure 5: Relative Gibbs free energies of the bcc and fcc structures with respect to

Figure 6 :

 6 Figure 6: Gibbs free energy of fcc and hcp at Earth's inner core conditions, relative to

Figure 7 :

 7 Figure 7: Compositional phase diagram for the Fe-Si system (at. %) at the inner core

  Figure 1

Corresponding author. Tel: +44 (0)20 7679

2425; Fax: +44 (0)20 7679 2685, Email address: a.cote@ucl.ac.uk

Acknowledgements

The work of ASC was supported by NERC grant NE/C519662/1. The computations were performed on the HECToR service.