text:
# Ab initio study of the 𝛽-tin→Imma→sh phase transitions in silicon and germanium

## I Introduction

Although the structural sequence of the high-pressure phase transitions in Si and Ge has been considered to be well known since two decades, a new high-pressure phase was found experimentally between the body-centered tetragonal structure with two basis atoms (which will be referred to the following as $`\beta `$-tin or BCT phase) and the simple-hexagonal (sh) structure in SiMcMahon1993 ; McMahon1994  and Ge.Nelmes1996 ; McMahon1996  After its space group this intermediate phase was called Imma phase and corresponds to a body-centered orthorhombic Bravais lattice with two atoms in the unit cell (which will be called in the following BCO). Thus, the experimental structural sequence is found to be cd $``$ $`\beta `$-tin $``$ Imma $``$ sh $`\mathrm{}`$, where cd indicates the cubic-diamond structure. A BCO structure was first theoretically proposed for Si by Needs and MartinNeeds1984  about twenty years ago, whereas Lewis and CohenLewis1994  predicted exactly the Imma phase for Ge a few years before the experimental confirmation. The first theoretical investigations for Si and Ge indicated for both, the $`\beta `$-tin$``$Imma and the Imma$``$sh case, a continuous second-order phase transition.Lewis1994 ; Lewis1993  At variance, the experiments show first-order phase transitions with a discontinuity of the volume in both substances and a hysteresis effect in Si.McMahon1994  For Ge the first-order transition was detected only for the $`\beta `$-tin$``$Imma transition,Nelmes1996  because the Imma$``$sh transition was expected to be beyond the experimentally accessible pressure range. The transition pressures in Ge are expected to be higher than the corresponding ones in Si, because of the strong repulsive character of the atomic potential in Ge due to the presence of $`d`$ electrons in the core.Yin1982  From the theoretical point of view, recent total-energy calculations show two first-order phase transitions $`\beta `$-tin$``$Imma and Imma$``$sh both for SiChristensen1999  and Ge.Ribeiro2000  A second-order phase transition was found for the $`\beta `$-tin$``$Imma case from elastic instabilities.Hebbache2001

First information on the order of the phase transition can be given by group-theoretical arguments:Boccara1968  If the order of the point group of one phase is one half of the order of the point group of the other phase, then the phase transition can be of second order. If the point group of one phase is one third of the order of the point group of the other phase, then the phase transition must be of first order. The order of the point group $`m\overline{3}m`$ of the cd structure (space group $`O_h^7`$, $`Fd\overline{3}m`$) is 48. Since the order of the point group $`4/mmm`$ of the $`\beta `$-tin phase (space group $`D_{4h}^{19}`$, $`I4_1/amd`$ ) is 16, the cd$``$ $`\beta `$-tin transition has to be of first order. The order of the point group $`mmm`$ of the Imma phase (space group $`D_{2h}^{28}`$, $`Imma`$) is 8, thus the $`\beta `$-tin$``$Imma transition can be a second-order one. Furthermore, the order of the point group $`6/mmm`$ of the sh structure (space group $`D_{6h}^1`$, $`P6/mmm`$) is 24, the sh$``$Imma and also the Imma$``$sh transitions are of first order. Following the group-theoretical analysis none of the available total-energy calculations have described correctly the order for the phase transitions of the whole sequence, neither for Si nor for Ge.

This paper is organised as follows: In Section II we give a short overview of the methods on which our calculations are based. In Section III we report the results for the equilibrium parameters of these phases and compare them with the experimental volume dependence. Transition pressure ranges are presented in Section IV. The order of the phase transitions is analysed in Section V, and the enthalpy barriers between the phases are determined in Section VI. All results are summed up and discussed in Section VII. The conclusions are summarised in Section VIII.

## II Method

Our calculations have been performed using the plane-wave pseudopotential approach to the density-functional theory (DFT)Hohenberg1964 ; Kohn1965  implemented within the Vienna ab-initio simulation package (VASP).Kresse1996  For the exchange-correlation potential we have used the generalised-gradient approximation (GGA) by Perdew and WangPerdew1992  for Si and Ge, and, because of poorer results, the local-density approximation (LDA)Perdew1981 ; Ceperley1980  for Si only. The interaction with the ion cores has been described by ultrasoft pseudopotentials.Vanderbilt1985 ; Kresse1994  The wave functions have been expanded in terms of plane waves up to a kinetic-energy cutoff of 270 eV (410 eV) for Si (Ge). This choice provided an error smaller than 0.5 kbar (0.2 kbar) for Si (Ge) to the pressure according to the Pulay stress.Pulay1980  The Pulay stress arises from using an incomplete plane-wave basis set and causes errors in the stress tensor and, correspondingly, in the pressure. Convergence also required a 18$`\times `$18$`\times `$18 (24$`\times `$24$`\times `$24) mesh of Monkhorst-Pack pointsMonkhorst1976  which yielded 864 (1962) k-points in the irreducible wedge of the Brillouin zone for Si (Ge). All high-pressure phases are metallic. Thus, we have used a Methfessel-Paxton smearingMethfessel1989  with a width of 0.2 eV. In order to determine the total energy for a given volume or structure we have minimised the energy with a conjugate-gradient algorithm.Press1986 ; Teter1989 ; Bylander1990  The equilibrium properties have been obtained by fitting the values of the energy as a function of the volume alternatively to the VinetVinet1986  and to the MurnaghanMurnaghan1944  equation of state.

## III Structural relaxation and equilibrium properties

In this section we present our results for the volume dependence of the structural parameters and compare them with the experimental ones. The structural parameters are more useful indicators of phase transitions than the often-used tangent construction to the equation-of-state curves, the difficulties of which will be illuminated.

We have focussed our attention to the structural sequence $`\beta `$-tin$``$Imma$``$sh. The calculations for the cd phase have been performed just for the sake of comparison. The most general structure which is compatible with all four phases is the body-centered orthorhombic cell with a basis of two atoms at (0,0,0) and (0,$`0.5b`$,$`\mathrm{\Delta }c`$) (BCO structure). The other structures are special cases: For the BCT structure ($`\beta `$-tin) one has the lattice constants $`a=b`$ and the internal parameter $`\mathrm{\Delta }=0.25`$; for the sh structure one has $`b=\sqrt{3}c`$ and $`\mathrm{\Delta }=0.5`$; for the cd phase one has $`b=a`$, $`c=\sqrt{2}a`$, and $`\mathrm{\Delta }=0.25`$. The pure Imma phase is realised as a BCO structure with $`abca`$ and $`0.25<\mathrm{\Delta }<0.5`$. The minimization of the total energy has been performed using the following procedure: For all phases we have started from the BCO structure and minimised the energy as a function of the structural parameters with the constraints appropriate to the various phases. In this way we have avoided an energy offset which arises from using different cells. For example, in the case of the cd structure we have imposed, according to the previous considerations, $`b=a`$, $`c=\sqrt{2}a`$, and $`\mathrm{\Delta }=0.25`$, and then we have minimised the total energy as a function of $`a`$. In the special case of the sh structure, we could keep only the value of $`\mathrm{\Delta }`$ fixed at 0.5, and we have relaxed the value of $`b/c`$ (and also $`c/a`$) by minimization of the energy; the structure which is obtained after relaxation will be indicated in the following as SH. Indeed, converged calculations should be yield $`b/c=\sqrt{3}`$ (SH $`=`$ sh) after relaxation. For the Imma phase, all lattice and internal parameters of the BCO structure were allowed to relax.

The relaxed parameters as a function of the volume are displayed in Figs. 1 and 2 for Si and Ge, respectively. The results from calculations performed withing the GGA and the LDA are quite similar for Si. Since it turned out, that the pressures calculated by GGA are more accurate in Si, only GGA calculations have been performed for Ge. As can be seen in Fig. 1, the relaxation of $`b/c`$ for the SH phase of Si does not match the ideal value. At variance, the calculation for Ge does reproduce this value. This indicates that for a full convergence a higher kinetic-energy cutoff and/or a larger grid of special points should have been used for Si. However, the condition $`b/c=\sqrt{3}`$ is almost fulfilled in the region of stability of the sh phase. Furthermore, the energy obtained from a calculation in the same cell with $`b=\sqrt{3}c`$ fixed and $`c/a`$ relaxed differs from the value obtained using the previous procedure by less than 0.7 meV. For a more detailed discussion see Ref. Gaal2003, .

The agreement between our data and those experimental ones which are not included in the figure is good: The reported value of $`c/a=0.552`$ for the $`\beta `$-tin phase of SiHu1984  is well within the range of our values. The experimental values of $`b/a=0.917`$ and $`c/a=0.530`$ at the volume $`V=24.43`$ Å3 of Ref. Duclos1990,  agree with our data. Furthermore, the measured value of 24.22 Å3 for the volume of the sh phase of SiHanfland1999  is also in the stability range we have found for this phase. For Ge in the BCT structure, a $`c/a`$ ratio of 0.548–0.554 was detected experimentallyJamieson1963 ; Baublitz1982 ; Quadri1983 ; Menoni1986 ; Vohra1986  corresponding to a volume region around 33 Å3; in this region, our $`c/a`$ ratio is 0.546–0.549. For the sh phase of Ge, a $`b/a`$ ratio of 0.93 was measuredVohra1986  which also agrees with our data.

Figures 1 and 2 show that the relaxation of the general BCO structure leads to the BCT structure for large volumes and to the SH structure for small ones. In the intermediate region, extending from about 26.5 to 28.0 Å3 for Si and 23.5 from 25.5 Å3 for Ge, we localise the stability range of the pure Imma phase. Whereas the transition between the Imma and the sh phase is characterised by a discontinuity in most of the structural parameters, the Imma and the $`\beta `$-tin phase are quite difficult to distinguish by an inspection of Figs. 1 and 2. A clearer picture is obtained from the electronic structure and the Fermi energy (see Appendix A for details). From this procedure, the stability range of the Imma phase is from 26.6 to 28.3 Å3 for Si (GGA), from 26.6 to 28.0 Å3 for Si (LDA), and from 23.7 to 25.4 Å3 for Ge.

The standard procedure for calculating the volume change and the pressure at the transition is the so-called common-tangent construction. It relies on the fact that, in the case of a first-order phase transition, two structures with the same pressure $`p`$ and enthalpy $`H`$ but different total energy $`E(V)`$ coexist. The volume derivative of the total energy is the negative pressure. Therefore, the slope of the common tangent to the two energy curves gives the pressure for the transition between the corresponding structures. Due to the fact that calculations are performed at a finite number of volumes, a fit of an analytical expression to a continuous energy curve is needed in order to obtain a continuous derivative.

In Fig. 3 we show the energy of the relaxed structures for Si as a function of volume. The volume region where the pure Imma phase is stable does not include any minimum of the total energy $`E`$. This makes a fit to any equation of state for the Imma phase slightly unstable, which results in an unreliable determination of the transition pressures and of the volume change through the common-tangent construction in this region. The results from a fit of the VinetVinet1986  and MurnaghanMurnaghan1944  equation of state are shown in Table 1 for Ge as well as a comparison of the LDA and GGA calculations for Si.

Typically, the GGA tends to overestimate the equilibrium volume whereas the LDA tends to underestimate it. For the cd phase of Si the experimental equilibrium volume at room temperature isMcMahon1993 ; Bergamin1999 ; Basile1994 ; Holloway1991 ; Windisch1990  40.05 Å3, which lies between the LDA and GGA ones, with a deviation of about 2$`\%`$. The experimental value for cd Ge isQuadri1983 ; Holloway1991 ; Yoshiasa1997 ; Cooper1962  between 45.00 and 45.31 Å3, with a deviation from our result of less than 6$`\%`$. Because the other phases exist only under pressure, experimental values do not exist to which our equilibrium parameters could be compared.

If the standard common-tangent construction is used in order to extract the transition parameters, an accurate determination of both, equilibrium volumes and curvature of the energy vs. volume curve, is needed. This applies even more so in the present case in which $`E(V)`$ curves and equilibrium volumes of the $`\beta `$-tin, Imma and sh phases are very close together. Even though the present calculation are numerically involved the resulting data of Table 1 are not precise enough. For example, in all calculations (LDA and GGA) for Si the fitted (Vinet and Murnaghan) energy curves of the Imma phase fall below the ones of the $`\beta `$-tin phase, even for large volumes. In these cases the Imma phase is more stable than the $`\beta `$-tin phase. Thus, within this construction either a phase transition from the cd to the $`\beta `$-tin phase can be observed, if one neglects the Imma phase, or a phase transition from the cd to the Imma phase can be observed without a locally stable $`\beta `$-tin phase. A direct transition from the $`\beta `$-tin to the Imma phase cannot be detected. Similar drawbacks are found for Ge. Furthermore, the method is applicable only to the investigation of first-order phase transitions. For all these reasons, it is very difficult to extract the reliable information on the transitions using the tangent construction. In the following section, we consider an alternative procedure for analysing the transitions which does not contain the drawbacks of the above one.

## IV Determination of the transition pressures and hysteresis effect

In this section we present results for the pressure dependence of the structural parameters in order to give an estimate of the transition pressures by considering the stability range of the Imma phase.

The pressure $`p`$ in a crystal can be directly calculated from the stress tensor $`𝝈`$. Within the BCO structure the off-diagonal components of the stress tensor vanish. Therefore, the most general form of the stress tensor of our structures is

$`𝝈=\left(\begin{array}{ccc}p_x& & \\ & p_y& \\ & & p_z\end{array}\right).`$ (4)
The results of the previous section correspond to hydrostatic pressure

$`p=p_x=p_y=p_z,`$ (5)
and, in our results of Section III, Eq. (5) is fulfilled within an error of 0.1 kbar for Si and 0.05 kbar for Ge, which is less than the pressure error due to the Pulay stress. For each of the relaxed BCO structures as described in Section III, we have calculated the hydrostatic pressure corresponding to the structural parameters as shown in Figs. 1 and 2. The pressure dependence of the structural parameters is presented in Figs. 4 and 5. There is a clear discontinuity of the parameters of the Imma phase and the sh phase. The values of the parameters of the pure Imma phase are highlighted as special values of the relaxed BCO structure. For both Si and Ge, the $`\beta `$-tin$``$Imma phase transition seems to be continuous. For Si, the relaxed parameters calculated within LDA and GGA are very similar, but the corresponding pressures differ as shown in Fig. 4. In contrast, the volume range of stability in both cases is similar (see Fig. 1). Usually, the GGA transition pressures are closer to the experimental ones than the LDA values. Therefore, we will focus in the following mainly on the GGA results.

The assignment of the phases from the data points in Figs. 4 and 5 are taken from Section III. In the case of the GGA calculation for Si we have found the following sequence with pressure increase: For the $`\beta `$-tin$``$Imma phase transition, the transition pressure is approximately between 104 (last relaxation of the system to a BCT structure) and 109 kbar (first relaxation of the system to a BCO structure with $`\mathrm{\Delta }0.25`$). The next transition is the one from the Imma to the sh phase. As can be seen in Figs. 4 and 5, there is an overlap of different structures for certain values of the pressure, which can be interpreted as a hysteresis effect. That means that the transition value for increasing and decreasing pressure is different. This effect was first observed first in the experiment for Si.McMahon1994  We estimate the hysteresis effect for Si as follows: The last appearance of the pure Imma phase ($`\mathrm{\Delta }0.5`$) without any SH structure at the same pressure is at 186 kbar, the first appearance of the sh phase is at 191 kbar (downstroke), the last appearance of the pure Imma phase is at 197 kbar and the first appearance of the sh structure without any parallel Imma phase is at 198 kbar (upstroke). The transition pressures estimated in this way are listed in Table 2.

Since we have calculated the lowest energy structure, the values in Table 2 give an estimate of the hysteresis effect. An alternative approach to the hysteresis effect will be described in Section VI.

## V Order of the phase transitions

The thermodynamical potential for a pressure-induced phase transition is the Gibbs free energy. Neglecting temperature effects, the Gibbs free energy is equal to the enthalpy $`H=E+pV`$. For calculating the enthalpy, we have used the energy $`E(V)`$ of the relaxed structures with a volume $`V`$ and the hydrostatic pressure $`p`$ from Eq. (5).

At a given pressure, the stable structure is the one with the lowest enthalpy. Therefore, the transition pressure $`p_t`$ for a transition $`A`$$``$$`B`$ is defined by $`H^A(p_t)=H^B(p_t)`$. For a first-order phase transition the enthalpy curves $`H^A(p)`$ and $`H^B(p)`$ of the two phases cross at $`p_t`$, and their derivatives $`\mathrm{d}H/\mathrm{d}p`$ at $`p_t`$ are different for the two phases. For a second-order transition, the second derivatives of the enthalpy curves at $`p_t`$ for the two phases are different, whereas the first derivatives are equal. In this case the corresponding enthalpy curves have a boundary point at $`p_t`$.

In the following, we analyse the transitions in order to extract the corresponding transition pressure and order. The enthalpy $`H(p)`$ and its first derivative are presented in Fig. 6 for Ge. The enthalpy curves have been reduced by the value of the BCT structure for the $`\beta `$-tin$``$Imma transition,

$`\mathrm{\Delta }H(p)=H(p)H^{\mathrm{BCT}}(p)`$ (6)
and by the value of the SH structure for the Imma$``$sh transition,

$`\mathrm{\Delta }H(p)=H(p)H^{\mathrm{SH}}(p).`$ (7)
The derivative indicated by black dots in the lower panels of Fig. 6 has been calculated by finite differences, the straight solid line is a linear fit of these data, obtained by using only points belonging to the numerically stable regime (790–1030 kbar). An integration of this linear function is fully consistent with the enthalpy values in the upper panels. The enthalpy differences of the points in the pressure range $`p<790`$ kbar stem from the numerical noise of our calculation.

For the Imma$``$sh transition (right panels of Fig. 6) the two enthalpy curves cross at the transition pressure $`p_t=1009`$ kbar. Therefore, this phase transition is clearly of first order. The change of the derivatives at the transition point is the volume change $`\mathrm{\Delta }V(p_t)`$ and can be directly read from the lower panel of Fig. 6.

For $`\beta `$-tin$``$Imma (left panels of Fig. 6) the enthalpy curves are very close to each other in the region where the phase transition is expected. From this point of view, transitions of both first and second order are conceivable. We have thus analysed our data assuming alternatively a first- and a second-order transition.

Assuming a second-order transition the transition pressure is obtained by looking at the crossing of the $`\mathrm{d}H/\mathrm{d}p`$ curves.

A precise value of the $`\beta `$-tin$``$Imma transition pressure assuming first order can not be determined by a crossing of the enthalpy curves due to the numerical accuracy of our calculation. We are only able to give a rough estimate of the pressure at which the transition can occur due to an extrapolation of the more reliable data points away from the transition pressure. The results for Si have been obtained in an analogous way.

For both Si and Ge, the values of $`p_t`$ and of the volume changes between the low- and the high-pressure phase are shown in Table 3.

For Si the volume changes $`\mathrm{\Delta }V`$ for an assumed first-order $`\beta `$-tin$``$Imma transition from LDA and from GGA and are very small compared to the changes for Imma$``$sh. Thus we conclude that the phase transition $`\beta `$-tin$``$Imma in Si most probably is of second order.

For Ge the volume changes for a first-order $`\beta `$-tin$``$Imma and Imma$``$sh transition are within the same range. The value of $`p_t=792`$ kbar is higher than the ones from the estimated range (750–763 kbar) of the previous section, and overestimating $`p_t`$ results in overestimating $`\mathrm{\Delta }V`$. A reason for the large value of $`p_t`$ can be found in the fact, that the last reliable data points are too far away from the region where the phase transition is expected to occur. Thus, it is not possible to decide whether the transition $`\beta `$-tin$``$Imma is of first or of second order. Because of the similarities of Si and Ge a second-order transition seems more probable. If this phase transition is a discontinuous one, it is at least weakly of first order.

Using the enthalpy method the transition pressure for cd$``$$`\beta `$-tin in Si is calculated to be larger than the one for the $`\beta `$-tin$``$Imma transition. These two phase transitions take place in a small pressure range and therefore numerical errors in the determination of the transition pressure may lead to a wrong sequence. See the discussion at the end of Section III.

## VI Energy surfaces and enthalpy barriers

In this section, we present the results for the energy surfaces for Ge (similar calculations have been done for Si) as well as for the enthalpy barriers for the cd$``$$`\beta `$-tin and Imma$``$sh transitions.

We have obtained the energy surface by calculating the energy corresponding to given values of the parameters $`V`$, $`b/a`$, and $`c/a`$ and a relaxed internal parameter $`\mathrm{\Delta }`$. On the energy surface the pressure is non-hydrostatic except one or two lines or a few points. Next, we consider the pressure defined in Eq. (5) in the non-hydrostatic case

$`p_0={\displaystyle \frac{1}{3}}\left(p_x+p_y+p_z\right)`$ (8)
where $`p_x`$, $`p_y`$, and $`p_z`$ are the components of the diagonal stress tensor, Eq.(4).

For every first-order phase transition, there exists an energy barrier which has to be overcome on the path from the one phase to the other. For most reactions, such as, e.g., adsorption processes, it is a barrier of the total energy $`E`$ which lies between two local minima in the corresponding energy space. For this kind of transition, the reaction path and the energy barrier between the two phases can be found by using, e.g., the nudged-elastic-band method,Mills1995 ; Jonsson1998  by which the lowest energy path between two local minima is detected. However, this method requires the existence of well-defined local energy minima corresponding to the two phases.

If the transition, for instance the one from the cd to the $`\beta `$-tin phase in Ge, would occur without any influence of pressure, the corresponding energy space would be the energy surface $`E(V,c/a)`$ drawn in Fig. 7. The reaction coordinate here is taken as $`c/a`$. As expected, there are two local minima and a saddle point between them. By symmetry, the condition $`p_x=p_y`$ is valid all over the energy surface. The energy minima and the saddle point lie on the contour line $`p_x=p_z`$, which indicates the hydrostatic condition. The energy barrier is then defined as the energy difference between the saddle point and the starting minimum.

However, in our case the phase transition is driven by the pressure, and therefore the associated thermodynamical potential is the Gibbs free energy or the enthalpy, if temperature effects are neglected. The enthalpy surface corresponding to the cd and $`\beta `$-tin phases for Ge is shown in Fig. 8. The enthalpy for each set of parameters ($`V,c/a`$) is calculated as $`H(V,c/a)=E(V,c/a)+Vp_0(V,c/a)`$. In the range displayed in Fig. 8 the enthalpy surface has no minimum. Therefore, the usual path methods are not applicable here. However, the enthalpy barrier can be calculated in a different way as follows: The phase transition between the $`A`$ and $`B`$ phases occurs at the pressure $`p_t`$ at which $`H^A(p_t)=H^B(p_t)`$. Hence, we have to look for the points along the isobars in Fig. 8 which have the same value of the enthalpy.

In order to visualise the enthalpy barrier, we draw in Fig. 9 the enthalpy at constant average pressure (with $`p_x=p_y`$) as a function of the reaction coordinate $`c/a`$. For convenience the enthalpy is reduced by the starting point with $`c/a=\sqrt{2}`$,

$`\mathrm{\Delta }H(c/a)=H(c/a)H(\sqrt{2}),`$ (9)
which determines here the cd structure. As can be seen in Fig. 9, there exists just one enthalpy line at constant pressure, which touches the zero line twice according to $`H^A(p)=H^B(p)`$. The corresponding pressure is the transition pressure $`p_t`$ and the height of the maximum between these minima is the height of the enthalpy barrier for this phase transition. Note that on the isobars at most three points refer to hydrostatic pressure, while the other points corresponds just to the condition $`p_x=p_y`$.

Similar pictures can be drawn for the cd$``$$`\beta `$-tin transition of silicon. The height of enthalpy barrier is found to be 423 meV for Ge within GGA and 515 and 508 meV for Si within GGA and LDA, respectively. In agreement with previous results,Gaal2002  the barrier for Ge is smaller than for Si.

For the cd$``$$`\beta `$-tin phase transition the determination of the enthalpy barrier is easy because only the volume and the $`c/a`$ ratio can vary, and the latter can be used as the reaction coordinate. For the $`\beta `$-tin$``$Imma$``$sh transition all the parameters $`c/a`$, $`b/a`$, $`b/c`$, and $`\mathrm{\Delta }`$ are changing (but only two of the ratios are independent). We have calculated the total energy for each set of parameter ($`V,c/a,b/a`$), the internal parameter $`\mathrm{\Delta }`$ has been relaxed for each set. Selected energy contour plots at constant volume are presented in Fig. 10. In these plots, the places where two of the diagonal components of the stress tensor are equal are explicitly shown by the dotted, dashed, and dash-dotted lines. Obviously, hydrostatic condition is indicated by the points where all of these lines cross. Different hydrostatic pressures at one volume (like at $`V=24`$ Å3) can lead to a multivalued pressure-volume relation.

In order to calculate the enthalpy barrier in analogy to the cd$``$$`\beta `$-tin transition, we present a contour plot of the enthalpy as a function of the volume $`V`$ and of the $`b/a`$ ratio in Fig. 11. Actually, the enthalpy depends also on the other variables $`c/a`$ and $`\mathrm{\Delta }`$. The contour plot of Fig. 11 has been obtained using the following procedure: For each value of the set ($`V,b/a`$), we have chosen $`c/a`$ so as to fulfill the condition $`p_z=p_x`$ (see Fig. 10), which is suggested because at the end we will look for the points which fulfill the hydrostatic condition. The value of $`\mathrm{\Delta }`$ has been taken so as to minimise the energy for a given choice of the other parameters. The two-dimensional presentation of the results as in Fig. 11 corresponds to the choice of $`b/a`$ as the reaction coordinate. Completely equivalent results can be obtained by choosing $`c/a`$ or $`\mathrm{\Delta }`$ instead. The hydrostatic condition in Fig. 11 is then realised by $`p_y=p_z`$.

The comparison with our previous results shows that most of the data shown in Fig. 5 lie on the hydrostatic line. Numerical instabilities lead to deviations (generally of less than 1 kbar) between the relaxed points and the line $`p_y=p_z`$. Nevertheless, the relaxed points show the same behaviour (to a lesser extent) as the hydrostatic line, which is shown schematically in the inset of Fig. 11.

The schematical behaviour of the hydrostatic curve between approximately 23.5 and 24 Å3 as shown in the inset of Fig. 11 can be used to discuss the hysteresis effect qualitatively. If one follows the hydrostatic line $`p_y=p_z`$ with increasing pressure $`p`$ (decreasing volume) starting from A, one arrives at the point C where the direction of the line is changing. From C to D, the pressure decreases with a simultanous decrease of the volume which is an instable situation. Thus, at C the phase transition must occur through a jump from C to E, which is at the same pressure. From E increasing pressure leads to F. Following the path in the opposite direction starting from F, i.e., by decreasing the pressure (volume increase), the direction of the curve changes at D, which results in a jump to B. As a consequence of that it is possible to have a higher transition pressure with increasing than with decreasing pressure. Therefore, an alternative estimate for the limit of the hysteresis effect can be given. For Si within GGA (LDA), we obtaine a transition pressure of 198 (133) kbar with increasing and of 169 (109) kbar with decreasing pressure. For Ge within GGA, the corresponding values are 1010 and 965 kbar.

A behaviour similar to the one depicted in the inset of Fig. 11 can be found also for the cd$``$$`\beta `$-tin phase transition (see Fig. 8). There, the low-pressure edge corresponding to the point C of the inset of Fig. 11 is at a negative pressure, which would lead for $`\beta `$-tin$``$cd to a very low or even negative transition pressure, the latter case being connected to the irreversibility of the phase transition. For the cd$``$$`\beta `$-tin transition the volume and the enthalpy as a function of hydrostatic pressure is presented in Fig. 12. The values have been extracted from Fig. 8 along the hydrostatic line $`p_x=p_z(=p_y)`$. The ideal cd structure ($`c/a=\sqrt{2}`$) has been reached within an error of 1% and therefore cd is noted as CD in analogy to the difference between sh and SH (see Section III). In order to discriminate the enthalpy curves of Fig. 12 we have subtracted a linear background to arrive at the curves $`\mathrm{\Delta }H`$. The solid and the dashed line mark the CD and the BCT structure, respectively. The dotted line indicate the values along the part of the hydrostatic line, which connects the CD and the BCT structure. The designated points correspond to the inset of Fig. 11.

The condition of local stability is fulfilled along all curves shown in Fig. 12. Thus, the volume vs. pressure curves are strictly monotonically decreasing and the enthalpy vs. pressure curves are convex. The condition of global stability requires the enthalpy to be minimal, and so the curves from A to T and from T to F mark the globally stable regime. The point of the phase transition in analogy to the Maxwell construction is noted as T. The pressure range of the coexistence regime is determined by the overlap of the globally stable and the instable regions, which is here from B (or D) to E (or C).

A similar picture could be drawn for the Imma$``$sh transition, but with less resolution. A reliable identification of the structures within the $`H(p)`$ curve analoguos to Fig. 12 requires a finer grid of data points for the interpolation and/or a higher convergence.

The enthalpy barrier for the cd$``$$`\beta `$-tin phase transition can be determined with a procedure like illustrated in Fig. 9. In Figs. 13 and 14, we show the enthalpy at given pressures calculated within GGA as a function of the reaction coordinate $`b/a`$, for Si and Ge, respectively.

The enthalpy in Figs. 13 and 14 is reduced by the value at $`b/a=1`$ (corresponding to the $`\beta `$-tin phase),

$`\mathrm{\Delta }H(b/a)=H(b/a)H(1).`$ (10)
The dotted line in Fig. 13 corresponds to a pressure where the enthalpy has a local minimum for the $`\beta `$-tin phase. At a pressure of about 130 kbar we find a continuous transition to the Imma phase. Above this pressure, a minimum appears at $`b/a<1`$. By further increasing $`p_0`$ a second minimum apprears at $`b/a0.92`$, corresponding to the sh phase. At the pressure of 194 kbar the two minima have the same enthalpy, therefore a transition can occur. The enthalpy barrier between the two phases is surprisingly small, namely 1.2 meV. At higher pressure the sh phase is stable. For the LDA calculation a transition pressure of 133 kbar and an enthalpy barrier of 1.3 meV are obtained. If one considers the enthalpy curves of Fig. 13 at the transition pressures determined in the Sections IV and V, one can give a further estimate of the enthalpy barriers, and, the values obtained following this procedure do not differ appreciably from the ones mentioned above.

As can be seen in Fig. 14 for Ge, a minimum of the enthalpy can be identified corresponding to a pure Imma phase, slightly below 820 kbar. However, within our resolution (about 0.1 meV) the pressure curves have only one minimum. Therefore, an enthalpy barrier between the Imma and the sh phase can not be identified. Indeed, this transition is of first order. Hence, the existing barrier must be smaller than 0.1 meV.

All the enthalpy barriers which we are able to obtain for the Imma$``$sh transition are much smaller than the thermal energy at room temperature (25 meV). Thus, experiments performed at room temperature should not be able to distinguish between a first- or a second-order phase transition.

## VII Discussion of the results for the phase transitions

The transition pressures for the cd$``$$`\beta `$-tin$``$Imma$``$sh phase transitions have been calculated with different methods in the previous sections. In Sections IV and VI an estimate of the hysteresis range was given by considering the relaxed values (R) and the hydrostatic line (HL) respectively. In Section V the transition pressure was determined via the enthalpy vs. pressure curves (EC). In Table 4 an overview of all our results for the transition pressure is given in comparison with the experimental results.

For Si the value of the pressure for the cd$``$$`\beta `$-tin transition calculated with GGA is within the range of the experimental values whereas the LDA calculation underestimates it. We could not reproduce the experimentally observed hysteresisMcMahon1994  for the $`\beta `$-tin$``$Imma transition in Si. For a first-order phase transitions, such as Imma$``$sh, a hysteresis could occur due to the behaviour of the hydrostatic line as in Fig. 11. Since our results support a second-order phase transition for $`\beta `$-tin$``$Imma rather than a first-order one, no hysteresis should be found in this case. The transition pressures mentioned in Ref. McMahon1994,  are obtained either from the first appearance of a new phase or from the complete transition of the sample from one phase to the other. Thus we estimate the experimental range in which the $`\beta `$-tin$``$Imma and Imma$``$sh transitions occur in analogy to the procedure explained in Section IV. Unfortunately, the calculated transition pressures for $`\beta `$-tin$``$Imma lie below the ones for the first cd$``$$`\beta `$-tin phase transition and underestimate the experimental ones.

The next phase transition to be considered is Imma$``$sh. Here a hysteresis was found within our calculations as well as in the experiment. The results for the transition pressures within LDA underestimate again the experimental ones, whereas the ones within GGA overestimate them. Not included in Table 4 are the experimental pressures for mixed (50:50) Imma-sh samples, which are 154 kbar for pressure increase and 159 kbar for pressure decreaseMcMahon1994 . This is the only case in which the transition to the (mixed) Imma phase occur at a lower pressure for an increasing than for a decreasing external pressure. All other experimental points show the same kind of hysteresis which is reproduced in our results. The difference of the pressures, where the Imma phase was measured with pressure increase and decrease is 5 kbar for the mixed samples and 11–12 kbar for the pure ones. From the experimental transition pressures the hysteresis is found to be between 3 and 14 kbar. Within our results the largest possible hysteresis is found to be around 29 (24) kbar and the smallest around 1–12 (2–7) kbar for GGA (LDA), which is in agreement with the experimental values. The transition pressures determined by the crossing of the enthalpy curves are found to be within the hysteresis range calculated here.

For Ge the transition pressure for the cd$``$$`\beta `$-tin transition is slightly smaller than the experimental one but the results for the $`\beta `$-tin$``$Imma transition match much better. In the experimentNelmes1996 ; McMahon1996  just the $`\beta `$-tin$``$Imma transition was examined, because the transition pressure for the Imma$``$sh transition was higher than the maximum accessible pressure of 810 kbar. In comparison to the experimental results of Si it is assumed that the Imma phase for Ge is stable from 750 to 850 kbar, which is compatible with the observation of a sh phase at 900 kbar.Vohra1986  From this point of view, the transition pressure calculated here is slightly too high.

There are many reasons for the discrepancy between our results and the experimental ones. On one hand, the experimental conditions like pressure profiles or measurement times are not perfectly controlled, whereas nonhydrostaic conditions affect the transition pressures; on the other hand, temperatue effects, which are neglected in our calculation, could be of some relevance. In fact, for Si the theoretical transition pressure decreases with increasing temperature.Gaal2001

## VIII Conclusions

We have presented a first-principles investigation of the pressure-induced phase transitions cd $``$ $`\beta `$-tin $``$ Imma $``$ sh in Si and Ge. Numerical accuracy has required these calculations to be performed within the same unit cell and the pressure to be calculated from the stress tensor instead of using any equation of state. We have determined the phase transitions by investigating the relaxed parameters, the enthalpy vs. pressure curves, and the hydrostatic condition. From the pressure dependence of the relaxed structure parameters the hysteresis effect does not seem to have been verified in a theoretical investigation before. An examination of the hydrostatic condition provides a thermodynamical explanation for the existence of this hysteresis. Further on, the orders of the phase transitions have been investigated with the helps of the curves of the enthalpy vs. pressure and of the derivative of the enthalpy vs. pressure. For Si the transitions Imma$``$sh and $`\beta `$-tin$``$Imma are found to be of first and of second order, respectively. For Ge the results support a second-order $`\beta `$-tin$``$Imma transition more than a first-order one, but in this case the order cannot be definitely determined. A calculation of the enthalpy barriers between the phases show a vanishing barrier between the $`\beta `$-tin and the Imma phase for Si, which also supports the second-order $`\beta `$-tin$``$Imma transition, whereas no barrier could be found for all $`\beta `$-tin$``$Imma$``$sh transitions. Finally, our alternative results are fully consistent with each other and show an acceptable agreement with the available experimental data.

## Acknowledgement

Support by the Heinrich Böll Stiftung, Germany, is gratefully acknowledged.

## Appendix A Fermi energy and band structure

We examined the volume dependence of the Fermi energy of BCT, BCO, and SH. In Fig. 15 we show the Fermi energy reduced by the BCO values for Ge,

$`\mathrm{\Delta }E_\mathrm{F}(V)=E_\mathrm{F}(V)E_\mathrm{F}^{\mathrm{BCO}}(V).`$ (11)
In the range of stability of the $`\beta `$-tin phase the Fermi energy of BCO is identical with the one of BCT. Similarly, the Fermi energy of BCO and SH are identical in the range of stability of the sh phase. Therefore, the regime where the Fermi energy of BCO is neither equal to the one of BCT nor to the one of SH is assumed here to be the range of stability of the pure Imma phase.

We analysed the volume dependence of the band structure for the BCT, BCO, and SH structure, too. The results for the lowest four bands at the $`\mathrm{\Gamma }`$-point for Ge are displayed in Fig. 16.

The energy values $`E_i`$ are reduced in such a way that the Fermi energy at any volume is set equal to zero,

$`\mathrm{\Delta }E(V)=E_i(V)E_\mathrm{F}(V).`$ (12)
In most of the volume range presented in Fig. 16 the band structure presents a linear dependence on volume. However, the energy gradients differ for each structure. Within the transition region between the Imma and the sh phase one observes a small discontinuity (here apparent at $`V=23.64`$Å3) which is attributed to the first-order phase transition. When starting from the $`\beta `$-tin phase, a deviation of the BCO band energies from the BCT energies can be observed, which increases continuously. Thereby, the degeneracy of the BCT bands is lifted, which indicates the appearance of the less symmetrical Imma phase. The behavior of the band structure for Si is very similar to the one of Ge.


format: text

text:
# Direct Detection of a CME-Associated Shock in LASCO White Light Images

## 1 Introduction

White light coronagraphs regularly observe fast, impulsive coronal mass ejections (CMEs) with speeds that exceed both the coronal sound speed ($`200`$ km/s) and the Alfvén speed ($`800`$ km/s at 4 R Mann et al., 1999). Therefore, CMEs are capable of driving wave perturbations and even shocks in the low corona (e.g., Hundhausen, 1988). However, the imaging of CME-associated shocks remains an observational challenge. There are two observational approaches to this problem.

The first approach relies on the study of shock proxies. The existence of shocks during mass ejections is supported by a large amount of indirect evidence, such as radio type-II observations (e.g., Cliver, Webb & Howard, 1999) and distant streamer deflections (Gosling et al., 1974; Michels et al., 1984; Sheeley, Hakala & Wang, 2000). The contribution of type-II observations to the the study of coronal CME shocks is limited by the lack of imaging of type-II sources, the frequent concurrence of CMEs and flares (flares can also drive shocks in the low corona), and the uncertainty in the CME initiation times. Although there are some promising recent results (Maia et al., 2000), the connection between type-IIs and CMEs will likely remain unclear for the near future. Observations of distant streamer deflections offer the best, indirect, evidence for white light shocks. Such observations have been more numerous in recent years thanks to the increased sensitivity and temporal coverage of the LASCO coronagraphs (Sheeley, Hakala & Wang, 2000).

The second approach to the problem is to establish shock signatures directly from the white light CME images. Initially, it was thought that the looplike front of many CMEs was the density enchancement from a fast MHD shock (e.g., Maxwell and Dryer, 1981; Steinolfson, 1985). But Sime, MacQueen & Hundhausen (1984) pointed out several discrepancies between the model predictions and the observations. They argued, for example, that many looplike CMEs propagate too slowly to form a fast shock. This led Hundhausen, Holzer & Low (1987) to suggest slow shocks as candidates for some CME fronts. Steinolfson and Hundhausen (1990a, b) simulated in some detail the appearance of both slow and intermediate shocks on coronagraph images. However, the applicability of the simulations to the CME analysis is hindered by the complexity of the ejected structures. Many CMEs have irregular fronts or no well-defined fronts at all and it is usually difficult to differentiate between coronal material, which is inherently looplike, and shock-related structures.

In fact, there has been only one published identification of a white light shock despite the observations of thousands of CMEs since the early 1970s. The observation of an unusual looplike CME with the Solar Maximum Mission (SMM) coronagraph (MacQueen et al., 1980) led Sime and Hundhausen (1987) to suggest that the loop front could be a fast shock wave. They argued that the high lateral expansion speed (800 km/s) of the CME and the absence of any deflections of coronal structures before the loop impinged on them were strong indicators of a shock. There were some problems with this interpretation, however, that arose from the limitations of the available observations. For example, the restricted field of view and sensitivity of the SMM images could not rule out the existence of a disturbance ahead of the observed CME nor could provide a reliable determination of the acceleration profile of the CME. One could only note that the loop did not appear to deccelerate during the streamer crossings contrary to theoretical (Odstrc̆il and Karlický, 2000) and observational (MacQueen and Fisher, 1983; Sheeley, Hakala & Wang, 2000) expectations. Finally, the nature of the loop front could not be established without observations in other wavelengths.

The high quality LASCO observations provide an excellent dataset for locating possible white light signatures. We searched the LASCO database for CME events which were simple enough to allow an unambiguous identification of a white light shock. The best case is an event on April 2, 1999, when a fast ejection with an exceptionally clear density enhancement along its flanks was observed. The availability of CME simulation codes and high quality datasets, in several wavelengths, provide us with the necessary tools to investigate the nature of the observed enchancement more thoroughly than it was possible in the past. We start with an overview of the available observations in § 2. The results from the analysis of the observations and the MHD modeling are presented in § 3 and we discuss their implications in § 4. We summarize our findings in § 5.

## 2 Observations

On April 2, 1999, the Large and Spectroscopic Coronagraph (LASCO, Brueckner et al., 1995) and the Extreme Ultraviolet Imaging Telescope (EIT, Delaboudiniére et al., 1995) aboard the Solar and Heliospheric Observatory (SOHO, Domingo, Fleck & Poland, 1995) observed a jet-like CME (jet-CME hereafter) at the northeast solar limb (Figure 1). Ejections of this type are characterized by their small widths ($`10^{}20^{}`$), limited latitudinal expansion and simple structure (Gilbert et al., 2001; Dobrzycka et al., 2003). We can easily identify the source region in the EIT images because part of the ejecta can be seen in absorption in the 195Å image (Figure 1). The jet-CME originated in NOAA active region 8507. An M1.1 soft X-ray flare was reported from the same location (Solar Geophysical Data Reports). The flare began at 8:03 UT, and peaked at 8:21 UT. The EIT images between 8:12 – 8:54 UT showed a series of dark spray ejections (Figures 1-2). By 8:30 UT, the front of the jet-CME appeared in the C2 field of view at a height of 3.5 R (Figure 2). The ejection maintained its narrow width ($`20^{}`$) out to 30 R. The filamentary structure of the ejecta clearly points to a filament eruption in association with the jet-CME. A filament was visible in the Hα image from the Meudon observatory taken at 8:17 UT. This is not a typical filament eruption because the event lacks the familiar loop-cavity-core configuration. This is apparent upon inspection of the C2 image, which illustrates a loop-cavity-core shape CME occurring over the western limb (Fig. 1). The western CME has a well-defined loop front followed by a cavity and what appears to be a core just over the C2 occulter. By comparison, the jet-CME event appears to consist of the main body of the CME without a loop or cavity preceeding the ejecta. We can think of two possible explanations for the lack of overlying streamer material in the jet-CME:

* The jet-CME occurred outside a streamer which is not an uncommon situation. In a study of a large sample of LASCO CMEs, Subramanian et al. (1999) found that 27% of the events were displaced from preexisting streamers.
* The overlying corona was disrupted by an earlier CME (hereafter CME1) that was first seen in C2 at 1:30 UT and had just exited the C2 field of view when the jet-CME followed in its wake. The apparent position angle of CME1 was $`91^{}`$, close to the position angle of the jet-CME ($`72^{}`$). CME1 was also wide enough ($`74^{}`$) to extend over the location of the jet-CME. We carefully inspected the LASCO-C2 images to look for signatures of coronal depletions after CME1. The only significant brightness depressions were localized over the location of CME1. However, this result does not preclude the possibility that most of the material overlying the jet-CME was located below $`2.3`$ R, the inner cutoff of the C2 coronagraph. In that case, no depletions would be visible in the LASCO images. Such a coronal configuration is also consistent with the lack of any obvious relation of AR8507 to a C2 streamer during its east limb passage. Therefore, the observations do not dismiss the possibility that CME1 carried away part of the streamer material over the site of the jet-CME, leading to the unusual appearance of the jet-CME.

The most intriguing aspect of this event, however, is the sharpness of the southern flank of the CME. A much fainter counterpart is barely visible along the northern flank. The CME bears an uncanny resemblance to a fast projectile and its associated bow shock. This similarity prompted our investigation on whether the sharp CME flank is actually the white light counterpart of a shock.

We searched the available datasets in other wavelengths for possible shock evidence in association with this CME. No metric type-II emission, the main proxy for coronal shocks, was reported by any radio spectrographs. Our own analysis of the Potsdam spectrograms revealed only a group of type-III emissions between 8:09-8:15 UT, which was most likely connected to the flaring in the active region. The Ultraviolet Coronagraph Spectrometer (UVCS, Kohl et al., 1995) was observing along the northeast limb between 8:11 and 9:44 UT. Unfortunately, the UVCS slit did not intercept the CME flanks but it intercepted part of the filament in the CME core. The spectra suggest that the core is an untwisting filament (A. Chiaravella, personal communication). The Transition Region and Coronal Explorer (TRACE, Handy et al., 1999) 171Å and the Yohkoh/Soft Xray Telescope (Tsuneta et al., 1991) images, during this event, do not show any coronal features or ejecta that could be identified with the sharp white light flank.

We conclude that these observations do not provide strong evidence either for a shock or for a coronal structure that could be associated with the white light feature. To establish the nature of the white light feature, therefore, we rely on the analysis of the LASCO data in conjuction with an MHD simulation.

## 3 Analysis of the White Light CME

### 3.1 LASCO Measurements

We analysed the kinematics of the event by constructing height-time plots along two radial positions using the EIT, C2 and C3 images. The first position angle (PA), marked as P1 in Figure 2, corresponds to the front of the filament and consequently of the CME (PA=80). The second position (P2) was taken at a random position along the southern CME flank (PA=100). The position angles were measured counterclockwise from the solar North Pole. Figure 3 demonstrates that the height-time curves can be fitted well by second degree polynomials. Both curves show deceleration. The average speeds in both locations (800-1000 km/s) are well above the median speed of 476 km/s for CMEs observed with LASCO (St. Cyr et al., 2000). They are comparable to shock speeds deduced from metric type-II bursts (e.g., Klassen et al., 2000) and are similar to the measured speeds for the July 6, 1980 white light shock candidate (Sime and Hundhausen, 1987). The derived speeds, therefore, are consistent with the existence of a shock. Note that, as all coronagraphic measurements, the above heights, speeds and decelerations are projected values in the plane of the sky. In this case, however, we are confident that the measured parameters are close to their true, radial values because the CME source region is very close to the solar limb and since most CMEs appear to move radially (Hundhausen, Burkepile & St. Cyr, 1994; St. Cyr et al., 1999, 2000) the jet-CME was most likely propagating close to the plane of the sky.

We measured the mass across the sharp feature and derived the corresponding density profile, as follows: First, the LASCO C2 and C3 images were corrected for vignetting, exposure time and other instrumental effects and were calibrated in units of excess brightness after subtracting a pre-event image. This is a standard procedure for the analysis of coronagraph CME observations (Poland et al., 1981) and, for the LASCO case, is described in some detail in Vourlidas et al. (2000). The resulting excess brightness images were converted to excess mass images under the usual assumptions: (i) the scattering electrons are concentrated on the plane of the sky and, (ii) the ejected material comprises a mix of completely ionized hydrogen and 10% helium. These images are shown in Figure 2. The visibility of the CME flank is enhanced somewhat by this procedure and it can be followed out to at least 20 R. We measured the mass profiles across the front at the same position angle (P2) where the height-time measurements were taken. To improve the signal-to-noise ratio, especially at the larger elongations, we averaged over a 10-wide swath along the flank. The result is a measurement of the line of sight density of the white light structure. To convert the profiles to volume density we made another assumption about the unknown depth of the feature along the line of sight. The sharpness of the feature suggests a limited extent along the line of sight. Therefore, we assumed its depth equal to its width ($`0.2`$ R), as measured in the C2 image at 8.30 UT. The resulting profiles represent estimates of excess density. We converted them to total density profiles by adding the pre-event (background) coronal density. This density was derived from the closest C2 polarization brightness (pB) image using the well-known pB inversion method (van de Hulst, 1950; Hayes, Vourlidas & Howard, 2001). The pB image was obtained at 21:00 UT on April 1, before the passage of the earlier CME discussed in the previous section. It is likely that the derived background density might be slightly higher than the actual density at the time of our event. The final density profiles are shown in Figure 4 where they are compared to the commonly-used SPM background model (Saito, Poland & Munro, 1977) coronal density profile. We see that the background density profile, derived from the C2 pB image, agrees very well with the SPM model profile. The density increase across the feature is very sharp in the C2 images and the profiles are very similar to the expected density profile across a shock. The profile retains its sharpness until 9:18 UT, becoming smoother at larger distances (not shown here). The density jump across the profile is about a factor 3 (at 8:30 UT) and strongly suggests that the sharp white light feature is actually a shock. However, one should note that, given the assumptions used to derive the density profiles, this density jump is only an estimate. To examine the viability of a shock at the flanks of the CME, we simulate the event using the measured speeds and background density profiles as constraints.

### 3.2 MHD Simulation

We start with a time-dependent, 2D plane-of-sky ideal MHD model. The model is described by the standard equations of MHD theory (e.g., Priest, 1982) with additional momentum and heating terms to accommodate the bimodal (fast, slow) solar wind (Wang et al., 1998; Wu et al., 2000). The complete set of equations is given in Wang et al. (1998). The model incorporates four constraints that are based on the LASCO observations. Namely, (i) the CME does not have a loop or a cavity preceeding its core, (ii) the CME center is located at 80 counterclockwise from the solar north pole, (iii) the event shows significant deceleration and (iv) the simulated density profiles should match the observed ones (Figure 4). Finally, we are not interested in the details of the initiation of the event and thus we use a simple initiation mechanism which is described later.

To accommodate the first constraint we need to consider the initial magnetic field topology. Simulations with closed field configurations and beta ratios of the order of unity tend to lead to loop-cavity-core CME morphologies with considerable latitudinal extent while low beta, open field configurations result in elongated, laterally confined ejections (see Figure 5 and Table 1 in Steinolfson et al., 1978). Therefore, closed field configurations are better suited for modeling CMEs similar to the west limb event in Figure 1 or to the loop-like CME described in Sime and Hundhausen (1987). Instead, the narrow width of the CME and the absence of the 3-part structure, led us to the choice of a radial magnetic field (Figure 5) with beta ratios ranging from 0.83 at the equator to 0.34 at the poles. Figure 5 also shows the initial solar wind velocity and density distribution at the position angle of the CME core (P1). The model density is in good agreement with the measured density. The background flow velocity is obtained from the MHD model. Note that the contrast between the CME center and polar background flow velocities is insignificant because we use a radial magnetic field configuration. To account for the observed deceleration, we used a spatially dependent polytropic index ($`\gamma =1.05`$ between 1-4 R, then linearly increasing to 1.45 at 30 R). This implies that the energy equation for the ideal MHD model is modified to include non-adiabatic processes (Wang et al., 1998). To reduce computation time, we assumed symmetrical conditions at the angular computational boundaries at 0 and 80 from solar north. In other words, the CME is considered symmetric relative to its center at P1 which is consistent with the appearance of the CME in the LASCO images. The lower radial boundary is the solar surface (1 R) and is prescribed according to the theory of characteristics (Wu and Wang, 1987). The initial conditions at the solar surface at locations P1 and P2 are: magnetic field 1.5 G and 1.2 G, number density $`5\times 10^7`$ cm-3 and $`10^8`$ cm-3, and temperature $`2.24\times 10^6`$ K and $`1.74\times 10^6`$ K, respectively. A detailed mathematical representation of the compatibility relations obtained by methods of characteristics in two dimensions is given in Suess, Wang, & Wu (1996), for example. In short, we specify that $`B_r`$ and $`p/\rho ^\gamma `$ are constants at the lower boundary (the solar surface), and $`p`$, $`B_\theta `$, $`\mathrm{v}_\mathrm{r}`$ and $`\mathrm{v}_\theta `$ are computed from the compatibility relations (Wu and Wang, 1987; Wang, 1992; Suess, Wang, & Wu, 1996). The linear extrapolation method is used to compute the upper radial computational boundary at 30 R which has non-reflecting boundary conditions. The system is in quasi-steady equilibrium in its initial state. The same modeling method was applied in the simulation of another event in Wu et al. (1999).

The sole objective of our simulation is the identification of the nature of the flanks of this particular jet-CME, not the processes of CME initiation. Thus, we do not seek elaborate initiation configurations that attempt to model CMEs in more realistic ways but impose a premium on computational time. Instead, we initiate the event with a simple pressure pulse. In any case, the nature of the CME flanks should not depend on the details of its initiation mechanism. We introduce a 20-wide pressure perturbation by raising the density but keeping the temperature identical to the local value at the solar surface, which is symmetric relative to the center of the event (P1). This density pulse is three times higher than the background and the maximum perturbed velocity at its center reaches 200 km/s at 200 seconds, maintains this value for about 5000 seconds, and then declines to the background velocity. This perturbation represents a mass flow that forms as an ejecta to produce the observed shock. The total computational time is 293 minutes.

Figure 4 demonstrates the excellent agreement between the simulated and measured density profiles. The simulated and measured height-time curves are compared in Figure 6. Overall the simulation reproduces the CME measurements well. We are also interested in reproducing the morphology of the event. In Figure 7, we assemble the observed images together with simulated pB images and the magnetic field configuration. The core of the CME is immediately recognized acting as the piston for the bow-shock feature. We see that the choice of radial field led our model to a good match of the narrow width of the CME.

#### 3.2.1 Analysis of Simulation Results

The next step is the analysis of the waves produced in the simulation. In an MHD medium, the information from a disturbance can propagate with one of two characteristic speeds; the slow-mode, $`V_s`$ or the fast-mode, $`V_f`$ speed that are given by

$$V_{f,s}^2=0.5\left[(V_\alpha ^2+C_s^2)\pm \sqrt{(V_\alpha ^2+C_s^2)^24V_\alpha ^2C_s^2\mathrm{cos}^2\theta }\right]$$
(1)
where $`V_\alpha `$ is the Alfvén speed ($`V_\alpha =B/\sqrt{4\pi \rho }`$), $`C_s`$ is the sound speed ($`C_s=(\gamma p/\rho )^{1/2}`$), and $`\theta `$ is the angle between the magnetic field and the direction of the wave propagation. For propagation parallel to the magnetic field ($`\theta =0`$), we find that $`V_f=V_\alpha `$ and $`V_s=C_s`$. In the case of propagation perpendicular to the magnetic field, only the fast-mode propagates with speed $`V_f=(V_\alpha ^2+C_s^2)^{1/2}`$.

When the wave is compressive, it may be observed in the coronagraph images if the density increase at the wave front is large enough. For sufficiently energetic drivers, the wave may steepen into a shock. Note that all MHD shocks are compressive shocks. From shock theory (Jeffrey and Taniuti, 1964), the existence of a MHD shock solution for propagation parallel to the field requires

$`V_f>U_{sh}>V_s,`$ $`\mathrm{slow}\mathrm{MHD}\mathrm{shock}`$
$`U_{sh}>V_f,`$ $`\mathrm{fast}\mathrm{MHD}\mathrm{shock}`$ (2)
and for propagation perpendicular to the magnetic field requires

$`U_{sh}>(V_\alpha ^2+C_s^2)^{1/2}`$
where $`U_{sh}`$ is the normal component of the relative shock speed

$`U_{sh}=(\stackrel{}{V}_{prop}\stackrel{}{V}_{SW})\stackrel{}{e}`$ (3)
where $`\stackrel{}{e}`$ is the normal to the shock front, $`V_{prop}`$ is the propagation speed of the disturbance and $`V_{SW}`$ is the solar wind speed. In other words, $`V_{prop}`$ is the usual speed derived from the coronagraph measurements and $`U_{sh}`$ is the speed of the disturbance in the solar wind frame. To derive the propagation speed from the simulation we follow Wu et al. (1996) and identify the locations of the shock ($`r_1,r_2,\mathrm{},r_n`$) at times $`t_1,t_2,\mathrm{},t_n`$. Then

$`V_{prop}={\displaystyle \frac{r_{i+1}r_i}{t_{i+1}t_i}}`$ (4)
The simulation results for $`U_{sh}`$, $`V_{f,s}`$, $`V_{SW}`$ and $`V_{prop}`$ are shown in Tables 1 and 2 for both the CME flank and front, respectively. The characteristic speeds are calculated at a point in the pre-shock region. We see that the wave velocity is larger than the fast-mode speed at both locations. Therefore, our simulation supports the existence of a shock at the flanks of the CME. But could this shock be visible in the white light images? We can check the visibility of the shock by estimating the density increase across the shock front as follows. For the fast-mode shock the Mach number is defined as

$`M_\alpha ={\displaystyle \frac{U_{sh}}{V_f}}`$ (5)
The maximum density increase across the shock front occurs for propagation along the magnetic field

$`{\displaystyle \frac{\rho _2}{\rho _1}}={\displaystyle \frac{(\gamma +1)M_\alpha ^2}{(\gamma 1)M_\alpha ^2+2}}`$ (6)
where $`\rho _1`$,$`\rho _2`$ are the upstream and downstream densities, respectively. For propagation at any other direction, the density compression is smaller than Eq.(6) and there is no exact formula for it. Only for very large Mach numbers, the density compression reduces to

$`{\displaystyle \frac{\rho _2}{\rho _1}}={\displaystyle \frac{\gamma +1}{\gamma 1}}`$ (7)
for any propagation direction. From Eqs (5) and (6) we can now calculate the maximum density compression for the shock. The results are shown in Tables 1-2 for the CME flank and front, respectively. It may be noted that we use the polytropic index $`\gamma `$, and not the ratio of specific heats, in computing the density compression ratios. The results are close to the measured density compression, taking into account that the actual shock does not propagate along the magnetic field. The actual compression and Mach numbers should be less than the predictions in Tables 1-2.

In summary, the simulation leads us to the following conclusions:

* The observed bow shaped feature at the flank of the CME is the density enhancement from a fast mode MHD shock.
* The shock strength measured by its Mach number remained high throughout its propagation in the C2 and C3 fields of view.
* The observed and model density enhancements (Fig. 4) show very good agreement. The measured density compression ratios are slightly lower than those predicted by the shock jump condition. This is expected because the shock is propagating at an oblique angle to the largely radial magnetic field, as an inspection of the images suggests.
* The deceleration of the shock and the front is simulated by using a spatially dependent polytropic index. We do not model the energy dissipation processes that result in the deceleration; instead, we use a spatially varying polytropic index as a proxy for the energy dissipation. We have used the value of the spatially varying polytropic index at the location of the shock in computing the density compression ratio. In the absence of an explicit energy equation for the flow and knowledge about radiation losses at the shock, we feel that this is a reasonable approach. If we did solve an explicit energy equation for the flow, we would use the ratio of specific heats in order to compute the ratios of various quantities across the shock. However, this too would not be correct if the shock were lossy; one would then have to employ a separate energy balance equation across the shock that takes these losses into account.

## 4 Discussion

As we pointed out in the introduction, it is generally difficult to differentiate between the white light signatures of shocks (or waves) and those of coronal ejecta such as loops and filaments. For the event we examined here, we can exclude the possibility that the flanks of the CME correspond to coronal structures for several reasons.

First, there is no visible ejecta ahead of the jet-CME. The absence of overlying material might be due to an earlier, larger CME event that entered the C2 field of view at 1:30 UT. The earlier CME propagated along the same path as our event and possibly carried away some of the overlying corona. It is likely that the coronal density above the east limb was still low when the second eruption occurred.

Second, there is no evidence for preexisting or erupted large coronal loops (on scales comparable to the white light front) in the observations of every available wavelength, from EUV to X-rays.

Third, the flanks of the CME diffused rapidly as they propagated outwards. This behaviour is inconsistent with that of a coronal structure. This is readily apparent from a comparison to the filamentary structures in the CME core that remain identifiable at large heliocentric distances.

Fourth, the southeastern streamer is clearly deflected when the CME flank reaches its location. We then see the deflection propagating to the adjacent streamers thus creating a ”front” traveling southwards with a slower speed than the CME front or its flanks (Figures 1-2). This behavior is expected from a shock on theoretical grounds (Odstrc̆il and Karlický, 2000; Van der Holst, Van Driel-Gesztelyi & Poedts, 2002). It is, also, quite common in LASCO CME images where it has been interpreted as an indirect signature of shock or waves caused by high-speed CMEs (Sheeley, Hakala & Wang, 2000).

From the above arguments, we deduce that the sharp white light feature at the flank of the CME is likely a wave. Whether it is actually a shock or not depends on the local coronal conditions. However, no direct measurements of the magnetic field or the coronal density exist. Thus we employed a simulation (§3.2) to assess the likelihood of a shock. As the high CME speeds suggest and the results of the simulation support, the CME flank is likely the signature of a fast mode MHD shock. The lack of metric type-II emission does not contradict our conclusion since metric type-II bursts do not always accompany fast CMEs and vice versa (e.g., Cliver, Webb & Howard, 1999; Gopalswamy et al., 2000).

To our knowledge, this event is the best case for a direct identification of a white light shock in a coronagraph image. It doubles the number of CME shock candidates to date (Sime and Hundhausen, 1987, this work). Clearly, this is a very small sample and does not permit any conclusive arguments on whether the CME fronts or other CME features seen in LASCO images are indeed shocks. It is important to extend this work with further studies. There are many fast LASCO CMEs that exhibit sharp features, along either their fronts or their flanks (Fig. 8). Based on the present work, it is likely that these features could be shock or fast mode wave signatures. This information might be proven useful in the analysis of these events by providing, for example, the possible locations of type-II sources.

Another important result, in our opinion, is the simultaneous observation of a streamer deflection and the shock wave responsible for it. Although, streamer deflections have always been interpreted as the indirect signature of a CME-associated disturbance, the “hard” evidence was missing, until now. It is apparent, that the disturbance propagates through the streamer at a lower speed than in the open corona. This is most likely due to the high streamer density but projection effects may come into play (Sheeley, Hakala & Wang, 2000) and would probably complicate the derivation of any physical parameters about the streamer (or the wave). Although it could yield interesting results on the properties of shocks and streamers, no such analysis has been carried out so far, to our knowledge.

Another noteworthy consideration is the scale of the event. We showed that apparently small spray-like ejections in the low corona can give rise to a CME with significant influence over a much larger area. Yet, the brightest part of the CME (its massive core) is constrained within a narrow width ($`<20^{}`$). If this was an earth-directed event, it would have remained inside the coronagraph occulters and hence would not have been detected. These arguments raise the question whether there exists a class of undetected earth-directed CMEs that can have some influence on the terrestrial space weather conditions. This question cannot be answered until observations of CMEs along the sun-earth line become available during the planned STEREO mission.

## 5 Conclusions

We have analysed a unique CME which exhibited a sharp feature at its flanks reminiscent of shock signature. From the analysis of the LASCO data and other available observations, we concluded that the white light flanks could not be coronal structures. We employed a MHD simulation to assess the possibility of the flanks being the density enhancement from a shock. We used the measured speeds and coronal densities as constraints to the simulation and showed that a fast MHD shock was formed at the front and the flank of the CME. We conclude that shocks could be directly observed in white light coronagraph images, under suitable conditions.

The event also provided us with the first direct observation of a streamer deflection by a fast mode shock, thus lending firm support to the interpretation of the commonly-seen streamer deflections as proxies to CME-induced coronal shocks. Although the appearance of this CME is rather unique, it is not fundamentally different that the ejections seen daily in the LASCO images. Therefore, the aforementioned results could be used in the search for shocks signatures in other events.

We thank the referee for helpful comments that improved the manuscript. We thank H. Aurass for providing the Potsdam radioheliograph data and M. Andrews for useful comments. A.V acknowledges support from NASA grants. S.T.W and A.H.W wish to acknowledge support by NSF grant ATM-0070385 and AFOSR grant F49620-00-0-0204. SOHO is an international collaboration between NASA and ESA. LASCO was constructed by a consortium of institutions: the Naval Research Laboratory (Washington, DC, USA), the Max-Planck-Institut fur Aeronomie (Katlenburg- Lindau, Germany), the Laboratoire d’Astronomie Spatiale (Marseille, France) and the University of Birmingham (Birmingham, UK).