text:
# Surprises from a deep ASCA spectrum of a Broad Absorption Line Quasar PHL 5200

## 1 Introduction

Broad absorption line quasars, in which the kinetic energy carried out in the absorbing outflow is a significant fraction of the bolometric luminosity of the quasar, offer a challenge to our understanding of the quasar energy budget (as suggested in Mathur, Elvis & Wilkes 1995; see also $`\mathrm{\S }11.3`$ in Krolik 1999). At the same time, they also offer new insights into the nuclear structure of quasars (Ogle 1998, Elvis 2001). X-ray observations are important in this investigation as they offer precise measurements of absorbing column densities. Combined UV and X-ray analysis is a powerful tool to understand the physical conditions in the absorbing gas (Mathur et al. 1998 and references there in). However, X-ray observations of BALQSOs have essentially resulted in non-detections (Green & Mathur 1996). While the available evidence clearly suggested that the observed X-ray weakness of BALQSOs is a result of strong X-ray absorption (Mathur et al. 2000 and references there in), in some objects absorption was not apparent (Gallagher et al. 1999). So the nature of X-ray weakness of BALQSOs and the column density of absorption if any, could not be determined with certainty for the lack of a good X-ray spectrum.

PHL 5200, the prototype BALQSO at z=1.98 (Burbidge 1968) was clearly detected in EXOSAT medium energy experiment (ME), but not by the low energy experiment (LE) (Singh, Westergaard & Schnopper 1987). To obtain consistency between ME and LE requires a column density of $`\genfrac{}{}{0pt}{}{_>}{^{}}`$1022 cm-2 at the source. The detection of PHL 5200 with EXOSAT made it an excellent candidate for detailed spectral study in X-rays. A 1994 ASCA observation of PHL 5200 yielded the first X-ray spectrum of a BALQSO (Mathur, Elvis & Singh 1995, MES95 here after), and remained the only one until this year. However, the quality of the ASCA spectrum was poor and the parameters of the fit were not well determined in the 17.7 ksec observation. The photon index of the best-fit power-law was uncertain to $`\pm 0.9`$ (all the errors are quoted to 90% confidence, unless noted otherwise). While excess absorption with NH =$`1.3_{1.1}^{+2.3}\times 10^{23}`$ cm-2 at the source provided a better fit, a model with only Galactic absorption was acceptable (see also Gallagher et al. 1999). A good signal to noise spectrum of a BALQSO was very much needed and PHL 5200 remained the best target. So we re-observed PHL 5200 with ASCA, to obtain better quality data and better constrain the parameters of an X-ray spectrum of a BALQSO. We find that the spectrum contains more surprises that affect the interpretation of the BAL phenomenon.

## 2 Observations and Data Analysis

### 2.1 Observations

We observed PHL 5200 with ASCA (Tanaka et al. 1994) on 1999 November 21. ASCA contains two sets of two detectors, SIS (Solid-state Imaging Spectrometer) and GIS (Gas Imaging Spectrometer). The effective exposure times in SIS0, SIS1, GIS2 and GIS3 were 83,522 seconds, 73,320 seconds, 93,385 seconds and 93,337 seconds respectively. The SIS was operated in 1CCD mode with the target in the standard 1CCD mode position. The GIS was operated in pulse height (PH) mode. The data were reduced and analyzed using FTOOLS and XSELECT in a standard manner (see ASCA Data Reduction Guide for details of data reduction).

### 2.2 Data Analysis

The source was clearly detected in all the four instruments with background subtracted net count rate of 8.6$`\pm 0.5\times 10^3`$, 4.8$`\pm 0.4\times 10^3`$, 5.6$`\pm 0.6\times 10^3`$, 10.9$`\pm 0.6\times 10^3`$ for SIS0, SIS1, GIS2 and GIS3 respectively. We extracted the source spectrum from a circular region centered on the source. For GIS2 and GIS3 we used the recommended 6 radius. For SIS0 and SIS1 we used 3 and 2 radii, respectively. Larger regions were not possible as the source was quite close to the chip gap. Background events were extracted from the the identical regions for each instruments from the deep background images supplied by the ASCA guest observed facility. For GIS, the latest (March 1995) response matrices were used, while for SIS, they were generated using FTOOL sisrmg. The spectra were then corrected for vignetting and the channels were grouped to contain at least 15 net counts per channel. In summary, the spectra for all instruments were extracted in a standard manner.

We also extracted background from different source-free regions on the same detector from our PHL 5200 observation. We found that the spectra extracted in this way did not alter our results in any significant way. For the rest of the paper we will discuss only the spectra with background extracted from the deep background observations.

Our observations were taken towards the end of ASCA mission, and the CCDs on the satellite were significantly degraded over time. There are problems with the charge transfer inefficiency (CTI), the residual dark distribution (RDD), and an unidentified additional loss of SIS low energy efficiency.111For ASCA Calibration Uncertainties see: http://legacy.gsfc.nasa.gov/docs/asca/cal\_probs.html For our observation in the 1CCD mode and for a weak source like PHL 5200, the latter is the most serious issue. The loss of the low energy response of the two SIS detectors can be parameterized with excess absorption (Yaqoob et al. 2000). This “NH correction” is different for SIS0 and SIS1 and is a function of time. We have taken into account all the calibration uncertainties carefully in our spectral extraction and analysis.

### 2.3 Spectral Analysis

The extracted spectra were analyzed using XSPEC version 11.0 (Arnaud, 1996). Data below 0.6 keV for SIS, below 0.9 keV for GIS and above 9.5 keV for all the instruments were ignored in spectral fits as these energy channels are not well calibrated. We first fitted the spectra from each instrument separately and found the results to be generally consistent. We then fitted all the spectra simultaneously leaving the normalizations free, to allow for small differences in absolute calibration.

As a first step we fitted the spectra with a simple power-law and Galactic absorption (with NH =$`4.8\times 10^{20}`$ cm-2, Stark et al. 1992). The “NH correction” discussed in section 2 was applied to the SIS spectra, which for our epoch of observation was $`6.8\times 10^{20}`$ cm-2 for SIS0 and $`9.8\times 10^{20}`$ cm-2 for SIS1, with $`\pm 10\%`$ uncertainty. The fit was not good, with $`\chi ^2=226.5`$ for 137 degrees of freedom. As excess absorption was indicated in the earlier spectrum of PHL 5200 (MES95) and in other BALQSOs (Mathur et al. 2000 and references there in), absorption at the source was added as a next step. The fit was better, but clearly showed negative residuals around one keV, and positive residuals at lower energies (Fig. 1, Fig. 2). Note that the behavior of the low energy residuals is exactly opposite to that due to the low energy calibration problem which produces spurious negative residuals. The implied recovery of the underlying power-law is a typical signature of partial covering by the absorber. So we fitted the spectrum with a partial covering model and the fit was significantly better ($`\chi ^2=182.5`$, 135 degrees of freedom). The resulting NH =$`5\pm 1\times 10^{23}`$ cm-2 and covering factor is $`0.9_{0.06}^{+0.05}`$. The confidence contours of the interesting parameters are plotted in figure 3. The scale over which figure 3 is plotted shows the range of $`3\sigma `$ uncertainty in the earlier observation (MES95). The column density at the source is greater than $`2\times 10^{23}`$ cm-2 and the power-law slope $`\alpha `$ is greater than 1.2 at 99% confidence (flux $`f(E)E^\alpha `$ where E is the energy).

The only residuals left after the partial covering model fit, are those above $`8`$ keV (observed frame), lying below the model values (Figure 2). Such residuals are not expected to be due to any calibration problems. They imply that a simple power-law continuum must turn over at higher energies. So we added a high energy cutoff to the model and re-fitted the data. This resulted in a much better fit with $`\chi _\nu ^2=1.18`$. This “best fit” model, however, resulted in an unusually low value of the high energy cut-off of E$`{}_{cuoff}{}^{}=18`$ keV (rest-frame) with e-folding energy of 0.4 eV. The typical high energy cutoff in quasar spectra is believed to be about 300 keV, though Yaqoob (2000) has drawn attention to the fact that this value is not well measured. We are aware of at least one another AGN with a very low cutoff value. In ESO 103-G35, Wilkes et al. (2001) found E$`{}_{cutoff}{}^{}=29\pm 10`$ keV. Low values of Ecuoff are likely to be important in understanding the spectrum of the cosmic X-ray background (Yaqoob, 2000). We note, however, that the negative residuals observed in our spectrum of PHL 5200 are significant only in the last few channels. So we caution against attaching too much significance to this result.

As discussed above, the source is highly absorbed with NH =$`5\times 10^{23}`$ cm-2. While the absorber is not Thomson thick, the effects of Thomson and Compton scattering may start becoming important at such column densities. In the analysis above, only photoelectric absorption is taken into account. The negative residuals observed at higher energies could be the effect of excess absorption due to Compton scattering, though for the observed column density it is not expected to be a significant effect (Matt et al. 1999). Nevertheless, we re-fitted our spectra with a partial covering model that fully incorporates Compton scattering in the absorber, as discussed in Matt et al. (1999). The resulting parameters are: NH =$`6\pm 2\times 10^{23}`$ cm-2, $`\alpha =1.9\pm 0.4`$, and covering fraction =$`0.94\pm 0.04`$. These are consistent with the model without Compton scattering, as expected.

For a spherical distribution of matter with NH =$`5\times 10^{23}`$ cm-2, a Fe-K emission line with equivalent width (EW) of 270 eV is expected (as measured against the observed continuum). We do not detect any line at the position of Fe-K$`\alpha `$ (rest energy 6.4 keV). The upper limit on the observed EW for a narrow line is 71 eV for SIS0 data and 102 eV for GIS3 data. The corresponding rest-frame EWs are 212 eV and 304 eV respectively. So we cannot put any meaningful constraint on the geometry of the absorber as viewed by the continuum source.

## 3 Discussion

The BALQSO PHL 5200 is clearly detected in all the instruments in this deep ASCA observation. With a total of 2610 “good” counts, this is the best X-ray spectrum of a BALQSO yet. The parameters of the spectral fit are well constrained; the X-ray flux is highly absorbed with N$`{}_{\mathrm{H}}{}^{}=5\pm 1\times 10^{23}`$ cm-2. While consistent with the conclusions of MES95, the earlier low signal to noise spectrum was also consistent with a model with no excess absorption (see also Gallagher et al. 1999). We unambiguously confirm large absorbing column density in PHL 5200. The X-ray weakness of BALQSOs thus is surely a result of absorption.

The absorber does not cover the continuum fully, with a covering fraction 90$`\pm 5`$%. This must be the reason for the apparent lack of excess absorption and flatter spectrum in the low signal to noise spectrum (see also Green et al. 2001). Partial covering of the absorber is consistent with the large optical polarization observed in this source (Goodrich & Miller 1995, Cohen et al. 1995) implying that there are at least two lines of sight to the nucleus: one direct, highly absorbed, and another scattered but unabsorbed. Gallagher et al. (1999) have discussed the possibility that partial covering by the absorber might be responsible for highly polarized BALQSOs to be relatively X-ray bright. Here we find that indeed, that is the best fit model to the PHL 5200 spectrum.

An associated absorber with high column density and partially covering the continuum is our preferred interpretation of the observed spectrum of PHL 5200. While physically motivated, this is not the only way to model the ASCA spectrum. For example, the low energy recovery of the absorbed spectrum may be be modeled as an unrelated soft-excess.222The nearby object, $`8^{\prime \prime }`$ away from the optical position of PHL 5200, is an early K star with no emission lines, and is therefore very unlikely to contribute to the X-ray flux of PHL 5200. The dip at $`1`$ keV is also a signature of a warm absorber. Note, however, that in PHL 5200 it corresponds to rest-frame energy of $`3`$ keV. Since absorption due to highly ionized sulphur or silicon occurs at these energies, such an interpretation might be plausible. However, it seems unlikely given that these are not among the most abundant elements.

The most surprising result of this observation is the steep spectrum, with X-ray power-law slope $`\alpha =1.7\pm 0.4`$. This might be an additional reason behind the unexpected non-detections of BALQSOs even in sensitive hard X-ray observations ($`\mathrm{\S }1`$). The mean ASCA slope for radio-quiet quasars at high redshifts is $`\alpha =0.67\pm 0.11`$, with a dispersion of $`\sigma =0.07`$ (Vignali et al. 1999, in the reshift interval between z=1.9 and 2.3). Our observations imply that the spectra of BALQSOs may not be exactly like other radio-quiet quasars, but are steeper. The observed, absorption corrected, 2–10 keV flux is $`3.7\times 10^{13}`$ ergs cm-2 s-1 (SIS0), consistent with the 1994 ASCA observation (MES95). So, the observed steep spectrum is unlikely to be a result of a short lived high-state. The only other BALQSO, PG 2112+059, for which a spectrum is available, has $`\alpha =0.98_{0.27}^{+0.4}`$ (Gallagher et al. 2001). As noted by Gallagher et al. , this is consistent with the mean quasar slope. However, it is also consistent with as steep as 1.38 at 90% confidence.

As discussed in section 2.2, CCDs on ASCA have been degraded over time. While we have been very careful in our analysis, there might be some unknown calibration uncertainties which have affected the results presented here. Future observations with Chandra and XMM-Newton will be useful in this respect. The discussion below assumes that these ASCA results are correct.

At low redshifts, Brandt, Mathur & Elvis (1997) found that the hard X-ray spectra of Seyfert galaxies are typically flatter than about $`\alpha =1.0`$. The Seyfert galaxies with steeper spectra are the narrow line Seyfert 1 galaxies (NLS1s). Brandt (2000) has discussed the analogy between low ionization BALQSOs and NLS1 galaxies. Mathur (2000) has discussed the analogy between BALQSOs and NLS1s further and has argued that NLS1 may be active galactic nuclei (AGNs) in the making. In this scenario, young radio-quiet AGNs are accreting at a high Eddington rate and have steep spectra. Over time, the accretion rate drops and the X-ray spectrum flattens. BALQSOs may be in that early evolutionary phase when the shroud surrounding the nuclear black-hole is being blown away and a quasar emerges (e.g. Fabian 1999, see also Becker et al. 2000). If the observed steep spectrum of PHL 5200 is a general property of BALQSOs, then it supports the evolutionary hypothesis of Mathur (2000) and further supports their analogy with NLS1 galaxies. On the other hand, PHL 5200 is one of the X-ray brightest BALQSOs, and so may not be representative of the BALQSO population. What is more likely, however, is that its orientation and cloud geometry fortuitously allow for substantial reflection of X-rays into our line of sight.

It has been an accepted paradigm that the orientation to the line of sight determines the appearance of an AGN. In addition to orientation, the age of a quasar also seems to play a role in its appearance. The steep X-ray spectra in NLS1s are a result of their high accretion rate, close to the Eddington limit (Pounds et al. 1995). BALQSOs, however, are much more luminous than NLS1s, and so must have far more massive black holes. Steep X-ray spectra in such systems may prove to be challenging to accretion theories.

We thank Tahir Yaqoob for help with ASCA calibration uncertainties, and Perry Berlind for obtaining the optical spectrum of the K star at the Fred Lawrence Whipple Observatory. This work is supported in part by NASA grants NAG5-9270 & NAG5-8913 (LTSA) to SM.


format: text

text:
# Free energy and molecular dynamics calculations for the cubic-tetragonal phase transition in zirconia

## I Introduction

A large class of advanced ceramics are solid solutions of zirconia (ZrO2) with cubic stabilising oxides like Y2O3, MgO or CeO, and are generally called stabilized zirconias. The long list of functional applications includes high-temperature devices, thermal barriers, and oxygen sensors. Moreover, partially stabilized zirconias represent a new generation of structural materials, by far the toughest ceramic oxides, strengthened by the mechanism called transformation toughening. The processing and service conditions of these materials involve phase transformations whose underlying physics is still a subject of controversy. One of these is the high-temperature cubic-tetragonal phase transition which is the subject of the present paper.

Zirconia is monoclinic ($`m`$) at low temperatures,  tetragonal ($`t`$) between 1400 and 2570 K,  and cubic ($`c`$) up to the melting point of 2980 K.  High-temperature X-ray experiments on stabilized zirconia revealed the existence of a $`ct`$ phase transition between 2300 and 2600 K,  depending on the atmosphere, but the mechanism of the transformation still has not been fully explained. The $`c`$ and $`t`$ unit cells are shown in Figure 1: note the characteristic tetragonal distortion of the oxygen sublattice in the $`t`$ phase.

It is not possible to quench to low temperature the $`c`$ and $`t`$ forms of pure zirconia, hence the experiments are difficult because of the high temperatures involved. Alternatively, the $`c`$ and $`t`$ structures may be stabilized at low temperatures by impurities. The available measurements are mostly done on stabilized samples. This simplifies the experimental procedure but complicates the interpretation of the results, because, besides the equilibrium $`t`$ phase, other metastable tetragonal structures are observed in stabilized crystals, denoted by $`t^{}`$ and $`t^{\prime \prime }`$. The former is the microstructure of a solid solution quenched from the field of stability of the $`c`$ phase into the biphasic $`c+t`$ one.  The $`t`$ and $`t^{}`$ forms are the same phase, they belong to the space group $`P4_2/nmc`$, but have different composition;  $`t^{}`$ is also called non transformable because it does not spontaneously transform to the $`m`$ phase. The $`t^{\prime \prime }`$ structure is observed in the ZrO2-ErO1.5  and ZrO2-Y2O3  systems, and has a cubic unit cell with the oxygen sublattice tetragonally distorted.

The microstructure of samples rapidly cooled from the $`c`$-phase region presents twinned domains separated by antiphase boundaries. The nature and composition of these domains is related to the phase transition mechanism and has been a subject of controversy. Originally they have been interpreted as the result of a diffusionless martensitic reaction.  Later, Heuer and Rühle  suggested that the transformation could be $`non`$-martensitic: homogeneous, massive, and displacive. Similarly, the observations of Lantieri et al.  were interpreted to mean that the $`ct^{}`$ transformation is diffusionless but non martensitic, and that the transformation always goes to completion. The same authors later proposed that the transition could be heterogeneous of the first-order with nucleation.  According to Sakuma,  the transition is instead second-order.

The temperature evolution of the tetragonality $`c/a`$ and of the anion sublattice distortion has been followed by Yashima et al.  in the ZrO2-ErO1.5 system. They showed that both order parameters depend continuously on the temperature and suggested that “the transition has the nature of a higher-order phase transition”.

Several attempts have been made in order to include this transformation in phenomenological theories. Hillert and Sakuma  expanded the free energy in terms of the defect concentration and assumed the transition to be second-order. Fan and Chen  used the time-dependent Ginzburg-Landau theory to expand the free energy of the transformation, treating it as a first-order one. The transformation was instead assumed to be second-order in the Landau energy expansion of Katamura and Sakuma.

The theoretical treatment of the $`ct`$ transition is simpler in stoichiometric zirconia: this is a case similar to the $`ct`$ phase transition in BaTiO3, where, according to symmetry considerations, the transformation could be either first or second order. A free energy Landau expansion for zirconia, involving the tetragonality of the cell only, without the distortion of the oxygen sublattice, inevitably predicts a first-order transformation.  But the inclusion of the latter in the Landau expansion opens the possibility for a second-order transition.  As already pointed out  the coupling between the order parameters may change the order of the transition from second to first.

In the case of BaTiO3 it has been possible to measure the order parameters very close to the transition temperature  and to establish the order and the mechanism of the transformation. Analogous experiments are difficult in pure zirconia because of the high transformation temperature. The neutron diffraction analysis of Aldebert and Traverse  provides the most complete thermomechanical description of pure $`c`$ and $`t`$ zirconia at high temperature. Aldebert and Traverse observe that: (i) The tetragonal distortion of the oxygen sublattice persists in the whole field of stability of the $`t`$ structure. (ii) The tetragonal distortion of the oxygen ions vanishes in the $`c`$ structure. (iii) The volume thermal expansion is linear and very close to isotropic up to near the transition point. As a consequence the $`c/a`$ ratio is almost temperature-independent over a wide range of temperatures, and sharply decreases near the transition temperature. (iv) The isotropic Debye-Waller factors of both species strongly increase before the transition temperature: the authors interpret it as a possible structural phenomenon anticipating the phase transition which could increase the ionic mobility.

The plan of the present paper is as follows. In Sec. II we introduce the main theoretical tools we have used to study the phase transition: the Landau theory of phase transformations, the thermodynamic integration technique, the constrained dynamics and the analysis of the order parameter fluctuations. The results are discussed in Sec. III. The phase transition mechanism was investigated using two sets of calculations. The first one, described in Sec. III A, is a traditional Molecular Dynamics (MD) analysis with which we observed the softening of a particular vibrational frequency. The second one, described in Sec. III B, combines the thermodynamic integration technique with constrained MD simulation to calculate the free energy surfaces around the phase transformation. We summarize the results in the final section.

## II Theory

### A Landau theory of the phase transition

#### 1 Order parameters

The Landau theory of phase transformations  describes the relationship between two crystal structures which share a common symmetry group $`𝑮`$0. The disappearance of a particular symmetry operation is quantitatively described by order parameters, which are zero in the high symmetry phase and become non-zero in the low symmetry one.

Our preliminary analysis  of the $`ct`$ phase transformation, based on 0 K calculations, showed that the transformation is driven by the distortion of the anion sublattice which is described by the primary order parameter $`\delta `$. This is a measure of the distance between each oxygen atom and the corresponding centrosymmetric position it occupied in the $`c`$ structure. The $`T=0`$ K calculations of certain phonon frequencies of the $`c`$ phase show that a frequency of vibration at the $`X`$ point of the BZ is imaginary.  This phonon, labelled $`X_2^{}`$, involves the oxygen sublattice only, and is shown in Figure 1. It transforms according to the $`A_2`$ irreducible representation of the little co-group of the $`X`$-point $`D_{4h}`$. The star of $`D_{4h}`$ contains three equivalent points, consequently the order parameter describing the tetragonal distortion has three components: $`\delta _x`$, $`\delta _y`$, and $`\delta _z`$.

In transforming to the $`t`$ phase, the primitive unit cell doubles, so that the phonon corresponding to $`X_2^{}`$ is at the $`\mathrm{\Gamma }`$ point, and is generally labelled $`A_{1g}`$. Nevertheless, in order to unify the description for both the $`c`$ and the $`t`$ structures, we will not use this convention and we will always refer to the soft mode as the $`X_2^{}`$ one, also in the $`t`$ phase.

Besides the tetragonal distortion of the oxygen sublattice, it is necessary to capture the change of the unit cell shape. This is done by introducing auxiliary order parameters $`\eta _i`$ defined in terms of the strain tensor $`\mathit{ϵ}`$. We decompose the six independent components of the strain tensor into irreducible representation of the $`O_h`$ cubic point group in Table I.

It was shown previously that, at equilibrium, each auxiliary order parameter is second-order in $`\delta `$ From now on, the order of the expansion terms will be expressed with respect to the order in $`\delta `$, therefore, as an example, a term like $`\delta _x^2\eta _1`$ is fourth order.

#### 2 Energy expansion

The Landau theory assumes that the appropriate thermodynamic potential of the crystal $`\mathrm{\Phi }`$ can be expanded in powers of the order parameters about the transition point. The Taylor expansion of $`\mathrm{\Phi }`$ must be invariant under the symmetry operations of the high symmetry phase. As a consequence, the allowed terms in the expansion have to be symmetry-invariants as well, and can be found using group theory. The terms in the energy expansion will be polynomials in the strains $`ϵ_i`$ and displacements $`\delta _i`$ of Table I. We constructed all the possible polynomials up to the sixth-order and symmetrized them with respect to the symmetry operations of the cubic point group $`O_h`$. The resulting invariants are shown in Table II.

This analysis showed that all the third-order invariants are identically zero, which is a necessary (but not sufficient) condition for a phase transition to be second-order. We already mentioned that the instability appears at the boundary of the BZ, therefore it halves the number of symmetry elements, and this is a further condition allowing the $`ct`$ phase transition to be second-order.

In order to keep the discussion as simple as possible, at this stage we assume the phase transition to be second-order, truncating the Taylor expansion at the fourth-order term in $`\delta `$. The possible importance of the higher order terms will be discussed later. The energy expansion, expressed in terms of the basis function defined in Table II, is as follows:

$`F`$ $`=`$ $`F_0+{\displaystyle \frac{a_2}{2}}A_2\left(\delta ^2\right)+{\displaystyle \underset{i=1}{\overset{2}{}}}{\displaystyle \frac{a_{4i}}{4}}A_{4i}\left(\delta ^4\right)`$ (2)
$`+{\displaystyle \underset{i=1}{\overset{3}{}}}b_iB_i(ϵ,\delta ^2)+{\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle \frac{c_i}{2}}C_i\left(ϵ^2\right)+𝒪\left(\delta ^6\right).`$
$`F_0`$ is the energy of the high-symmetry phase and is a function of the hydrostatic strain $`\eta _1=Tr(`$$`\mathit{ϵ}`$$`)`$. The choice of the reference volume fixes $`F_0`$ and the expansion coefficients $`a_2,\mathrm{},c_3`$. In the present case, the energy $`F`$ was expanded about the minimum of the energy-volume curve for the $`c`$ structure predicted by the SC-TB calculations.

### B Free energy calculation

Free energy surfaces may be calculated directly from MD simulations in terms of ensemble averages by using the thermodynamic integration technique.  Here we briefly describe how this method was applied to zirconia.

The thermodynamic integration method allows us to calculate free energy differences between a reference state, for which the internal energy $`U_0`$ is known, and another state of the same system, with internal energy $`U`$. The idea is to relate the two structures with a switching parameter $`\gamma `$ which is zero in the reference state and non-zero otherwise. The free energy variation in the infinitesimal change $`d\gamma `$ may be calculated using standard statistical mechanics:

$$dF=\frac{U(\gamma )}{\gamma }_{\overline{\gamma }}d\gamma $$
(3)
This is equivalent to the reversible work done for the structural modification described by $`d\gamma `$, implicitly assumed to be adiabatic. By $`\mathrm{}`$ we indicate the ensemble average, which has to be calculated at a constant value of $`\gamma =\overline{\gamma }`$. The free energy difference can be obtained by integrating the previous equation. In the general case, $`U(\gamma )`$ is not known. The common strategy is to perform several constrained MD simulations at different values of $`\gamma `$ and then integrate Eq.(3) numerically. Many calculations may be necessary in order to integrate Eq.(3) with sufficient precision.

A knowledge of the functional form of the energy would greatly simplify this procedure, reducing the number of calculations and allowing the analytic integration of Eq.(3). The Landau theory in combination with MD simulations can provide such useful information. In order to apply this formalism, it is necessary to define the thermodynamic variables of Eq.(2) from a MD run at finite temperature. Statistical mechanics allows us to calculate the order parameters by averaging the corresponding time-dependent ones over all the available atomic configurations.

The primitive $`c`$ cell is unstable with respect to 3 modes of vibration whose frequency is degenerate. The instability appears at the $`X`$, $`Y`$, and $`Z`$ points of the primitive BZ and the corresponding eigenmodes distort the anion sublattice along the $`x,y`$, and $`z`$ directions. In the following we consider a supercell which is not the primitive one, and those points, originally at the border of the BZ, are folded in at the $`\mathrm{\Gamma }`$ point. The eigenvectors are therefore real.

Let us denote by $`𝐮`$ the atomic displacements from a perfect site of the high symmetry phase. We expand $`𝐮`$ in normal coordinates using the notation of Maradudin et al. :

$$𝐮(\kappa )=\frac{1}{\sqrt{M_\kappa }}\underset{j}{}𝐞(\kappa |j)Q(j).$$
(4)
$`\kappa `$ and $`M_\kappa `$ label the atoms and their mass in the cell, $`𝐞(\kappa |j)`$ is the eigenvector $`j`$ at the $`\mathrm{\Gamma }`$ point of the BZ, and $`Q(j)`$ is the corresponding normal coordinate. We will denote by $`|\alpha )`$, where $`\alpha =x,y,z`$, the indexes $`j`$ describing the soft modes.

Given a general atomic configuration at time $`t`$ (we now include the time in the notation), we define the time-dependent order parameter as the average displacement along $`X_2^{}`$ of the $`r_\mathrm{O}`$ oxygen atoms of the cell:

$`\delta _\alpha (t)`$ $``$ $`{\displaystyle \frac{1}{\sqrt{r_\mathrm{O}}}}{\displaystyle \underset{\kappa }{}}𝐮(\kappa ,t)𝐞(\kappa |\alpha )`$ (5)
$``$ $`{\displaystyle \frac{Q(\alpha ,t)}{\sqrt{r_\mathrm{O}M_\mathrm{O}}}}.`$ (6)
The time averages of these quantities, $`\overline{\delta }_\alpha `$, are the experimentally measurable order parameters which we now take as thermodynamic variables.

The factor $`U(t)/\delta _\alpha `$ entering in Eq.(3) can now be calculated at each time step by applying the definition of $`\delta `$ given in Eq.(5), and by using the chain rule:

$$\frac{U(t)}{\delta _\alpha (t)}=\sqrt{r_\mathrm{O}M_\mathrm{O}}\underset{\kappa }{}\left(\frac{U(t)}{𝐮(\kappa ,t)}\right)\frac{𝐮(\kappa ,t)}{Q(\alpha ,t)}.$$
(7)
Noting that the eigenvectors are orthonormal and that the first term of the sum is the force $`𝐅`$ acting on the atoms we end up with the following expression:

$$\frac{U(t)}{\delta _\alpha (t)}=\sqrt{r_\mathrm{O}}\underset{\kappa }{}𝐅(\kappa ,t)𝐞(\kappa |\alpha ).$$
(8)
Therefore the free energy gradient is calculated from the time average of the atomic forces projected along the $`X_2^{}`$ mode of vibration:

$$\frac{dF}{d\delta _\alpha }=\sqrt{r_\mathrm{O}}\underset{\kappa }{}𝐅(\kappa ,t)𝐞(\kappa |\alpha )_{\overline{\delta }_\alpha }.$$
(9)
Note that the above average has to be taken on an ensemble with a constant value of order parameter $`\overline{\delta }_\alpha `$, i.e. it is necessary to constrain the order parameters during the MD simulations.

### C Constraining the order parameters

The dynamics of canonical and microcanonical ensembles with fixed cell shape automatically constrain the auxiliary order parameters. On the contrary, in order to constrain the dynamics of the primary order parameters and then integrate Eq.(9) from the results of the MD simulation, it is necessary to modify the Lagrangian of the system.

The goal is to obtain an equation of motion describing the time evolution of a system with a fixed order parameter $`\overline{\delta }`$. This is done by extending the Lagrangian of the unconstrained system $`^u`$:

$$=^u\underset{\alpha }{}\lambda _\alpha \sigma _\alpha .$$
(10)
The superscript $`u`$ stands for unconstrained, the $`\lambda `$’s are the Lagrange multipliers to be calculated and the $`\sigma `$’s are the functions describing the constraints. Three of them are needed, one for each direction $`\alpha `$ of the tetragonal distortion:

$$\sigma _\alpha (t)=\delta _\alpha (t)\overline{\delta }_\alpha =0.$$
(11)
The Lagrangian of the constrained system is obtained from equations (10) and (11), and the corresponding equations of motion are:

$$M_\mathrm{O}\ddot{𝐮}_\alpha (\kappa ,t)=𝐅_\alpha (\kappa ,t)\frac{\lambda _\alpha (t)}{\sqrt{r_\mathrm{O}}}𝐞(\kappa |\alpha ).$$
(12)
Where $`𝐮_\alpha `$ and $`𝐅_\alpha `$ are the $`\alpha =x,y,z`$ components of the displacement and of the force. The orthonormality of the normal modes of vibration decouples the equations along the three crystallographic directions, simplifying the implementation of the method. Moreover, since the tetragonal distortion involves the anion sublattice only, we need apply the above modified equation only to the $`r_\mathrm{O}`$ oxygen atoms. In general, the Lagrange multipliers have to be found numerically, but in the present case (decoupled crystallographic directions and linear constraints) an analytical solution does exist.

The expression of the Lagrange multipliers may be found as follow: (i) Advance the atomic positions with a fake unconstrained MD step. (ii) Use these unconstrained coordinates to find the multipliers which exactly satisfy the constraints. (iii) Use these values of $`\lambda `$’s to perform the true MD step which satisfies the constraining equations by construction. Here we specify this procedure for the leapfrog Verlet algorithm.

Given a set of atomic positions $`𝐮(t)`$ which satisfy the constraining equations (11), the fake step involves solving the equation of motion corresponding to the Lagrangian $`^u`$. By doing so, the set of unconstrained coordinates $`𝐮^u(t+\mathrm{\Delta }t)`$ is obtained. These are related to the constrained atomic positions $`𝐮(t+\mathrm{\Delta }t)`$ as follows:

$$𝐮_\alpha (\kappa ,t+\mathrm{\Delta }t)=𝐮_\alpha ^u(\kappa ,t+\mathrm{\Delta }t)\frac{\lambda _\alpha (t)\mathrm{\Delta }t^2}{\sqrt{r_\mathrm{O}}M_\mathrm{O}}𝐞(\kappa |\alpha ).$$
(13)
Applying the definition (5) to these coordinates, a similar relationship may be found for the order parameters.

$$\delta _\alpha (t+\mathrm{\Delta }t)=\delta _\alpha ^u(t+\mathrm{\Delta }t)\frac{\lambda _\alpha (t)\mathrm{\Delta }t^2}{r_\mathrm{O}M_\mathrm{O}}$$
(14)
The analytic solution of the Lagrange multipliers is obtained by imposing the constraining equations $`\sigma (t+\mathrm{\Delta }t)=0`$ and then solving the resulting linear equation in $`\lambda `$:

$$\lambda _\alpha (t)=\frac{r_\mathrm{O}M_\mathrm{O}}{\mathrm{\Delta }t^2}\left[\delta _\alpha ^u(t+\mathrm{\Delta }t)\overline{\delta }_\alpha \right].$$
(15)
The substitution of (15) in (13) gives the constrained coordinates at $`t+\mathrm{\Delta }t`$ in the NVE ensemble:

$`𝐮_\alpha (\kappa ,t+\mathrm{\Delta }t)`$ $`=`$ $`𝐮_\alpha ^u(\kappa ,t+\mathrm{\Delta }t)`$ (17)
$`\sqrt{r_\mathrm{O}}\left[\delta _\alpha ^u(t+\mathrm{\Delta }t)\overline{\delta }_\alpha \right]𝐞(\kappa |\alpha ).`$
It is important to notice that using this method, the expressions for the multipliers are functions of both the integrating scheme and any other additional constraints, such as thermostats. The simple case of the leapfrog Verlet algorithm described here has to be slightly modified in order to include the Nosé-Hoover thermostat.  The same procedure may be repeated for the NVT ensemble and the resulting equations of motion are:

$`\ddot{𝐮}_\alpha (\kappa ,t)`$ $`=`$ $`{\displaystyle \frac{𝐅_\alpha (\kappa ,t)}{M_O}}`$ (19)
$`\sqrt{r_\mathrm{O}}\left[{\displaystyle \frac{\delta _\alpha (t+\mathrm{\Delta }t)\overline{\delta }_\alpha }{\mathrm{\Delta }t^2}}\left(1+\xi (t){\displaystyle \frac{\mathrm{\Delta }t}{2}}\right)\right]𝐞(\kappa |\alpha ),`$
where $`\xi `$ is the thermostat variable.

### D Fluctuations

The fluctuations of the instantaneous order parameter $`\delta _\alpha (t)`$ were used to calculate the frequency of a particular vibration directly from the MD run. The central point of this analysis is the calculation of the fluctuation correlation function spectrum:

$$S_\alpha (\nu )=e^{i2\pi \nu t}\delta _\alpha (t=0)\delta _\alpha (t)𝑑t.$$
(20)
The above dynamic form factor is known to exhibit two important features,  a temperature-dependent resonant peak at $`\overline{\nu }`$, and an additional central peak at $`\nu =0`$. The relative magnitude of the two peaks depends on the transformation mechanism and on the temperature. This has been proved for phase transition mechanisms as different as order-disorder and displacive.  Therefore, without loss of generality, following Padlewski et al. , the power spectrum (20) can be modelled as a superposition of two peaks with the following functional form:

$$S_\alpha (\nu )=\frac{2AB}{B^2+\nu ^2}+\frac{CD}{D^2+\left(\nu \overline{\nu }_\alpha \right)^2},$$
(21)
where $`A`$, $`B`$, $`C`$, $`D`$ and $`\overline{\nu }_\alpha `$ are parameters to be fitted to the calculations. The analytical form of the time-dependent correlation function $`S(t)=\delta _\alpha (t=0)\delta _\alpha (t)`$ may be found by substituting (21) in (20) and passing into the time domain with an inverse Fourier transform:

$$S_\alpha (t)=Ae^{Bt}+Ce^{Dt}\mathrm{cos}(2\pi \overline{\nu }_\alpha t).$$
(22)
The time-dependent order parameter is calculated from the MD atomic positions. The time correlation function of $`\delta _\alpha (t)`$ is then obtained using the multiple time-origin method  and fitted to Eq.(22). The fitting procedure provides both the time correlation function and its Fourier transform.

## III MD Simulations

The polarizable self-consistent Tight-Binding model  was used to perform two sets of MD simulations. In the first one, standard MD calculations were used to investigate the temperature dependence of the order parameters and to follow the softening of the $`X_2^{}`$ mode of vibration up to the transition point. In the second one, we combined the constrained MD simulations and the thermodynamic integration technique to calculate the free energy of the $`ct`$ phase transition, to study the nature of the order parameter fluctuations and to explain the high temperature stability of the $`c`$ phase.

### A Standard MD simulations

#### 1 Softening of a vibrational frequency

The time evolution of a system of 96 particles with periodic boundary conditions has been followed at temperatures between 300 and 2200 K. The lattice parameter of the simulation cell were the correspondent experimental values of Aldebert and Traverse, which, where necessary, have been linearly extrapolated at lower temperatures. During each MD run, the temperature has been constrained with a Nosé-Hoover thermostat,  and the equations of motion have been integrated for not less then 5 ps with a typical time step of 5 fs. Near the transition point, the time step has been reduced to 2.5 fs and the total simulation time has been increased to 15 ps.

The cell size was constrained by the relatively high number of MD runs necessary to follow the phase transition. In total, we simulated the time evolution of 96 particles for more than 120 ps. As discussed later, the cell size does not change our qualitative description of the phase transition, and the 324-atom unit cell would have just implied a heavier computational effort, without adding further information to the physical picture provided by the smaller cell.

We started the simulations from the crystallographic positions of the tetragonal phase and equilibrated the system at the temperature of 300 K. This temperature is well inside the field of stability of the $`m`$ phase, however, during the MD simulations, the system remained in the $`t`$ phase because of the existence of an energy barrier between the two structures. We calculated the vibrational frequencies of the $`t`$ structure by diagonalising the dynamical matrix at the origin and at the borders of the Brillouin Zone (BZ) along the (100), (110), and (111) directions. This analysis showed that all the vibrational frequencies are real and that the $`t`$ phase does not spontaneously distort towards the $`m`$ structure.

In this set of MD simulations we followed the approach of Padlewski et al.  described in Section II D, focusing on the instantaneous order parameters $`\delta _\alpha (t)`$ \[see Eq. (5)\] which fluctuate about the mean value $`\overline{\delta }_\alpha `$. Figure 2 shows the typical time evolution of the primary order parameters for the MD run at T=700 K. Figure 3 shows the fluctuation autocorrelation function $`S(\nu )`$ and the corresponding frequency spectrum $`S(t)`$ for the MD run at 700 K: the arrow points at $`\overline{\nu }_z`$, the frequency which softens. It can be seen that the $`x`$ and $`y`$ components, corresponding to the transverse optical frequencies, are degenerate.

On increasing the temperature, the softening of the frequency $`\overline{\nu }_z`$ is evident from the dynamic form factor, where the resonant peak shifts. At the same time, the primary order parameter decreases continuously (Figure 4), as experimentally observed in the similar system ZrO2-12%ErO1.5 The calculated temperature dependence of the macroscopic order parameter $`\overline{\delta }_z`$ and of the corresponding vibrational frequency, shown in Figure 5, was then interpreted using the Landau theory. We found that the critical exponent for this phase transition is $`\beta `$=0.35. According to the same theory, the critical exponent $`\beta ^{}`$ for the auxiliary order parameters ($`\eta _2,\eta _3`$) describing the tetragonality of the cell is bigger then $`\beta `$. Therefore the $`c/a`$ ratio should depend more strongly on the temperature than $`\delta `$.

As the transition temperature is approached, the decrease of the order parameter and of the corresponding frequency is accompanied by an increase in the order parameter fluctuations, which theoretically diverge at $`T_c`$ for a second order phase transition. As a result, it was not possible to follow the complete softening of the frequency: there is a temperature window about $`T_c`$ where, even though long MD simulations allow one to evaluate the average order parameters, it is not possible to calculate the frequency. In this temperature range, the frequency $`\overline{\nu }`$ is so low that the corresponding peak in the dynamic form factor $`S(\nu )`$ merges with the central peak and it is not possible to separate them.

The theoretical transition temperature of $``$1800 K is $``$30% lower than the experimental value of $``$2600.  This may be explained by noting that the first-principles calculations underestimate the energy difference between the $`c`$ and $`t`$ structures $`\mathrm{\Delta }U^{tc}`$, which determines $`T_c`$. This is the ab initio energy barrier between the minima of the double well, which was used to parametrize the SC-TB model. In particular the SC-TB results underestimate the experimental $`\mathrm{\Delta }U^{tc}`$ by 30%, which is consistent with the underestimate of the transition temperature.

According to the renormalized phonon group theory,  $`\overline{\nu }^2`$ depends linearly on the temperature in both the regions $`T<T_c`$ and $`T>T_c`$, and the correspondent slopes are related by the following relationship:

$$R=\left(\frac{d\overline{\nu }^2}{dT}\right)_{T<T_c}/\left(\frac{d\overline{\nu }^2}{dT}\right)_{T>T_c}=2.$$
(23)
However, our simulations at $`T>T_c`$ suggest that the $`ct`$ phase transition in zirconia has a different behaviour from the ideal case described by Eq.(23), because no frequency was observed above $`T_c`$.

The exploration of the high temperature region of the zirconia phase diagram has been carried out in two stages. As a first attempt, we continued the MD simulations on the system described above, simply increasing the temperature. This has been done up to 2200 K. The time autocorrelation function (22) of these simulations exponentially decayed without showing any structure. As a result, the central peak dominated the corresponding dynamic form factor, and therefore it was not possible to isolate the resonant peak at $`\overline{\nu }`$ from the central one. A possible explanation of this may be proposed by noticing that, according to Eq.(23), for $`T>T_c`$ the slope of $`(\frac{d\overline{\nu }^2}{dT})`$ is half that for $`T<T_c`$. This means that the temperature window around $`T_c`$, in which it is not possible to calculate the frequency, extends more in the high temperature field than in the low temperature one. Probably 2200 K is still in the region of disturbance of the transition point.

In order to verify if the frequency does eventually increase in the high temperature region, we studied a similar system with the same properties of that one described above but with a lower transition temperature. The idea is based on the following argument. It is well established that the relative energetics of the two phases is governed by a double well in the potential energy that depends on volume and $`c/a`$ We studied in detail its dependence,  which is also captured by the Landau expansion (2). Both the hydrostatic and tetragonal strains modify the double well in the same way: the smaller the volume (or the $`c/a`$ ratio), the smaller the energy difference and therefore the smaller the transition temperature. Incidentally, this is connected to the ab initio underestimate of the energy barrier, which is calculated with the structural parameters corresponding to 0 K.

By exploiting this property of the energy surfaces, we made a new set of MD simulations aimed to explore the temperature range $`T>T_c`$. In these calculations the volume was chosen to lower the transition temperature to $``$1300 K and the cell was kept cubic ($`c/a=1`$) even in the low temperature region. As expected, for $`T<T_c`$, the linear softening of $`\overline{\nu }^2`$ (Figure 5 set b) was obtained for this system as well. The slope was slightly different because in the previous simulations the thermal expansion of the cell was included in the description, while in this case the volume was fixed to the initial value. Because of this we could calculate the frequency up to within $``$200 K of the transition temperature. The temperature was then increased to 2000 K. Surprisingly, even in this case, there was no structure in the autocorrelation function $`S(t)`$ and the expected hardening of the frequency was not observed. This suggests that in the $`c`$ phase the motion of the oxygen sublattice along the $`X_2^{}`$ mode of vibration is, in terms of $`S(t)`$, uncorrelated. This behaviour will be clarified by the free energy surfaces described in Section III B.

### B MD simulations at constant $`𝜹`$

In our previous papers  we restricted the analysis of the 0 K energy surface to one tetragonal invariant only. By doing so, we defined a simplified version of the energy expansion (2) involving the strain and one component of the primary order parameter. We then fitted the correspondent coefficients, $`a_2`$, $`a_{41}`$, $`b_1`$, $`b_2`$, $`c_1`$, and $`c_2`$, to the results of total energy calculations. We also showed how the coupling between the primary and auxiliary order parameters could create a critical point where the transformation becomes first order. Here that analysis is extended by exploring the topology of the energy surface in the whole $`\delta `$ domain, and by following its temperature evolution through the phase transition into the field of stability of the $`c`$ phase. This sheds light on the mechanism of the phase transformation and on the high-temperature stability of the $`c`$ phase.

The following results were obtained using a 12-atom unit cell with different $`c/a`$ ratios (1, 1.01, 1.02) at the 0 K theoretical equilibrium volume of the $`c`$ structure. Preliminary unconstrained MD simulations were done to explore the effect of the cell size on the physical picture of the phase transformation described in the previous Section. Even in this small system, the frequency $`\overline{\nu }_z^2`$ depends linearly on the temperature and the predicted transition temperature is of $`1600`$ K. The effect of using a small cell is to shift $`T_c`$ to higher temperatures. This is consistent with the physical picture proposed in Section III A: the autocorrelation function $`S(t)`$ measures the degree of correlation between the motion of the oxygen atoms along the $`X_2^{}`$ mode of vibration. We described how the temperature acts on $`S(t)`$ by reducing the correlation until this is completely lost above $`T_c`$, where the corresponding frequency is soft, and where the structure is $`c`$. The small cell size and the periodic boundary conditions force the motion of atoms in adjacent cells to be correlated, and therefore counteract the effect of the temperature on $`S(t)`$. As a result, in the small system, higher temperatures are needed to observe the complete softening of the frequency $`\overline{\nu }_z`$.

The 12-atom and 96-atom supercells have the same temperature dependence of $`\delta _z`$ and $`\overline{\nu }_z`$ but the corresponding curves are shifted to different temperatures. We can therefore conclude that the phase transformation mechanism is the same in the two systems. We shall calculate the free energy of the transition from the MD simulations of the small cell, and assume that the resulting qualitative physical picture applies to bigger cell sizes.

Before exploring the free energy temperature dependence, it is useful to simplify the complete energy expansion (2) by neglecting the order parameters which are unlikely to play an important role in the phase transition. The transformation between the $`c`$ and $`t`$ structures does not distort the cell shape as described by the order parameters ($`\eta _4,\eta _5,\eta _6`$). It is therefore reasonable to neglect them in the discussion of the following results. Moreover, even though the transformation between the $`c`$ and the $`t`$ structure does involve a change in the volume, the energetic contribution of the associated order parameter $`\eta _1`$ is well understood and has already been discussed. Apart from the 0 K case, we will not consider the terms $`B_3`$ and $`C_3`$ in the energy expansion. However, their possible influence on the character of the phase transition in terms of softening of the corresponding elastic constant, will be discussed a posteriori in the final Section III B 3.

#### 1 Topology of the 0 K surface

We start our analysis with the primary order parameter. Two sets of calculations on a stress-free cubic unit cell were used to fit the coefficients $`a_2`$, $`a_{41}`$, and $`a_{42}`$. These have been determined by distorting the oxygen sublattice along $`\delta 00`$ and along $`\delta \delta 0`$. We plot the resulting energy surface, which we take as starting point of our analysis, as a function of two tetragonal invariants in Figure 6. In this simple case, because of the cubic cell, the three components of the primary order parameter are equivalent.

The same set of calculations was then done on a tetragonal cell ($`c/a=1.01`$), by which we determined the parameters $`b_2`$ and $`c_2`$. The latter is proportional to the elastic constant $`C^{}`$. The transferability of the parameters was then checked by redoing the calculations for a different tetragonal cell ($`c/a=1.02`$): Figure 8 shows that the same set of coefficients fit the results for this cell as well. If $`z`$ is the tetragonal axis, the tetragonality of the unit cell shortens the average interatomic distances in the transverse $`x,y`$ plane and lengthens them along the tetragonal axis. As a consequence, the energy surface section in the transverse plane $`\delta _x,\delta _y`$, shown in Figure 7 (a), is similar to the reference one of Figure 6 but shallower and tighter, while it is deeper and broader along $`\delta _z`$ \[Figure 7 (b)\].

Finally, the same procedure described above was used to fit the remaining coefficients $`b_1`$ and $`c_1`$ by distorting the cell with respect to the order parameter $`\eta _1`$ defined in Table I.

In conclusion, the static calculations show that the 0 K energy surfaces can be captured by Taylor expansion up to fourth-order and are therefore completely defined by the set of coefficients given in Table III.

#### 2 Free energy surfaces

The MD simulations were carried out in the temperature range from 50 K to 2000 K, constraining the primary and secondary order parameters. Let us first focus on the results for the cubic cell, commenting later on the effect of the $`c/a`$ ratio. The explorations along the directions $`\delta 00`$ and $`\delta \delta 0`$ fully determine the free energy surfaces to fourth order, therefore we constrained the order parameters along these directions from 0 to 0.7 a.u., using the dynamics described in Section II C. The quantity defined in Eq.(8) was accumulated during the MD run and its time average provided the ensemble average required in Eq.(9). The analytical form (2) of the energy surface was then differentiated along the corresponding direction and fitted to the results of the simulations. For this particular case, the fit provided both the free energy gradient and the free energy itself. This is because we chose the reference energy as the top of the double well for a cubic crystal. The integration of Eq.(9) provides the energy difference $`\mathrm{\Delta }F`$:

$$\mathrm{\Delta }F=F(ϵ,\overline{\delta })F(ϵ,\overline{\delta }=0)=_0^{\overline{\delta }}\frac{U}{\delta }_{\overline{\delta }}𝑑\delta .$$
(24)
We arbitrarily set to zero the integration constant for the cubic cell $`F(ϵ=0,\overline{\delta }=0)`$. If the cell is tetragonal the same constant is $`c_2C_2(ϵ^2)/2`$. Therefore, the fit of the free energy gradient to the MD simulations of the cubic cell provided the coefficients $`a_2`$, $`a_{41}`$, $`a_{42}`$ and $`b_2`$ at the corresponding temperature.

The fit to the computed results along the $`\delta 00`$ and the corresponding free energy profiles obtained from Eq.(24) are shown in Figures 9 and 10 respectively. The corresponding expansion coefficients are included in Table III. Also at high temperature the fourth order energy expansion well describes the calculated free energy gradients. The temperature acts on the double well by gradually reducing the energy difference between the distorted and undistorted structure. As expected, the energy surface is very flat near the critical temperature, but, surprisingly, it remains quite flat even at higher temperatures, well inside the field of stability of the cubic phase. The energy surfaces above the transition temperature are highly anharmonic.

The unconstrained MD simulations described in the previous Section generated the temperature dependence of the $`X_2^{}`$ vibration frequency (Figure 5) up to the transition point only. Above the critical temperature it was impossible to calculate the frequency $`\overline{\nu }_z`$, the dynamic form factor $`S(\nu )`$ being dominated by a wide central peak. The soft frequency $`\overline{\nu }_z`$, together with the large fluctuations of the order parameters $`\delta _z`$, suggests a disordered dynamics of the oxygen sublattice along the $`X_2^{}`$ mode of vibration. The anharmonicity of the energy surfaces explains this behaviour. It may also explain the fact that the optical-phonon branches are not experimentally observed  in cubic stabilized zirconia.

As in other perfect fluorite structures,  the vibrational motion of the anions in $`c`$ zirconia appears to be anharmonic. This is consistent with the neutron powder diffraction experiments of Kisi and Yuxiang  on cubic stabilized zirconia (ZrO2-9.4%Y2O3). They measured the temperature dependence of the Debye-Waller factor and proposed different models to fit the data. Both a simple Debye model and a Debye model plus a static disorder component provided a poor fit of the data. A radical improvement of the fit was obtained when an isotropic anharmonic vibration of both species was included in the description.

We may further investigate the nature of the double well temperature dependence by splitting the free energy into its energetic and entropic contributions. The time average of the internal energy $`U`$ from each constrained MD simulations is plotted in Figure 11. It is clear that the double well in the internal energy is present even at $`T>T_c`$. The internal energy double well is relatively insensitive to the temperature and the 0 K coefficients of Table III provide an excellent fit of the finite temperature results for both the low and high temperature structures. From the definition of $`F`$, the difference between the free energy calculated with Eq.(24) and the internal energy gives the entropic contribution $`TS`$ plotted in Figure 12. The high temperature stability of the $`c`$ phase is therefore ensured by this entropic term which change the shape of the energy surface from double to single welled.

#### 3 Coupling to the elastic strains

The fit of the Landau energy expansion (2) to the calculation results has allowed us to follow the temperature evolution of the free energy surface through the phase transition. The fitting coefficients for each temperature are included in Table III. We can see that at the critical temperature of $`1600`$ K, when $`a_2`$ goes to zero, the fourth order coefficient $`a_4`$ is positive. Therefore, in this cell, the phase transition is diffusionless displacive and second order.

If the thermodynamic potential $`F`$ is expanded in terms of elastic strains as well as $`\delta `$, the coupling between the primary order parameter and the strain in (2) renormalizes the fourth-order coefficient  and could make it small or negative near the transition point. As a result, the transformation may become first-order. This could happen if the temperature reduces an elastic constant $`c_i`$.

In order to see this, let us consider a tetragonal cell at the reference volume, whose tetragonal axis is $`z`$ and where the oxygen sublattice is distorted along the order parameter direction $`00\delta _z`$. With these restrictions the energy expansion (2) has a simplified form:

$$F=F_0+\frac{a_2}{2}\delta _z^2+\frac{a_{41}}{4}\delta _z^4+b_2\mathrm{\hspace{0.17em}\hspace{0.17em}2}\delta _z^2\eta _2+\frac{c_2}{2}\eta _2^2+𝒪\left(\delta _L^6\right).$$
(25)
At equilibrium the two order parameters $`\delta _z`$ and $`\eta _2`$ are not independent and the relationship between the two may be found by imposing the equilibrium condition:

$$\frac{F}{\eta _2}=0\eta _2=\frac{2b_2}{c_2}\delta _z^2.$$
(26)
A similar procedure may be repeated for the hydrostatic strain $`ϵ_{xx}+ϵ_{yy}+ϵ_{zz}=\eta _1`$ and the substitution of Eq.(26) back in Eq.(25), together with the corresponding one for $`\eta _1`$, clarifies the combined effect of the coefficients $`b_i`$ and $`c_i`$ on the renormalization of the fourth order coefficient:

$$F=F_0+\frac{a_2}{2}\delta _z^2+\left(\frac{a_{41}}{4}\frac{2b_2^2}{c_2}\frac{b_1^2}{2c_1}\right)\delta _z^4+𝒪\left(\delta _z^6\right).$$
(27)
It may be verified that the above formulation is independent of the initial choice of the tetragonal axis.

The coefficients $`c_1`$ and $`c_2`$ are proportional to the bulk modulus and to the elastic constant $`C^{}`$, respectively. If one or both of these quantities significantly decrease with temperature, the correspondent term in Eq.(27) may dominate the sign of the fourth-order coefficient, making it very small or negative. At 0 K the coupling terms $`2b_2^2/c_2`$ and $`b_1^2/2c_1`$ reduce the fourth-order coefficient by $`16`$% and $`4`$% respectively. Because of this difference we focused our attention on the elastic constant $`C^{}`$. Experimentally, the high temperature data  do not show any anomalous temperature dependence of the elastic constants, with a general decrease of $`1520`$% between 300 and 1700 K. If this is valid for pure zirconia as well, we may anticipate that the degree of softening of the elastic constants will not affect the character of the phase transition.

The calculations described above for the cubic cell have been repeated for two tetragonal cells with $`c/a=1.01`$ and $`1.02`$. The constrained MD simulations along the $`\delta 00`$ direction allowed the fit of the coupling coefficient $`b_2`$ to each temperature. Therefore it has been possible to obtain the temperature and $`c/a`$ dependent free energy curve along this direction. However, in this case we are interested in the relative position of the free energy surfaces for the different cells, depending on the integration constant $`F(ϵ,\delta =0)=c_2C_2(ϵ^2)/2`$. In principle one could calculate $`c_2`$ by a simple thermodynamic integration, but this would require monitoring the stress, which is not yet implemented in the current program. Therefore we decided to continue the analysis by choosing an extreme scenario, which would be a strong temperature dependence of $`C^{}`$. With a large safety margin with respect to the experimental values, we linearly reduced this elastic constant by $`75`$ % from 200 to 50 MPa between 0 and 2000 K.

With this assumption, the sections of the free energy surfaces at different temperatures between 500 and 2000 K are plotted in Figure 13 as a function of $`\delta `$ and $`c/a`$. As experimentally observed,  the tetragonal cell with $`c/a=1.02`$ is thermodynamically stable up to very near $`T_c`$ where the minimum configuration goes from $`(c/a=1.02,\delta =0.4)`$ to $`(c/a=1,\delta =0)`$ quickly but continuously. We believe that the sudden change of the order parameters, due to the flat energy surfaces near the transition point and therefore to the anharmonicity of the material, may explain the fact that this phase transition has been considered to be first order in the early studies.

It is interesting to note that both the $`a_{41}`$ and the coupling coefficient $`b_2`$ decrease with temperature (Table III). Because of this, even with the large postulated softening of $`C^{}`$, the renormalization of the fourth-order coefficient in Eq.(27) does not increase with temperature. From the data of Table III, we find that at $`T_c`$, the term $`2b_2^2/c_2`$ is still only $`10`$% of $`a_{41}/4`$, less than in the 0 K case. These results show that even a large softening of the elastic constant $`C^{}`$ does not change the character of the phase transition, which remains displacive second order.

Figure 13 shows also that, above the transition temperature, the minimum energy corresponds to the cubic cell and therefore the results about the high temperature structural stability of the $`c`$ phase discussed in the previous Section remain valid.

## IV Conclusions

The $`ct`$ phase transition of pure stoichiometric zirconia has been studied from different theoretical perspectives. Both symmetry arguments and the lattice dynamical analysis suggested that this transformation might be second-order. In the 0 K perfect $`t`$ structure, the primary order parameter $`\delta `$ has a clear definition in terms of the displacement of the oxygen atoms away from their centrosymmetric position, which define a zone boundary phonon $`X_2^{}`$. But the same definition cannot be directly applied to a finite temperature atomic configuration. Instead, we defined the more general macroscopic thermodynamic variable as an ensemble average of the displacements of the O atoms projected onto the $`X_2^{}`$ zone boundary phonon coordinate. The corresponding frequency of vibration was then calculated directly from the MD simulations by applying standard statistical mechanics techniques.

The temperature evolution of both the equilibrium order parameter $`\delta `$ and of the corresponding frequency $`\nu `$ was then followed during the MD run. The results of these simulations have been interpreted with the Landau theory: the critical exponent of 0.35 fitted to the results is between the extreme values of 0.5 and 0.25 corresponding to a second-order and to a tricritical phase transition.

Approaching $`T_c`$ from the field of stability of the $`t`$ phase, the order parameter $`\delta `$ gradually decreases up to very near the transition point, as in a displacive second-order phase transformation. In a temperature window about $`T_c`$, the large fluctuations of the order parameter degrade the quality of the averages achievable from MD simulations, but we were able to observe $`\nu ^2`$ decreasing linearly with $`T`$ by some 70 %. In contrast to the prediction of the Landau theory, no increase of that frequency in the field of stability of the $`c`$ phase was observed. At high temperatures the dynamics of the oxygen sublattice revealed a high degree of mobility of the anions and a low correlation between their motion along the $`X_2^{}`$ vibrational mode.

In order to clarify these observations, we calculated the free energy surfaces relating the two structures at different temperatures, combining constrained MD simulations and the thermodynamic integration technique. These calculations showed that the high temperature stability of the $`c`$ structure is due to the entropic contribution $`T\mathrm{\Delta }S`$ and not to a variation of the internal energy profile. Moreover, we showed that the energy surfaces of the $`c`$ phase are highly anharmonic, not only at the transition temperature, but also well above $`T_c`$. This confirms and may explain in terms of thermodynamic quantities the absence of the optical modes of vibration in the experimental spectra of Liu et al. and the postulated uncorrelated “fluidlike” motion of the oxygens about their centrosymmetric position. Similarly, the experimentally observed increase of the ionic conductivity and the “structural phenomenon” mentioned by Aldebert and Traverse  may possibly be connected with the soft dynamics of the oxygen sublattice caused by the flat energy surfaces.

Our analysis revealed the peculiar character of the $`ct`$ transformation in zirconia. Approaching $`T_c`$ from below, the variation of the free energy surfaces seems to support the thesis of a $`tc`$ displacive second-order phase transition. On the contrary, approaching $`T_c`$ from above, it is not possible to follow the softening of a particular phonon mode. The entropic term $`T\mathrm{\Delta }S`$ eliminates the double well in the free energy, and frees the atoms to move along that mode without an energy cost and without a well defined vibrational frequency. This is more akin to the high temperature phase in an order-disorder phase transformation.

###### Acknowledgements.

SF is grateful for support from the European Science Foundation, Forbairt and the British Council, and for discussions with Nigel Marks. ATP and MWF are grateful to the EPSRC for funding under Grants No. L66908 and No. L08380. This work has been supported by the European Communities HCM Network “Electronic Structure Calculations of Materials Properties and Processes for Industry and Basic Science” under grant No. ERBFMRXCT980178.


text:
# CONVECTIVE MOTION IN A VIBRATED GRANULAR LAYER

## Abstract

Experimental results are presented for a vertically shaken granular layer. In the range of accelerations explored, the layer develops a convective motion in the form of one or more rolls. The velocity of the grains near the wall has been measured. It grows linearly with the acceleration, then the growing rate slows down. A rescaling with the amplitude of the wall velocity and the height of the granular layer makes all data collapse in a single curve. This can provide insights on the mechanism driving the motion.
PACS.:
45.70.Qj Pattern formation
45.70.-n Granular systems

Granular materials display different behaviors depending on the circumstances: they can be assimilated to solids, liquids or gases . Yet the features and the origin of their behavior are far from understood and are the subject of current research. The focus of this work is on the motion that arises in a layer of noncohesive granular material when submitted to a vertical vibration. Such layers are known to develop large scale collective motions . Once a certain threshold in the acceleration is reached (close to $`𝐠`$, the acceleration of gravity) grains begin to move orderly. This movement takes a form in many respects similar to that obtained in convecting fluids heated from below (the Rayleigh-Bénard instability ). Other instabilities, such as parametric wave patterns and the formation of a heap, have also been reported for granular layers . Here we are dealing only with the closed flow of granular material, at velocities much smaller than the maximum speed of the vibrating container.

This convective flow is usually downward near the wall, and the grains rise at the center of the container . For a given frequency, if the acceleration is increased there is a bifurcation towards a pattern of rolls: the velocity oscillates spatially along one preferred direction . This symmetry-breaking bifurcation is known to display strong hysteresis. Further augment of the acceleration leads to a disordered state.

Some mechanisms have been proposed to explain this phenomenon. The most obvious is the friction of the grains with the walls, as suggested by Knight et al. . Other effect concerns the Reynolds dilatancy notion: a tightly packed granular material must increase its volume if it has to change shape. During each cycle, the material undergoes a process involving a compaction and a dilatation that has been proposed to be at the origin of the motion . And yet another analysis —based on the assumption that a “granular temperature” drives the flow, as in thermal convection— has been recently proposed . The role of the surrounding gas has also been described in the case of the formation of a heap : it acts as a lubricant; but it is generally admitted that it is not at the origin of the motion. Notice that all these explanations are not mutually exclusive. While they may all contribute to the convection, it is difficult to ascertain their relative significance in a particular situation.

A measurement of the velocity can clearly help to understand the underlying physics. But carrying out such a measurement in a non invasive way is complicated. The use of trace particles is possible but painful and time-consuming. Magnetic resonance imaging has been employed to obtain velocity profiles along a cross section of the container , but the setup is entangled by the fact that an electromagnetic shaker cannot be placed inside the machine. We have arranged an experiment that allows us to obtain the local velocities at the wall of the container. While the nature of the flow is three dimensional, the velocity at the boundaries provides a good way to define it.

The experimental setup consists of a cylindrical box with the lateral wall and top made of polycarbonate, and the base made of aluminum. The diameter of the box is 52 mm. Care was taken to avoid electrostatic effects by sprinkling the box with an antistatic spray and changing the beads whenever they are seen to stick to the wall. The box is attached to the shaft of a shaker. The electromagnetic shaker (TiraVib 52110) is able to deliver a sinusoidal vibration to the box with a distortion smaller than 1%, and its stiffness in the transversal direction is such that the residual acceleration is less than 0.05 $`𝐠`$ typically. The shaker is driven by an amplifier that is in turn commanded by a function generator, from where the frequency and amplitude are set. Two accelerometers are attached to the box to measure acceleration in the vertical and one horizontal direction. The acceleration is acquired by a digital oscilloscope. Both the function generator and the oscilloscope are linked to a PC computer. An image acquisition system is also employed. It consists of a standard CCD camera with a macro lens. Images are recorded on a commercial VCR for subsequent treatment.

The procedure to measure the velocities involves the acquisition of the bead traces under proper illumination. A spatiotemporal diagram is obtained by registering a vertical line of pixels at regular time intervals and stacking them. The angle of the traces yields the bead velocities. In order to measure the average angle we perform a Fourier transform of the image and find the position of the peaks. All errors combined yield to an uncertainty of about 10% at most.

We have observed that the speed of the beads near the walls depends on height. They move faster in the upper zone and their velocity decreases as they approach the bottom, where the grains slow down and are entrained into the bulk —the velocity has a radial component to close the flow. Therefore we have taken our measurements in the upper third of the granular layer. We checked that the bead velocity there does not depend on height, their radial velocity being negligible. In that zone it is possible to track a bead and to obtain a measurement of the downward velocity. Nevertheless, unavoidably the beads suffer small displacements in the azimuthal direction due to the rearrangements of the layer, and so the traces have a small dispersion. The source of error coming from the dispersion of the velocities adds to the fact that near the bottom the beads leap downward when another bead beneath is absorbed into the bulk.

Several different granular materials have been used. The results we are reporting here have been obtained with glass beads of $`0.5\pm 0.1`$ mm in diameter. Nevertheless, the same phenomena are also observed in sand, for example, with grains approximately that large but much more variable in size and form. Several layer depths were explored in the range $`20<H<70`$, where $`H`$ is the layer depth normalized by the bead diameter. In addition, the measurements shown here have been carried out at a frequency $`f`$=110 Hz, although in a wide range (from about 70 Hz to 120 Hz) the behavior was qualitatively similar. The non-dimensional acceleration $`\mathrm{\Gamma }=A\omega ^2/\text{g}`$ (where $`A`$ is the amplitude and $`\omega =2\pi f`$ the angular frequency) for which measurements were made ranges from about $`0.5`$ to $`7`$. For $`\mathrm{\Gamma }>7`$ the movement becomes rather disordered to make meaningful measurements.

Let us show the sequence of events for $`H=20`$, in the understanding that a similar picture emerges for other layer depths. When $`\mathrm{\Gamma }`$ is small, the beads begin to move uniformly down the wall. Usually a small heap (much shorter than the height of the layer) is formed. If the box is well leveled (we do it to within $`\pm `$ 0.5 mrad) the heap appears at the center of the layer. As $`\mathrm{\Gamma }`$ is increased, the heap decreases in height and by $`\mathrm{\Gamma }3`$ the layer surface is almost flat; it just curves down near the wall. At the same time, the velocity grows larger and larger. At this stage, the value of the velocity near the wall does not depend on the azimuthal position along the wall.

A further increase of acceleration produces a bifurcation towards another pattern: two rolls appear, breaking the circular symmetry of the previous planform (see Fig. 1). The acceleration threshold shows little —if any— dependence on the frequency in the range explored; it is about $`\mathrm{\Gamma }=6`$. This transition, as remarked by Aoki and coworkers , presents hysteresis: once the rolls have formed, if acceleration is decreased they revert to the previous pattern at a lower value of $`\mathrm{\Gamma }`$. The orientation of the rolls in the circular box from one run to another is seemingly random.

A profile of the velocities along the azimuthal direction is shown in Fig. 2 before and after the bifurcation. The spatial modulation of the velocity field is clearly noticeable. The projection of the velocities onto the vertical cross section of the cylinder is well fitted by a sine. The resemblance of these two patterns with those found in small aspect ratio Bénard convection is striking. When a fluid layer open to the atmosphere heated from below begins to move, fluid goes up at the center and descends near the wall. If the temperature at the bottom is increased, a pair of rolls appear forming a pattern quite similar to the one found in granular materials . The spatial patterns obtained bear a remarkable similarity. This strongly suggest that the convective motion in a shaken granular layer may be governed, at least in principle, by the symmetry restrictions imposed by the container. In order to provide the complete boundary conditions, the velocity field at the bottom of the container is needed. Work is in progress to obtain it following the same method.

In the range of frequencies and heights explored, we have always observed this transition from a circular pattern to two rolls. Further increase of the acceleration (explored only for some values of $`H`$) leads to an increase in the number of rolls, as described in previous experiments . The rolls are superimposed to the downward flow at the wall (see Fig. 2 for $`\mathrm{\Gamma }=6.10`$; note that the mean velocity is not zero).

In order to investigate the dependence of the velocity $`𝐯`$ near the wall with the acceleration, we will focus in the first pattern, the one that appears for small $`\mathrm{\Gamma }`$. This planform is concordant with the cross-section profiles obtained by Knight et al. . The velocity field takes the form of a torus, the grains rising at the center and going down near the walls. As it is not dependent of the azimuthal coordinate, a single value of the velocity determines the field. The dependence of $`𝐯`$ on the acceleration is shown in Fig. 3a for two layer depths: $`H=20`$ and $`H=40`$. The velocity is larger for thicker layers, and it grows monotonously with $`\mathrm{\Gamma }`$. The relationship between $`𝐯`$ and $`\mathrm{\Gamma }`$ near the threshold hints at an imperfect transcritical bifurcation (the growth of $`𝐯`$ is almost linear, as seen in Fig. 3b). The numerical results of Ramírez et al., who propose a “thermal convection” mechanism, predict a supercritical bifurcation towards their first convective state. Nevertheless, the first convective state they found is banned in the experiment for symmetry reasons, and besides their prediction is made for mean values of the velocity.

It is noticeable that before the transition from one to several rolls (for $`\mathrm{\Gamma }<6`$) two different regimes exist with two distinct slopes (Fig. 3 b). The changing of slope is more evident for the series corresponding to $`H=20`$. This behavior is more clearly displayed if the velocity $`𝐯`$ is normalized –dividing it by the amplitude of the velocity of the container wall $`𝐯_𝐰`$. In doing this we are assuming that somehow the movement of the wall is driving the flow. The normalized velocity $`𝐯/v_w`$ is shown in Fig. 4 as a function of $`\mathrm{\Gamma }`$. Clearly all series now give the same behavior: a first stage where $`𝐯/v_w`$ grows with $`\mathrm{\Gamma }`$ and a second one where it remains almost constant. This saturation strongly suggests an effect related to friction, a force that depends on the velocity. The friction could either be internal (among the grains themselves) or external (between grains and wall). It is difficult to make an educated guess to decide which is the case; additional work is needed to further clarify this point. The mechanisms invoked by Knight et al.  and by Ramírez et al.  are both compatible in principle with such a fact.

An additional rescaling can be done dividing the normalized velocity by $`H`$. After doing this, all the curves collapse in a single one (Fig. 5). This dependency with $`H`$ is revealing, but again it is difficult to decide whether it is the friction with walls or the “granular temperature” the driving force. The first one should yield a scaling with $`H`$ for a fixed cross section, as it is the case. The second one could also give a scaling with $`H`$, although the reason why it should scale *linearly* is —to our knowledge— unexplained.

In summary, we have presented the results of an experiment in which the velocity of the grains near the walls has been measured locally. This allows a quantitative description of the convecting patterns. Moreover, the behavior of the velocity shows two different behaviors as the acceleration $`\mathrm{\Gamma }`$ is increased; this strongly suggests that a mechanism associated with friction is involved. The scaling of the velocity with $`𝐯_𝐰`$ and $`H`$ favors the case of a driving force associated with wall friction, although “thermal convection” cannot be excluded.

This work has been partially funded by Project PIUNA from the Universidad de Navarra, Project PB98-0208 from DGYCIT (Spanish Ministerio de Educación) and Project HPRN-CT-2000-00158 from the European Community. We thank S. Boccaletti, W. González, H. Mancini and D. Valladares for their useful comments. JLI acknowledges AECI for financial support, and IZ thanks Asociación de Amigos de la Universidad de Navarra for a grant.

Figure Captions

* Figure 1 The surface of the convecting granular material under vertical vibration ($`\mathrm{\Gamma }=6.10`$, $`f=110Hz`$, $`H=20`$). The lateral illumination reveals the prominences of the relief, showing a planform of two rolls that protrude slightly at their axes. The surface of the pattern consisting of a single toroidal roll is almost featureless, except for a small heap at small accelerations.
* Figure 2 Velocity profiles near the wall as a function of the azimuthal angle $`\theta `$: circles, $`\mathrm{\Gamma }=2.01`$; triangles, $`\mathrm{\Gamma }=3.07`$; squares, $`\mathrm{\Gamma }=6.10`$. (This last series is noisy but it has the virtue that by chance the maximum and minimum of $`𝐯`$ are seen). Note that only a portion of the complete circumference is accessible (about 160). All the velocities are downward. It can be seen that $`𝐯`$ is constant for small $`\mathrm{\Gamma }`$ –the convective pattern is a single toroidal roll– while for larger $`\mathrm{\Gamma }`$ an spatial oscillation is revealed –two rolls are formed, as shown in Fig. 1. The measurements were taken at $`f=110Hz`$ and $`H=20`$.
* Figure 3 The dependence of $`𝐯`$ on $`\mathrm{\Gamma }`$ hints to an imperfect transcritical bifurcation (a), where a linear growth is clearly apparent for small acceleration (b). For a larger range of accelerations (but still corresponding only to the first pattern, i.e. the toroidal convective roll, for which $`𝐯`$ does not depend on the azimuthal coordinate) two linear regimes are seen. In the second one, the slope is smaller. Squares correspond to $`H=40`$ and triangles to $`H=20`$. The data shown correspond to several runs, some of them increasing $`\mathrm{\Gamma }`$ and others decreasing it: no hysteresis is found in $`𝐯`$. Measurements were made at $`f=110Hz`$. The fits are just guides for the eye.
* Figure 4 The normalized velocity $`𝐯/𝐯_𝐰`$, where $`𝐯_𝐰`$ is the maximum velocity of the wall, as a function of $`\mathrm{\Gamma }`$. Squares: $`H=20`$; circles: $`H=80/3`$; down triangles: $`H=100/3`$; up triangles: $`H=30`$; diamonds: $`H=140/3`$. The fits (logistic functions) are guides for the eye.
* Figure 5 The rescaled velocity $`\frac{𝐯}{𝐯_𝐰𝐇}`$ as a function of $`\mathrm{\Gamma }`$. Data are the same than in Fig. 4, all combined. The logistic fit is a guide for the eye.