text:
# Critical behaviour of thermopower and conductivity at the metal-insulator transition in high-mobility Si-MOSFETs

## I Introduction

Scaling theory of non-interacting, disordered, electron gases predicts that no metal-insulator transition (MIT) occurs in 2 dimensions as temperature $`T0`$. Nevertheless, what appears to be an MIT has been observed (at finite, though low $`T`$), first in $`n`$Si-MOSFETs and more recently in many other 2-dimensional (2D) hole and electron gases. In the particular case of Si-MOSFETs, the transition is most clearly visible in high-mobility samples, roughly $`\mu 1`$ m2/V s. As the electron density, $`n`$, is varied, there is a particular value, $`n_0`$, above or below which the resistivity $`\rho `$ shows metallic or insulating temperature dependence respectively. For the present purposes we will use as a working definition that negative $`d\rho /dT`$ indicates an ‘insulator’, and positive $`d\rho /dT`$ at the lowest temperatures we can reach corresponds to a ‘metal’ (possible deviations from this definition and the consequences will be mentioned later). At $`n>n_0`$ and not too close to $`n_0`$, metallic behaviour is visible over a wide range of $`T`$, roughly $`T<0.5T_F`$ where $`T_F`$ is the Fermi temperature. The decrease of $`\rho `$ in the metallic state for high mobility Si-MOSFETs samples is typically a few orders of magnitude larger than can be accounted for by electron-phonon scattering. There is no accepted explanation of this behaviour as yet. On the insulating side $`\rho `$ is found to increase exponentially as $`T`$ is reduced, typically showing good agreement with variable range hopping models.

Most previous work on these systems has focused on $`\rho `$, though measurements of the compressibility have also appeared recently. The present paper presents experimental data on thermopower, $`S`$, and conductivity, $`\sigma =1/\rho `$, at temperatures down to 0.3 K. When we extract the diffusion thermopower, $`S^d`$, we find that it diverges when approaching $`n_0`$. At $`n_0`$ there is an abrupt change in behaviour of the thermopower with lower densities showing an upturn in $`S`$ as $`T`$ is decreased. We also attribute this to $`S^d`$.

In addition we find that $`\sigma `$ exhibits critical behaviour around $`n_0`$ which appears to be largely consistent with the behaviour of $`S^d`$. Earlier, a scaling behaviour was described for the temperature dependence of $`\rho (T)`$ over a temperature range $`(0.050.3)T_F`$. In contrast, the present experiments are concerned with the limiting, low temperature behaviour of $`\sigma `$. As already mentioned, $`\rho `$ decreases with temperature when $`n>n_0`$ but it saturates at low temperatures. (We are unable to say what happens below 0.3 K and return to this point later). Our low temperature saturation values of $`\sigma `$, say $`\sigma _0`$, show a power-law critical behaviour as a function of $`(n/n_01)`$. Similar critical behaviour of $`\sigma `$ around $`n_0`$ has previously been reported in 2D p-type GaAs and lower mobility Si-MOSFETs. The observed characteristics of $`S^d`$ and $`\sigma _0`$ are strongly reminiscent of those expected for an Anderson MIT in 3D but such a transition should not occur in 2D.

We have reported some of our preliminary results, particularly on the thermopower, in conference form. The present paper extends the thermopwer results to a wider range of temperatures and densities, and extends the data into regions of high sample resistivities which were inaccessible to us earlier. It is interesting to note that, unlike resistivity, $`S`$ is measured with no current flow in the system and hence no self heating occur. We also present the results of systematic studies on $`\rho `$ on a number of samples; these were primarily taken to augment our understanding of the thermopower, but they are also interesting in their own right.

## II Sample and Experiments

The main sample used for the present $`\rho `$ and $`S`$ measurements (Sample 1) is the same as that described in a previous paper and the general techniques used to measure $`S`$ can also be found there. This sample has $`n_0=1.01\times 10^{15}`$ m-2 (as defined as above) and a peak mobility $`\mu =1.82`$ m2/Vs at $`T=0.3`$ K. $`S`$ and $`\rho `$ have been measured as a function of $`T`$, down to about 0.3 K, at many different values of $`n`$. We have also analyzed independent $`\rho (T,n)`$ data for two other samples over the same range of $`T`$, Sample 2 from the same wafer with $`n_0=0.96\times 10^{15}`$ m-2 and with peak $`\mu =1.96`$ m2/Vs (again at $`T=0.3`$ K) and Sample 3  with $`n_0=0.96\times 10^{15}`$ m-2 and peak $`\mu =3.6`$ m2/Vs. All connections to the 2D gas had isolation resistance $`>50`$G$`\mathrm{\Omega }`$ and all leads into the cryostat were well shielded and filtered against rf interference. For the sample leads, commercial, bulkhead, rf rejection filters were used covering a range 1 MHz - 10 GHz, and these were augmented by extra, series inductors to extend the range down to about 50 kHz.

Measurements were made at dc. With the sample in the metallic state ($`n>n_0`$), a Keithley 182 digital voltmeter usually gave the best compromise of input bias current, input impedance and noise. As previously, we obtained spurious voltages from the sample even with no temperature gradient or current present. We traced this to the input bias current of the voltmeter. With any amplifier, the input bias current combined with the source impedance produces spurious (offset) voltages which must be separated from the true signals. Spurious voltages have been noticed previously by many people and probably have a number of causes, but we verified the consistency of the explanation in the present case by measuring the 2-terminal resistance of the sample over a wide range of conditions, as well as measuring the input bias current of the amplifier. The bias current of the input amplifier of our Keithley 182 was $`15`$ pA for sources up to about 1 M$`\mathrm{\Omega }`$ but it increased to $`30`$ pA at higher source impedances.

When the system was metallic the offset voltages were always rather small and not usually too $`T`$ dependent, especially at lower temperatures where the thermoelectric signals are also small, so they were not a serious problem. For example, at $`n=2.1\times 10^{15}`$m-2, the sample is metallic and has a resistance of $`20`$ k$`\mathrm{\Omega }`$ at low temperatures. This gives an offset voltage of $`30`$ nV. Offsets were taken into account by taking the thermoelectric signal to be the difference between the measured voltages with the temperature gradient present and not present, ensuring that the mean sample temperature, and thus resistance, was the same in both cases.

On the other hand, when the sample was in the insulating state its resistance was much higher and strongly $`T`$dependent. In this case the offset voltage rapidly overwhelmed the thermoelectric signal as the temperature was lowered. To reduce this problem, an amplifier with input bias current $`<1`$ pA and input impedance $`>10^{12}\mathrm{\Omega }`$ was used. With some averaging it had a resolution of 0.1 $`\mu `$V for source impedances of less than a few hundred k$`\mathrm{\Omega }`$ rising to about 1 $`\mu `$V at 10-20 M$`\mathrm{\Omega }`$. However, when used with the sample there was always an apparent offset current of $`4`$ pA independent of sample voltage, i.e., offset signals were produced corresponding to a dc current through the sample of this magnitude. We were unable to determine the origin of this current. It was unaffected by changing the rf filters in the leads, or by rerouting some of the wiring, particularly that carrying the gate voltage. When the offset signals become large, one would like to have some check that the method of extracting the signal used above was reliable. Fortunately, it was found that if the sample was grounded at different contacts, the offset current always flowed in a direction towards ground. Thus by measuring the thermoelectric signals with the source or drain grounded in turn, the offset voltage reversed sign but the thermoelectric voltage did not. This gave two independent measurements and is the origin of the error bars on our data at low temperatures with the sample in the insulating state (see Fig. (4) later). The width of the error bars typically increases very rapidly as $`T`$ is lowered because of the strong variation of $`\rho `$.

Finally we provide some information on the temperature difference $`\mathrm{\Delta }T`$ across the sample at various temperatures so that the signal voltages $`\mathrm{\Delta }V=S\mathrm{\Delta }T`$ may be estimated under various conditions. The electric field in the sample may also be estimated given the length was 2.5 mm. Between 0.7 and 7 K, $`\mathrm{\Delta }T100`$ mK. Below 0.7 K $`\mathrm{\Delta }T`$ dropped steadily until by 0.4 K it was about 30 mK and at 0.34 K about 15 mK.

## III Results and Discussion

In the metallic region $`n`$, as obtained from the Hall constant, was found to be a linear function of gate voltage $`V_g`$, and it is believed to follow approximately the same dependence in the insulating region, at least close to $`n_0`$. On different cooldowns from room temperature there were small shifts in $`n`$ at fixed $`V_g`$ of up to $`0.04\times 10^{15}`$ m-2, but if the sample was held at 77 K between cooldowns, as was done for most of these measurements, the density and other properties were very reproducible. In general the absolute uncertainty in $`n`$ should be no worse than $`\pm 0.1\times 10^{15}`$ m-2 in the region of $`n_0`$ and 2% at high $`n`$; the reason for the larger uncertainty near $`n_0`$ is the fact that the Hall resistivity shows non-linearities of up to about 10% there.

Because our analysis of thermopower makes extensive use of $`\rho `$ data, examples of the temperature dependence of $`\rho `$ on both the metallic and insulating side are shown in Figs. 1 and 2. They are very similar to those seen in previous work.

Well into the insulating regime (roughly $`n0.65\times 10^{15}`$m-2) $`\rho `$ can be well represented by an equation developed for variable range hopping (VRH) through localized states in the presence of a Coulomb gap at the Fermi energy $`E_F`$ (referred to as Efros-Shklovskii VRH),

$$\rho =\rho _c\mathrm{exp}(T_{ES}/T)^{1/2},$$
(1)
where $`\rho _c`$ and $`T_{ES}`$ are constants that depend on $`n`$ but not $`T`$. Examples are given in Fig. 1 for Sample 1. The values that we obtain for the constants are similar with those determined by others. $`T_{ES}`$ varies from $`35`$ K at $`n=0.650\times 10^{15}`$m-2 to $`90`$ K at $`0.48\times 10^{15}`$m-2. If the exponent is left as a variable we find an average value of $`0.50\pm 0.05`$ over this same range of $`n`$. We find poorer fits assuming Mott VRH (no Coulomb gap at $`E_F`$) which corresponds to a power law of 1/3 in the above expression.

For $`n0.70\times 10^{15}`$m-2, a simple activated behaviour

$$\rho =\rho _c\mathrm{exp}(T_A/T)$$
(2)
provides a very good fit to the data below about 1 K. Thus these results suggest a transition from simple activated behaviour to VRH as $`n`$ decreases.

In the metallic region ($`n>n_0`$) we have fitted our data to the equation

$$\rho =\rho _0+\rho _1\mathrm{exp}\left[(T_0/T)^q\right]$$
(3)
where $`\rho _0`$, $`\rho _1`$, $`T_0`$ and $`q`$ are constants that depend on $`n`$ but not $`T`$. Examples are shown in Fig. 2 for Sample 2. For $`n1.5\times 10^{15}`$m-2, excellent fits were obtained by fixing $`q=0.5`$. For lower $`n`$, the values of $`q`$ (if this is left as a variable) increased and approached 2.0 at $`n1.0\times 10^{15}`$m-2. The model curves fit the data successfully over the range 0.3 K (our lowest $`T`$) to about $`0.3T_F`$. ($`T_F=7.3n`$ K with $`n`$ in $`10^{15}`$ m-2 units). Above this region the curves usually exhibit a maximum. $`T_0`$ also provides a natural scale for the data and, since $`T_0`$ decreases as $`n_0`$ is approached (Fig. 2 indicates $`T_0`$ for the lowest densities), a lowering of the upper temperature limit is also consistent with this. The reduced range of fitting at low $`n`$ means that, $`\rho _1`$, $`T_0`$ and $`q`$ are not well known here but, even so, $`\rho _0`$ retains a relatively low error.

Figure 3 shows the results on $`\sigma _0=1/\rho _0`$ as a function of $`n`$. For all three samples $`\sigma _0`$ follow the critical behaviour

$$\sigma _0=\sigma _m+\sigma _s\left(\frac{n}{n_0}1\right)^\nu .$$
(4)
The solid lines are the best fits with the following parameters, with $`\sigma `$ in units of $`e^2/h`$ (the data points at the nominal values of $`n_0`$ are $`\sigma `$(0.3K) and were not used in the fitting). For sample 2, $`\sigma _m=0.0\pm 0.3`$, $`\sigma _s=13.6\pm 0.7`$ and $`\nu =0.94\pm 0.03`$. A good fit was not possible over the whole range of $`n`$ and so was restricted to $`n<3.0\times 10^{15}`$ m-2 (a range of $`(n/n_01)`$ of about 0.03 to 3). Sample 1 has an identical behaviour within experimental error. The higher-mobility Sample 3 also follows the same equation with $`\sigma _m=0.36\pm 0.15`$, $`\sigma _s=34\pm 5`$ and $`\nu =1.39\pm 0.05`$. We have a smaller range of $`n`$ for this sample, but the fit still covers a range $`(n/n_01)0.020.4`$. These results suggest $`\nu `$ increases with peak mobility but clearly more data on a variety of samples are required. The values of $`\sigma _m`$ for Sample 1 and 2 are consistent with zero within experimental uncertainty. The fits are good over the whole range except for the last points near $`n=n_0`$ which do not fall on the curves. This is not necessarily unexpected because a value of $`\sigma =0`$ is not possible at a finite temperature due to thermal broadening of the electronic distribution. Small errors in $`n_0`$ will have a similar effect. For Sample 3, $`\sigma _m`$ may be finite. However, if $`n_0`$ is allowed to decrease from $`0.96`$ to about $`0.93\times 10^{15}`$ m-2, a fit which is essentially indistinguishable over the range of the data can also be obtained with $`\sigma _m=0.0\pm 0.2`$, $`\sigma _s=32\pm 5`$ and $`\nu =1.48\pm 0.05`$; both fits are actually shown on Fig. 3. A small discrepancy in $`n_0`$ could easily arise from the identification of the critical density for the MIT with that density, $`n_0`$, where $`d\rho /dT`$ changes sign, a procedure which has no firm physical foundation; in addition there are uncertainties and possibly non-linearities in the determination of $`n`$ from gate voltage values as mentioned above.

The critical behaviour described by Eq. (4) with $`\sigma _m=0`$ is formally the same as that expected for a (continuous) Anderson transition with a ‘mobility edge’ at $`n_0`$, whereas a finite $`\sigma _m`$ would correspond to a (discontinuous) Mott-Anderson transition; neither transition should arise in a non-interacting 2D gas. Our $`\rho `$ data in the insulating region are also consistent with this scenario in the sense that we observe activated behaviour for $`n`$ below $`n_0`$. Similar critical behaviour, usually with $`\sigma _m`$ consistent with zero, has been seen in many 3D systems, typically with values of $`\nu `$ in the range $`0.51.3`$. There are also two previously reported cases related to 2D. Hanein et al. have made a similar analysis to the one above for a 2D hole gas in GaAs and found a linear relation between $`\sigma _0`$ and $`n`$, but with a finite $`\sigma _m`$. Feng et al. also found critical behaviour in $`\sigma _0`$ in a low mobility Si-MOSFET though we note that their data all show negative $`d\rho (T)/dT`$, even in the range which they identify with ‘metallic’ behaviour. They found $`\sigma (T)`$ increased with $`T`$ and fitted a $`T^2`$ dependence. This behaviour was attributed to local momemts and is so different from that found here that the two cases may not be closely related. Recently  a model based on percolation of non-interacting electrons through local quantum point contacts also led to zero-temperature conductivity consistent with our observations.

Note that the data presented in Fig. 3 were obtained from extrapolating $`\sigma `$ to $`T=0`$ by a procedure which focusses only on the ‘strong’ exponential, $`T`$dependence of Eq. (3). It ignores any other weaker dependences that may be present below $`T=0.3`$K, including those due to weak localization and screening.

We now turn to the thermopower data. A selection of data on $`S`$ as a function of $`T`$ is shown in Fig. 4. From our previous work we know that the diffusion thermopower, $`S^d`$, is almost zero at $`n=8.5\times 10^{15}`$ m-2 and one sees only phonon drag, $`S^g`$ which varies approximately as $`T^6`$ at the lowest temperatures. The fact that $`S^d`$ is very low at high $`n`$ was predicted by Karavolas and Butcher who attributed it to a particular combination of scattering mechanisms. Indeed, when $`nn_0`$ the behaviour of both $`S^d`$ and $`S^g`$ appears to be in generally good agreement with theory.

The $`T^6`$ dependence of $`S^g`$ obtains only in the Bloch limit, i.e., when $`q2k_F`$ where $`q`$ is the average phonon wave number and $`k_F`$ the Fermi wave number. At $`n=8.5\times 10^{15}`$ m-2 the condition $`q2k_F`$ is satisfied below about 1.4 K. Above this temperature one sees a gradual decrease in the exponent to roughly $`S^gT^3`$ by 4 K, but no simple power law is expected in this region.

As $`n`$ decreases, $`S`$ begins to show two distinct regions with different $`T`$ dependences. At $`T<1`$ K, $`S`$ has a much weaker, approximately linear, $`T`$-dependence indicative of $`S^d`$ becoming dominant. For $`n<n_0`$, this low-$`T`$ dependence, which is characteristic of ordinary metals, is replaced by an upturn in $`S`$; we will return to this feature later.

Concentrating first on the metallic region, we analysed the data over the range 0.3-1.5 K to obtain $`S^d`$ by assuming

$$S=S^d+S^g=\alpha T+\beta T^s$$
(5)
where $`\alpha (=S^d/T`$) and $`\beta `$ are constants that depend on $`n`$. For the 2 highest densities we obtained good straight lines by plotting $`S/T`$ as a function of $`T^{s1}`$ with $`s=6`$ as expected. At lower $`n`$, proportionately lower temperatures are required to reach the Bloch limit and so we no longer expect to see $`S^gT^6`$ in our temperature range. Hence we simply used $`s`$ as a variable to give the best fit. For $`n=1.95\times 10^{15}`$ m-2 we found $`s=4`$, and for $`n=0.971.5\times 10^{15}`$ m-2, $`s=5`$. That $`s`$ seems to increase again at the lowest $`n`$ is unexpected, but more data are required before this is taken to be an experimental fact. We have no reason to believe that $`S^g`$ is not well-behaved in this temperature, density and resistivity range, but we must caution that our understanding of $`S^g`$ in these circumstances is not well founded as yet. Given this it is important to note that the value of $`\alpha `$ depends only weakly on the choice of $`s`$. Indeed, providing $`S^g`$ has a much stronger $`T`$ dependence than linear, then simply taking the measured values of $`S/T`$ at 0.3 K to be equal to $`\alpha `$ gives values in good agreement with those obtained using Eq. (5).

Fig. 5 shows that $`\alpha `$ as a function of $`n`$ appears to diverge as $`nn_0`$. One would expect such a divergence if $`E_F`$ approaches a band edge or a gap in the DOS giving $`S^d1/n`$. Das Sarma et al. have suggested that a similar mechanism, carrier freezeout, is responsible for the MIT in 2DEGs. However, the present results are inconsistent with these explanations because Hall data, both our own and those in Ref., show that in the vicinity of $`n_0`$ the mobile carrier density is not approaching zero but equals $`n`$ within a few percent; also any temperature dependence of $`n`$ at fixed $`V_g`$ is very small in the range $`0.34.2`$ K. Given that theory successfully predicted $`S^d`$ to be zero at high $`n`$ as mentioned above, it is possible that the present results have an explanation in terms of the specific scattering processes important at low $`n`$, but there is no theory available for this situation at present.

However, Eq. (4) also implies a divergence of $`S^d`$. With the assumption of a constant density of states (DOS), Eq. (4) is consistent with

$$\sigma (E_F)=\sigma _m+\sigma _s\left(\frac{E_F}{E_c}1\right)^\nu .$$
(6)
Again, with $`\sigma _m=0`$ this is formally equivalent to an Anderson transition with $`E_c`$ being the mobility edge. Taking $`\sigma _m=0`$, the use of the Mott relation $`S^d=(\pi ^2k_B^2T/3e)(\mathrm{ln}\sigma /E)_{E_F}`$ with Eq. (6) then gives

$$S^d=\frac{\nu \pi ^2k_B^2T}{3e(E_FE_c)}.$$
(7)
The use of the Mott relation, and hence this result, is valid only if $`(E_FE_c)/k_BT1`$; in the opposite limit $`S^d`$ tends to a constant ($`228\mu `$V/K in 3D). The saturation is a direct consequence of the fact that the contribution to the conductivity of electrons below $`E_c`$ is zero in this model. The relevant integral for $`S^d`$ has a weighting factor $`\sigma (E)(EE_F)`$ which, under normal conditions, leads to $`S^d`$ being a measure of the derivative of $`\sigma (E)`$ with respect to energy, and hence to the Mott relation. This is consistent with Eq. (7). However, when $`E_F=E_c`$ only electrons above $`E_c`$ contribute to the integral, which leads to a constant value of $`S^d`$. Numerical calculations show that the approximation of Eq. (7) gives a magnitude roughly a factor of 2 too large when $`(E_FE_c)/k_BT2`$ (i.e. because of the approach to saturation). For our sample at $`T=0.4`$ K this corresponds to $`(nn_0)/n_0=\mathrm{\Delta }`$ (say) $`0.11`$, using the ideal DOS, $`g_0`$, with an effective mass of $`0.19m_0`$. If the present model was indeed appropriate, it suggests that we should see saturation of $`S^d`$ as $`T`$ increases for the samples with $`n=0.97`$ and $`1.06\times 10^{15}`$ m-2 (for higher $`n`$ any saturation would be masked by the rapidly increasing contribution from $`S^g`$). We are unable to see any obvious evidence for saturation at these densities and Eq. (5) still seems to be valid there. However, we should recall that $`S^g`$ is not yet understood at low $`n`$ and high $`\rho `$ and there remains the possibility that an unexpected behaviour of $`S^g`$ could mask important features in $`S^d`$.

Proceeding with the model, we have fitted Eq. (7) to the measured $`\alpha `$, allowing for the possibility of saturation as $`n_0`$ is approached. To simulate a saturation we add $`\mathrm{\Delta }`$ in the denominator (but allow it to be a variable when determining the best fit to the data) and, rewriting Eq. (7) in terms of $`n`$, we have

$$\alpha =S^d/T=K/\sqrt{\mathrm{\Delta }^2+(\frac{n}{n_0}1)^2}$$
(8)
where $`K=\nu \pi ^2k_B^2/(3eE_c)`$ is a constant expected to be about 32 $`\mu `$V/K2 for Sample 1, again using $`g_0`$. If $`\sigma _m`$ is finite in Eq. (6), then the Mott relation suggests that it will contribute to the denominator of Eq. (7), also softening the divergence at $`n=n_0`$. However, the experimental $`\sigma _m`$ is so small that this is probably negligible compared to the finite-$`T`$ effect considered here.

The best fit of the data to Eq. (8) is shown by the solid line in Fig. 5 and uses $`\mathrm{\Delta }=0.15\pm 0.01`$ and $`K=(9.5\pm 1.5)\mu `$V/K2. Only data with $`n<4\times 10^{15}`$ m-2 are fitted, a range consistent with that used for $`\sigma _0`$. $`\mathrm{\Delta }`$ is consistent with that expected from the argument above, but $`K`$ is too small by a factor of about 3. As with $`\sigma `$, the fit can be improved if $`n_0`$ is decreased, perhaps by $`0.05\times 10^{15}`$ m-2, in which case $`\mathrm{\Delta }`$ decreases and may even become zero.

In spite of the relatively good correspondence between the model and the data, we emphasize that we should be cautious in necessarily concluding that the model is fundamentally correct. We are comparing our results on a 2D system with a theoretical model of an Anderson MIT believed to be appropriate only to non-interacting electrons in 3D. At this point it is difficult to know whether the situation in 2D might be significantly changed by including interactions along with the strong disorder; this is a complex and ongoing theoretical problem (e.g. see Refs.  and references therein). Some progress has been made on calculating $`S^d`$ with the inclusion of weak interactions and disorder. We are, however, unable to explain the observed strong $`n`$ dependence in terms of the calculations which predict universal values. The calculated corrections are logarithmic in $`T`$. These would be the equivalent of similar corrections in $`\rho `$ and, if they exist, would require much lower temperatures than we have available and would be difficult to detect. On the other hand, these facts mean that the analysis of $`\alpha `$, $`\rho _0`$, and their relationship as zero-temperature quantities, is a self-consistent procedure which has a semi-classical meaning. Nevertheless, whether or not $`\sigma `$ goes sharply to zero near $`n_0`$ as would be required in the Anderson model, the observed change in behaviour of $`S^d`$ as $`n`$ is varied around $`n_0`$ points to a fundamental change in transport taking place in this region.

We should finally mention that we can also represent the data reasonably well over the same range using the simple expression $`\alpha =56/n^{2.5}\mu `$V/K2, with $`n`$ in units of $`10^{15}`$ m-2, but this has no obvious physical explanation; in particular, it does not have the form that we would expect for $`S^d`$ approaching a band edge at $`n=0`$, i.e. $`\alpha 1/n`$.

Our data in the insulating regime also show a behaviour qualitatively consistent with activation across an energy or mobility gap. In such cases we expect an upturn of $`S^d`$ at low temperatures according to

$$S^d=\frac{k_B}{e}\left(A+\frac{E_cE_F}{k_BT}\right)$$
(9)
where $`A`$ is a constant of order unity. In the particular case of the Anderson model $`Ak_B/e`$ corresponds to the saturation value when $`E_c=E_F`$ noted above in our discussion of the metallic region. When $`n`$ is close to, but less than, $`n_0`$, we see activated behaviour in $`\rho `$ (see Eq. (2)) which is in accord with this. We do not see a high temperature saturation of $`S`$ as Eq. (9) predicts, but this is not surprisingly because it will be masked by other contributions to $`S`$, particularly $`S^g`$ (see also below).

There might also be another contribution to $`S^d`$ from VRH. Demishev et al. have recently demonstrated the existence of such a component in $`S^d`$ using a 3D GaSb sample. The $`T`$ dependence of our $`\rho `$ data and other previously published data well into the insulating region are consistent with Efros-Shklovskii VRH. (When two or more conduction mechanisms are present, the appropriate $`S^d`$ are weighted by their contributions to $`\sigma `$). For this mechanism one expects $`S^d`$ to be a constant given by $`S^d=(k_B/e)(k_BT_{ES}/C)(\mathrm{ln}g(E)/E)_{E_F}`$ where $`T_{ES}`$ is obtained from the temperature dependence of $`\rho `$ in Eq. (1), $`g(E)`$ is the background DOS, i.e, on which the Coulomb gap is superimposed, and $`C`$ a dimensionless constant $`6`$. If we take $`(\mathrm{ln}g(E)/E)_{E_F}1/E_F`$ (implying that $`E_F`$ may be in the tail of the DOS) and again using $`g_0`$ to estimate $`E_F`$, we find that the calculated $`S^d`$ are typically a factor of two smaller than the values of $`S`$ at the observed minima. The argument is not significantly changed if Mott VRH is assumed. In this case $`S^dT^{1/3}`$ but the magnitudes calculated for $`S^d`$ are similar. The references give the relevant details.

In the insulating regime, one expects a diverging localization or correlation length as the critical point is approached. Given this, one might also question whether the thermopower results are in the linear region close to the critical point. However, on this side of the transition we are never very close to the critical density. The closest point is at about $`n_0n0.5\times 10^{14}`$m-2. The correlation length estimated in Ref.  under these conditions is about 20$`\mu `$m. Assuming a magnitude of the disorder potential of order 1 K gives a threshold electric field for non-linear effects of order 80 mV/m whereas the thermoelectric electric field is about 1 mV/m over the range of 0.7-7 K. Thus non-linear effects do not appear to be a problem for the present data.

To put the present results in perspective, as far as we are aware the only previous work which attempted to follow $`S^d`$ into the region of 2D electron localization was that of Burns and Chaikin on thin films of Pd and PdAu. They found an upturn of $`S^d`$ in the strong localization region but no divergence at higher conductivities. The authors attributed their results to the opening of a Mott-Hubbard gap. In 3D, Lauinger and Baumann observed critical behaviour of $`\sigma `$ and a divergence of $`S^d`$ for metallic AuSb films, but the magnitude of the latter was 2 orders of magnitude smaller than seen here. Other experiments on bulk SiP and NbSi saw no divergence of $`S^d`$ on the metallic side.

We close this section by making a few comments about $`S`$ at higher $`T`$. $`S^g`$ appears to be dominant in this region because the strong $`T`$ dependence of $`S`$ (roughly $`T^2`$ to $`T^3`$) is inconsistent with any other mechanism. As we have already indicated, little is known about the behaviour of $`S^g`$ where $`n`$ is low and $`\rho `$ is high, and these are just the conditions that pertain around $`n_0`$. Thus, although $`S^g`$ should be present on the metallic side, its precise form is not known with any certainty.

When conduction is dominated by VRH, we expect $`S^g=0`$ because drag is based on the conservation of crystal momentum for electron-phonon scattering which will not hold for transitions between localized electron states. Thus, $`S^g`$ on the insulator side should only exist when carrier excitation to delocalized states occurs, though we emphasize that there are no calculations appropriate to these conditions. Our $`\rho `$ data are consistent with activated behaviour for $`n`$ just below $`n_0`$, but deeper into the insulating state VRH becomes dominant and so we would expect to see a strong diminution of $`S^g`$ in this region. Figure 6 shows experimental data on $`S`$ as a function of $`n`$ at fixed temperatures of 2 K, 3 K and 4 K. We see that $`S`$ rises as $`n`$ decreases but crosses $`n_0`$ smoothly, i.e., we no longer see divergent behaviour of $`S`$ at $`n_0`$ as we did for $`S^d`$. Indeed there is no feature that indicates that anything unusual occurs at $`n_0`$. This behaviour is consistent with activated conduction below $`n_0`$. However, we continue to see a rise in $`S`$, and presumably $`S^g`$, to the lowest densities, which is not expected from our argument above. This can only be understood if significant activated conduction is also present at all densities. In this light, the fact that our $`\rho `$ data at lowest $`n`$ in Fig. 2, (and the data of others) appear to follow the Efros-Shklovskii VRH model so well over a limited range of temperatures appears to be somewhat coincidental.

## IV Summary and Conclusions

Generally speaking, the main focus of attention with regard to the MIT in 2D systems has been whether the metallic behaviour is the result of a transition to a new state induced by the strong disorder and possibly electron-electron interactions, or whether it is the result of a conventional physical mechanism that has yet to be unambiguously identified.

On the metallic side, our analysis has mainly concentrated on data at low $`T`$, essentially below the region of rapid $`T`$ dependence of $`\rho `$. Our main result is that the critical behaviour of $`\sigma _0`$ and $`S^d`$ of Si-MOSFETs in this region are surprisingly consistent with equations that are formally equivalent to those describing a 3D Anderson MIT. On the insulating side, the behaviour of $`\rho `$ and $`S^d`$ over a wider range of $`T`$ are also consistent with such a scenario. However, not all features that we see can be understood in this way. In particular, although the Anderson model exhibits scaling behaviour in both $`\sigma (T)`$ and $`S^d(T)`$, it does not appear to be able to explain the large increase of $`\rho `$ with $`T`$ in the metallic region observed in Si-MOSFETs.

An Anderson transition can arise purely from disorder in 3D but general scaling results predict that the equivalent transition should not occur in 2D. Still, given that many features of the data do mimic an Anderson transition, a key question in the present work would then be whether a ‘mobility edge’ actually exists, in particular at low or zero temperature. If we look for more conventional explanations of the data, a point to bear in mind is that $`k_Fl0.3`$ at $`n_0`$, where $`l`$ is the electron mean free path estimated from the conductivity, which implies very strong disorder. Under these conditions it seems plausible that a sufficiently rapid drop in $`\rho `$ with decreasing $`n`$ could be mistaken for a mobility edge, but it might actually be caused by a smooth, though very rapid, transition from weak to strong (exponential) localization of the carriers in keeping with the scaling model. If this is the explanation, then it remains to be shown how this rapid drop can mimic a critical behaviour over a relatively wide range of $`n`$.

In our analysis, the apparent divergence of $`S^d`$ as $`n_0`$ is approached from above is basically a reflection of the rapid drop in $`\sigma _0`$, a result that might remain valid even if no Anderson transition occurred. An unambiguous indicator of a mobility edge would be saturation of $`S^d`$ as a function of temperature. Our data near $`n_0`$ show no visible indication of this, but the presence of phonon drag is a complicating factor here and might mask any saturation. Clearly it would be an advantage to suppress phonon drag and reveal $`S^d`$ over a wider range of $`T`$, but this is difficult to do in 2D. In principle, it could be done by the use of very thin substrates. We do not see any features in the overall behaviour of $`S`$ which appear to correlate with the strong $`T`$ dependence of $`\rho `$ in the metallic region.

Finally, an important point that must be stressed is that to reliably identify any observed critical behaviour with a MIT requires data in the zero $`T`$ limit. Although our analysis of $`\sigma _0`$ and $`S^d`$ is based on an extrapolation to zero $`T`$, the actual data extend only to 0.3K. If the strong drop of $`\rho `$ is caused by conventional physics, then we are examining the low temperature limit of this mechanism, but not necessarily the true low temperature limit of the system. This is true of practically all the experimental data published so far.

In spite of these reservations, it is clear that the observed critical behaviour in both $`\sigma _0`$ and $`S^d`$ indicates an unknown but interesting physics of strongly disordered and interacting systems. The results open up a new window on these systems that further constrains any theoretical model proposed to explain the MIT, whether such a transition be a quantum property or simply the result of a conventional physical mechanism.

## V Acknowledgements

We acknowledge the support of the NSERC of Canada, and from INTAS, NSF, RFBR, Programs ‘Physics of nanostructures’, ‘Statistical physics’, ‘Integration’ and ‘the State support of the leading scientific schools’. We are also grateful for the very useful comments by Rudolf Römer and Cristine Villagonzalo at Technische Universität, Chemnitz.


format: text

text:
# Molecular dynamic simulation of a homogeneous bcc →hcp transition

## I Introduction.

The topic of precursor phenomena in martensitic transformations has been the subject of numerous studies over the past two decades. Martensitic systems undergo first-order transitions which are characterized by a lattice strain (sometimes accompanied by a “shuffling” of the atoms associated with a particular phonon mode). Experimentally, it has been observed from various diffraction measurements that streaking of the Bragg spots occur in the parent phase close to the transformation, indicating an incipient instability towards the formation of the martensitic phase. There is also evidence that interesting anomalous dynamical behavior occurs concurrently: the experimentally observed phonon peaks become very broad with enhanced scattering occurring at low frequencies. More puzzling, the intensities and shapes of these peaks exhibit unexpected variations between different Brillouin zones. While most of the experimental data were obtained from alloy systems where the effects of stoichiometry variations may be of importance, similar behavior have also been observed in pure metals.

In this paper, we explore the relationship between anomalous phonon behavior and the phase transformation process through molecular dynamics simulations of Zr as a prototypical system, exhibiting a bcc to hcp transformation. By monitoring the behavior of the dynamic structure factor of the system, we examine the microscopic origin of the anomalous phonon behavior in these systems. We also study the dynamics of microstructure formation during the transformation. To clarify the problem, we will concentrate on homogeneous systems without defects. The effects of inhomogeneity on the dynamical behavior of the transition will be left to a later study.

Scattering lineshapes in first-order structural phase transformations have been previously studied using various analytic treatments and simulations. Such studies, for both first-order and second-order transformations, have focused on the one-phonon contribution to the scattering, which does not exhibit a Brillouin zone dependence, even in a full anharmonic treatment. Furthermore, most of these studies have focussed on models that do not exhibit a true elastic transformation: in the high temperature phase, the linearized phonon frequencies behave as $`\omega ^2\omega _0^2+c^2k^2`$ as $`k0`$, rather than $`\omega ck`$. Such models are more appropriate for transformations such as $`bcc\omega `$, characterized by a “freezing in” of a phonon mode with a finite frequency (and wavevector), rather than the development of a static amplitude of an elastic degree of freedom that occurs in many Martensitic transformations (including the $`bcchcp`$ transformation).

Another issue with previous work examining the dynamic response of these systems is that most work examines regions of the product (low temperature) phase embedded in the high temperature matrix. This product phase is hypothesized to be due to either defects, or due to long-lived “heterophase” fluctuations from the parent phase. Such studies are not surprising, as long-lived excitations are necessary to produce “central peaks” in the scattering. However, in first-order transformations, the evidence appears to indicate that there is no clear central peak separate from the phonon peaks (and no long-lived fluctuations into the product phase), especially in systems with little disorder.

In a previous paper, which we shall refer to as paper I, we demonstrated that in a homogeneous system, fluctuations in the high-temperature bcc phase can produce “anomalous” scattering, with significantly different scattering in equivalent Brillouin zones. Neutron scattering reveals different line shapes between different Brillouin zones. Thermal diffuse scattering (energy integrated scattering), observed in both X-ray studies and neutron diffraction experiments , is not featureless, as might be expected from usual arguments concerning multiple phonon scattering, but instead shows characteristics that imply a connection with the incipient Martensitic transition. That such characteristic scattering can occur in the absence of static or long-lived fluctuations into the product phase is somewhat surprising, yet is borne out by experimental work. In paper I, we argued that these are due to fluctuations toward, but not into, the low temperature phase.

More recent molecular dynamics simulations have examined the formation of microstructure during a homogeneous bcc$``$hcp transformation, utilizing a simple EAM model for zirconium, as well as the effect of stress on the microstructure. These simulations demonstrate a number of key points concerning the formation of microstructure. Of particular importance relevant to the work presented here, they have demonstrated that fcc regions form at the junction of three hcp variants sharing a common $`11\overline{2}0`$ direction, rather than forming a triple point between twin boundaries, and that the twin regions contain basal plane stacking faults that maintain the 60 orientation between hcp variants. The transformation is interpreted to occur via the Nishiyama-Wassermann mechanism (described in Section III below), rather than the Burgers mechanism.

Here, we present similar simulations, with the goal of examining the dynamics during the bcc$``$hcp transition, with a focus on the aspects relevant to the neutron and x-ray scattering experiments. Again, we use molecular dynamics as a tool, to examine the microscopic processes associated with the dynamics. We examine the dynamics in real space, to examine the growth process and the final microstructure. The crystallography of the resulting twinned structure corresponds to experimental results, and those of Pinsook. With respect to the microstructural development, the most significant difference between the work presented here and that of Ref.  is the choice of simulation cell orientation: unlike their simulations, we do not select out one of the bcc directions to be along one of our simulation axes, as that favors hcp variants with $`11\overline{2}0`$ along this direction.

We also examine the time-dependent structure factor $`S(𝐐,t)`$, which is the equivalent of time-dependent x-ray scattering during the transition. By doing this, we directly demonstrate that the asymmetry in the scattering between different Brillouin zones is strongest during the transition, and is essentially due to the formation and motion of new Bragg peaks during the transition.

## II Details of simulations.

We performed the simulations using the same EAM potential as in the previous paper, and described in detail in Ref. . This potential has a stable bcc phase, with phonons in good agreement with experimental results. The difference between the cohesive energies between the fcc and the hcp structures is small, compared with ab initio results, as is common for these types of potentials. The EAM potentials for Ti and Zr are prone to this problem, due to the lack of directional bonding in these descriptions of these materials. Another important caveat is that this potential does not accurately reproduce the $`\omega `$-phase properties of Zr or Ti: within this potential, the $`\omega `$-phase structure is mechanically unstable. As a result, dynamics associated with large-amplitude fluctuations towards this structure (corresponding to the bcc q=$`\frac{2}{3}[111]`$ longitudinal phonon) are not accurately reproduced. In Zr and Ti (and especially in certain alloys, including ZrNb), the large amplitude fluctuations of this will be anharmonic and possibly strongly non-linear, associated with the corresponding $`bcc\omega `$ transformation. This latter transition, while not studied in this work, is an important topic. The work presented here is essentially a model system, examining the dynamics of the bcc$``$hcp transformation in the simple case where no competing transformation occurs.

The simulations were performed using two different cubic system sizes: $`18\times 18\times 18`$ conventional bcc unit cells (11664 atoms), and $`36\times 36\times 36`$ unit cells (93312 atoms). Except as noted, all results presented will be those of the larger system. No qualitative changes occurred as a result of the change in the system size. Periodic boundary conditions were used in all cases. In order to minimally interfere with the dynamics during the transition, the simulations were performed under constant energy conditions, at an energy corresponding to the average energy for the bcc phase at 1500 K. During the bcc$``$hcp transition, the system lowered its potential energy, as can be seen in fig. 1. With constant energy simulations, this produces a slight increase in temperature. For both the $`18\times 18\times 18`$ system and the $`36\times 36\times 36`$ system, this increase was about 40 K. This is a minor change in temperature, and we do not believe that this affects the dynamics of the transition significantly. We have also used constant volume simulations, with a bcc lattice constant of 3.68 Å as in paper I. The transformed structure does not develop any significant macroscopic stress. The simulations were performed for 40,000 time steps, sufficient to observe the transition and subsequent equilibration. This will be demonstrated in the following section.

## III Development of microstructure during the transition.

We begin with a review of the crystallography of the bcc$``$ hcp transition. This is most commonly stated in terms of the Burgers relationship,

$$(110)_{\mathrm{𝑏𝑐𝑐}}(0001)_{\mathrm{ℎ𝑐𝑝}}[\overline{1}11]_{\mathrm{𝑏𝑐𝑐}}[11\overline{2}0]_{\mathrm{ℎ𝑐𝑝}}.$$
(1)
The transformation occurs by “shuffling” of alternate $`(110)`$ planes in the bcc lattice, simultaneous with a shearing of the $`(1\overline{1}2)`$ planes along the $`[\overline{1}11]`$ direction. (In this paper, we will commonly use three indices for directions and planes in the bcc lattice, and four for directions and planes in the hcp lattice. We will indicate the lattice structure explicitly where we wish to emphasize the lattice correspondence, but otherwise will omit these labels.) The $`(1\overline{1}2)`$ direction then becomes a prism plane in the hcp lattice. With these dynamics in mind, we may rewrite the Burgers relationship as

$$(1\overline{1}2)_{\mathrm{𝑏𝑐𝑐}}(1\overline{1}00)_{\mathrm{ℎ𝑐𝑝}};[\overline{1}11]_{\mathrm{𝑏𝑐𝑐}}[11\overline{2}0]_{\mathrm{ℎ𝑐𝑝}}.$$
(2)
The shearing of the lattice along $`\{1\overline{1}2\}`$ planes plays an important role in understanding the anomalous scattering, as will be shown in section 4. Note that if the $`[11\overline{2}0]_{\mathrm{ℎ𝑐𝑝}}`$ direction is along the $`[111]_{\mathrm{𝑏𝑐𝑐}}`$ direction, then there are three different $`\{110\}_{\mathrm{𝑏𝑐𝑐}}`$ planes that can become hcp basal planes, namely the $`(10\overline{1})`$, $`(01\overline{1})`$ and $`(1\overline{1}0)`$ planes. With four equivalent $`111_{\mathrm{𝑏𝑐𝑐}}`$ directions in the bcc lattice, there are 12 hcp variants that may occur. In our simulated microstructure (described below), the system spontaneously selects out one of the $`111_{\mathrm{𝑏𝑐𝑐}}`$ directions, and forms three hcp variants from the associated $`\{110\}_{\mathrm{𝑏𝑐𝑐}}`$ planes.

We now establish that there is a transition during our simulations. The most direct approach is to examine the potential energy vs. time during the transition. This is shown in fig. 1, for both the $`18\times 18\times 18`$ and $`36\times 36\times 36`$ system sizes. As can be seen in the figure, both systems spend considerable amounts of time (10,000–20,000 time steps) fluctuating about a value near -6.032 eV. This corresponds to the value obtained at 1500 K for the bcc phase. Both systems then have a rapid drop in the potential energy. This drop continues for some time, then levels off near the value appropriate for the hcp system. Thus, the period during which the potential energy changes systematically serves to approximately identify when the transformational dynamics occur.

Before describing the dynamics of microstructure formation, we present the final microstructure in the system. This is shown for the $`36\times 36\times 36`$ system in fig. 2. Again, results for the smaller system are qualitatively the same. The figure shows the view down the $`[1\overline{1}\overline{1}]`$ direction. There are three different hcp variants present; the $`[1\overline{1}\overline{1}]`$ direction forms a $`11\overline{2}0`$ direction in each of these variants, in accordance with the Burgers relationship. In the projection shown in fig. 2, the system appears to be very well ordered. Each of the variants forms a “needle-like” domain, extended indefinitely (via the periodic boundary conditions) along the $`[1\overline{1}\overline{1}]`$ direction. By coloring the atoms according to their potential energy, we have highlighted the twin boundaries between different hcp variants. These boundaries form along hcp $`10\overline{1}1`$ planes, in agreement with experimental observations. As in previous simulations,, some of the variants have numerous basal plane stacking faults, allowing for the average angle between hcp domains to be 60. However, specific twin boundaries deviate noticeably from this value.

We note that in the figure, it is possible to travel from variant II to variant III without passing through such a twin boundary. In fact, there is an fcc region formed where these variants intersect, as indicated in the schematic in the figure. This region is stable throughout the simulation, and is observed in both the $`18\times 18\times 18`$ and $`36\times 36\times 36`$ systems. It is not required to minimize the potential energy of the final configuration: a lower energy configuration would arise from transforming the variant III regions to variant I, and the fcc regions to variant II. This would eliminate the total amount of twin boundary, and would also eliminate the unfavorable fcc region. We note that the fcc-hcp energy difference is small in this potential, relative to ab-initio calculations of pure Zr, suggesting that fcc formation in Zr may be decreased or eliminated in the experimental system. However, some observations of fcc Ti have been made in association with the bcc$``$hcp transformation, as discussed in ref. .

We now examine the time evolution of the system, as it progresses from the bcc to hcp phases. This is shown for the $`18\times 18\times 18`$ system in fig. 3. We focus on the smaller sized system for clarity; the results for the $`36\times 36\times 36`$ system are quite similar. In the figure, we show a slice of the system perpendicular to the  direction. The atomic positions have been averaged over 2000 time steps, in order to eliminate the high frequency atomic motions, and reveal the slower dynamics associated with the transformation. The atoms have been colored according to their time-averaged potential energy. Again, this highlights the development of the microstructure.

The first frame of fig. 3, shows the results for time steps 18,000-20,000. Before this, there is little difference between the structure and that of a perfect bcc lattice. This is consistent with fig. 1, which shows the potential energy of the $`18\times 18\times 18`$ system beginning to decrease near step 20,000. In the first frame, the start of the transformation can be seen just left of the center of the system. The vertical rows of atoms show considerable bending around this region, indicating the large accommodation strain required for the nucleus. By the second frame (steps 20,000-22,000) the transforming region is quite clear. In steps 22,000-24,000, the transformed region is starting to develop sharp boundaries, as it extends along the $`[111]`$ direction. This transformation progresses and the boundaries develop through step 26,000; by steps 26,000-28,000, the transformation is nearly complete, with only the final stages of boundary development remaining. After step 30,000 (not shown), the structure is essentially static.

## IV Simulation results for the structure factor.

One of our main objectives is to understand the origin of the pre-martensitic behavior as seen in neutron and x-ray scattering experiments. One of the most notable anomalous behaviors is the difference in thermal diffuse scattering between the so-called “$`\omega `$” points, with strong scattering located at $`𝐐=\frac{1}{3}774`$, and significantly less scattering at $`𝐐=\frac{1}{3}558`$. These correspond to equivalent points in the Brillouin zone, with equivalent magnitudes. Furthermore, in a simple one-phonon picture, the dependence of the scattering on the phonon polarizations are also identical. Therefore, from the one-phonon picture, the scattering at these points should be identical. These modes correspond to the $`\frac{2}{3}111`$ longitudinal phonon, which is the mode that leads to the development of the $`\omega `$-phase. This connection between the anisotropy and possible precursors is very suggestive. However, as shown in paper I and elaborated on below, a significant difference between these Q vectors is related to the development of the hcp phase, and the dynamics of the transformation to this phase.

We begin by showing the structure factor $`S(𝐐)`$ calculated for a smaller system, $`12\times 12\times 12`$ unit cells, which does not transform. The results for scattering in the $`(1\overline{1}0)`$ plane are shown in fig. 4. This has been calculated in the same manner as the results in paper I. This figure shows that within our simulations, the anomalous behavior is reproduced, with the scattering near $`𝐐=\frac{1}{3}774`$ significantly greater than that near $`𝐐=\frac{1}{3}558`$. In our potential, the $`\omega `$ phase is higher in potential energy than the bcc phase, and is mechanically unstable. Therefore, the results we show here can not be attributed to $`\omega `$ phase fluctuations. The anomalous scattering can be seen as a ridge extending from the  Bragg peak down to the $`\frac{1}{2}[552]`$ N-point, where a small peak in the scattering occurs.

We note that in these simulations of the bcc phase, there is no evidence of long-lived fluctuations. The pair correlation function $`g(r)`$ exhibits clear peaks associated with the bcc phase, but none associated with the hcp phase. However, some short-lived fluctuations toward the hcp phase are observed, and are correlated with anomalous scattering (as described in paper I). This is in accordance with both experiments on the high temperature phase and theory. Thus, the anomalous scattering can not be ascribed to previously postulated long-lived fluctuations into the hcp phase.

In order to explore the scattering behavior, we examine here the time-dependent scattering at these “strong” and “weak” $`\omega `$ points in the Brillouin zone. As paper I demonstrated, the difference in scattering between these two points can not be found within the one-phonon contribution, and the full scattering function must be calculated. Moreover, as there are 12 equivalent wave vectors corresponding to each of these points, we must follow the dynamics of each point.

Figure 5 shows $`S(𝐐,t)`$ for $`𝐐=\frac{1}{3}774`$, for the $`18\times 18\times 18`$ system. As indicated previously, the transformation occurs essentially between steps 20,000 and 30,000. In the figure, we see that for three of the wave vectors, intense scattering occurs during the transformation. The scattering is as large as 20 times that of the other wave vectors, or of the same wave vectors before and after the transformation. For the other nine wave vectors, the scattering shows no distinctive behavior before, during, or after the transformation. Figure 6 shows $`S(𝐐,t)`$ for $`𝐐=\frac{1}{3}558`$. In this case, there is no distinctive behavior during the transformation. For three of the wave vectors, the scattering increases by roughly a factor of five once the transformation is complete.

We can understand these results by examining the full scattering $`S(𝐐)`$ in the $`(1\overline{1}0)`$ scattering plane as a function of time. Rather than perform the full calculation, we show the scattering calculated from the time-averaged positions, such as those shown in fig. 3. We show these calculations in figs. 7 and 8. Note that these figures contain essentially static information about the microstructure. Before the transformation is initiated (steps 14,000-16,000), the structure factor shows the Bragg peaks expected for the bcc lattice, plus some weak scattering near the N-point phonons at $`𝐐=[\frac{5}{2}\frac{5}{2}1]`$ and at $`𝐐=[\frac{5}{2}\frac{5}{2}3]`$. Some smearing of the Bragg peaks along the $`[11\overline{2}]`$ direction is evident. In steps 16,000-18,000, the scattering at the N-points has intensified.

The most dramatic change to the scattering (and the most important to understanding the anomalous scattering) is seen in steps 18,000-20,000. The N-point phonons are now developing into static Bragg peaks, and the Bragg peaks are now shifted from their bcc lattice positions. The shearing along the $`[11\overline{2}]`$ direction is evident. This shear brings scattering towards the $`𝐐=\frac{1}{3}774`$ scattering vector, while leaving little scattering in the $`𝐐=\frac{1}{3}558`$ direction. This is maintained through step 26,000 - essentially, until the transformation is nearly complete. The Bragg peaks now are clearly arranged in a hexagonal pattern. Thus, the anomalous large scattering near $`𝐐=\frac{1}{3}774`$ is not from phonon scattering from $`\frac{2}{3}[111]`$ phonons, but from essentially static displacements associated with N-point phonons combined with shears in the $`[11\overline{2}]`$ direction. These are exactly the displacements needed for the formation of the hcp phase. From step 26,000 on, the primary development is the narrowing of the peaks along the $`[111]`$ direction, corresponding to the formation of long, needle-like domains in this direction as seen in the real space pictures in fig. 3. The Bragg peaks also settle into their final positions, with little remnant scattering in the $`𝐐=\frac{1}{3}774`$ or the $`𝐐=\frac{1}{3}558`$ directions.

We also note that during the course of the transformation, a Bragg peak develops near $`𝐐=\frac{4}{3}[111]`$. This is a result of the differing hcp variants (and the fcc region, as well) developing nearest neighbor lattice vectors along the  direction. For real systems (i.e. Ti and Zr) the presence of this developing Bragg peak clouds the issue of interpreting the dynamic behavior of the high temperature phase at this wave vector, as anomalies can be associated either with the $`bcchcp`$ or the $`bcc\omega `$ transformation.

## V Discussion.

The main purpose of this paper has been to explore the dynamics of the bcc$``$hcp transformation, and to relate this to the anomalous diffuse scattering observed in the bcc phase. By examining the time-dependent scattering during the scattering, we have shown that the Brillouin-zone dependence of the scattering is greatest during the transformation, and is due to the formation and shearing of the Bragg peaks as they form the hcp reciprocal lattice. The peaks move through the $`𝐐=\frac{1}{3}774`$ region of reciprocal space, creating intense scattering at this wave vector, but their final position is not at this point. There is no such movement through the $`𝐐=\frac{1}{3}558`$ point; hence the difference in scattering at these points. A one-phonon picture, taken from the bcc lattice, would predict identical scattering at these points, even in a fully anharmonic theory. Previous work has relied upon this approximation.

We may now clearly explain the anomalous diffuse scattering along the ridge between the $`𝐐=[222]`$ bcc Bragg peak and the N-point $`𝐐=\frac{1}{2}[552]`$ position. This scattering is due to regions of the bcc lattice fluctuating towards the hcp lattice, with the result that intensity develops at the N-point, and shifts along the $`[\overline{1}\overline{1}2]`$ direction. These are precisely the modes required to take the bcc lattice towards the hcp lattice, as in Eq. 2. This interpretation is consistent with previous simulations of the bcc lattice.

We have also examined the microstructure development during the transformation. In the beginning of the simulations, there are four equivalent $`111`$ directions in the bcc lattice. (Previous simulations have chosen periodic cells that break the symmetry between these directions.) The final microstructure corresponds to domains of three different hcp variants, sharing a common $`11\overline{2}0`$ direction along one of the original $`111`$ directions. These domains are extended along this direction as well (corresponding to “needle-like” domains). In addition, there is a domain in the metastable fcc lattice. The final domain structure does not minimize the potential energy; a lower energy microstructure with two hcp domains would satisfy the constraints imposed by the cubic simulation cell. Thus, the domain structure must arise out of the dynamics of the transformation process.

The hcp variants are separated by $`\{10\overline{1}1\}`$ twin boundaries, in accordance with experimental observations. These boundaries lie approximately along the bcc $`\{112\}`$ planes, as shown in fig. 2. However, in some regions, these boundaries are not flat, caused by basal plane stacking faults within the hcp variants. Pinsook has discussed the formation of these stacking faults in terms of “transformation plasticity”,, and argues that these are necessary for the boundaries to lie parallel to the bcc lattice’s $`\{112\}`$ planes, which intersect at 60 degrees, instead of the slightly different angles formed by an ideal $`\{10\overline{1}1\}`$ hcp twin boundary (55.9 for an ideal $`c/a`$ ration). However, in our simulations, the boundaries often significantly deviate from being parallel to the bcc $`\{112\}`$ planes, indicating that more stacking faults occur than would be indicated by these geometric arguments. In both our simulations and those in ref. , the hcp basal plane stacking fault energy is lower than expected for Zr, which makes comparison with the experimental system difficult.

The dynamics of the transformation correspond to a localized nucleation event, followed by growth of the domains. This is seen quite clearly in fig. 3. Such a picture is inconsistent with recent simulations of a dislocation-nucleated bcc$``$hcp transformation, in which a bcc $`𝐛=100\{011\}`$ dislocation nucleates the transition via a long wavelength distortion that evolves into a set of twin boundaries. However, this latter work is dependent on the particular dislocation; similar work with a $`𝐛=1/2110\{1\overline{1}0\}`$ dislocation showed different dynamics. This work is also different from that of Pinsook and Ackland, who observed plate-like geometry very quickly from a quenched simulation. In this latter work, the transition is homogeneous; unlike the work presented here, however, one of the $`111_{\mathrm{𝑏𝑐𝑐}}`$ directions is singled out in the simulation cell. The effect of this on the dynamics is unclear. They also observed an fcc region, but in their simulations, this was a transient phase that was eliminated during the subsequent microstructural evolution. The presence of this region is likely to be sensitive to the choice of interatomic potential; for the work presented here (and, to a lesser extent, the potential used in Refs. ), the difference between the cohesive energies of the fcc and hcp phases is small compared to the actual value in Ti and Zr, making the appearance here more likely.

The present work provides a clear demonstration of the relationship between the dynamics of the bcc$``$hcp transformation, and the anomalous scattering observed in the parent phase. The observed zone-dependence does not require defects or static transformed (or distorted) regions in the lattice; this is clear from the simulations presented here and previously,, and from experimental studies of pure metals. This is an important step in understanding the experimental results on the pre-transition dynamics, as measured through scattering experiments. Further work is required to understand how these dynamics are affected by defects, compositional fluctuations in alloy systems, and the competing $`\mathrm{𝑏𝑐𝑐}\omega `$ transformation.

###### Acknowledgements.

JRM would like to thank Yuri Gornostyrev and M. Katsnelson for useful discussions. Part of this work was made possible by the Scalable Computing Laboratory, which is funded by Iowa State University and Ames Laboratory. Ames Laboratory is operated for the U. S. Department of Energy by Iowa State University under Contract No. W-7405-Eng-82. This work was sponsored by the Director for Energy Research, Office of Basic Energy Sciences, by the High Performance Computing and Communications Initiative, by the Division of Materials Sciences, Office of Basic Energy Sciences, U. S. Department of Energy.