text:
# Density Induced Interchange of Anisotropy Axes at Half-Filled High Landau Levels

## Abstract

We observe density induced 90 rotations of the anisotropy axes in transport measurements at half-filled high Landau levels in the two dimensional electron system, where stripe states are proposed ($`\nu `$=9/2, 11/2, etc). Using a field effect transistor, we find the transition density to be $`2.9\times 10^{11}`$cm-2 at $`\nu `$=9/2. Hysteresis is observed in the vicinity of the transition. We construct a phase boundary in the filling factor-magnetic field plane in the regime $`4.4<\nu <4.6`$. An in-plane magnetic field applied along either anisotropy axis always stabilizes the low density orientation of the stripes.

The physics of electrons in higher Landau levels has recently drawn much attention due to the discovery of large transport anisotropies at half fillings ($`\nu =9/2,11/2`$ etc) in high quality two dimensional electron gases (2DEG) in GaAsLilly et al. (1999a); Du et al. (1999). At these filling factors, longitudinal magneto resistances Rxx of samples grown on (001)-oriented GaAs show strong maxima in the $`1\overline{1}0`$ direction and deep minima in the $`110`$ direction. The ratio of the resistance values can exceed 1000:1 in high quality samples. This large anisotropy is viewed as strong evidence for the formation of a unidirectional Charge Density Wave or so called “stripe state” around these filling factorsKoulakov et al. (1996); Fogler et al. (1996); Moessner and Chalker (1996); Fog . Calculations based on the Hartree-Fock approximation show the stripes to form as the result of the competition between the short-range attractive exchange interaction and the long-range Coulomb repulsion. Recently there appeared calculations beyond the Hartree-Fock approximationFradkin and Kivelson (1999); Fertig (1999); Rezayi et al. (1999); MacDonald and Fisher (2000); Shibata and Yoshioka (2001). Some of them propose the existence of liquid crystalline states with stripe ordering and broken rotational symmetryFradkin and Kivelson (1999); MacDonald and Fisher (2000). Employing an edge state transport mechanism, the “easy”(low resistance) axis is along the stripes whereas the “hard” (high resistance) axis is perpendicular to the stripes. This identification is supported by experimentsWillett et al. (2001a).

In spite of the strong evidence supporting the formation of a stripe state, the origin of its preferred orientation remains poorly understoodFil (2000); Takhtamirov and Volkov ; Rosenow and Scheidl (2001). Theory does not a priori provide a preferred direction for the stripes. However, all previous experiments have identified the $`110`$ direction as the “easy” axis. It is often assumed that anisotropic imperfections at GaAs/AlGaAs interface act as the native symmetry breaking potential. Anisotropic roughness has been observed in GaAs layers grown by Molecular Beam EpitaxyOrme et al. (1994). Two groups have investigated the correlation between the surface roughness and the orientation of the stripes which form a few thousand $`\AA `$ below the surface. Atomic Force Microscope images reveal the presence of roughness elongated along both the $`110`$ and the $`1\overline{1}0`$ directions. However, its correlation with the orientation of the stripes remains controversialWillett et al. (2001b); Cooper et al. (2001). Theoretical studies on the influence of a simple periodic modulation on the orientation of the stripes suggest that they prefer to align *perpendicular* to a weak potentialYoshioka ; Aoyama et al. . However, a parallel alignment may occur in the strong potential limitAoyama et al. .

Previous experiments that identified the $`110`$ direction as the “easy” axis were limited to fixed densities in the range of $`1.53.0\times 10^{11}`$cm-2Du et al. (1999); Lilly et al. (1999a); Cooper et al. (2001); Willett et al. (2001b). Furthermore, the correlation between the orientation of the stripes and the native symmetry breaking potential is difficult to deduce since the latter is likely to vary from specimen to specimen. Employing a tunable density Heterostructure Insulated Gate Field Effect Transistor (HIGFET) in our experiments, we have the unique opportunity to access a wide density regime and to change the 2DEG density in situ without affecting other parameters. Our data show that as the density of the 2DEG is raised above $`2.9\times 10^{11}`$cm-2, the “easy” axis rotates from the $`110`$ direction to the $`1\overline{1}0`$ direction. This result demonstrates that the pinning direction of the stripes is *not* unique.

Our HIGFET consists of a 600$`\mu `$m square mesa with edges along the $`110`$ and the $`1\overline{1}0`$ directions. The structure of the specimen consists of a (001) GaAs substrate, overgrown by MBE with 0.5$`\mu `$m of GaAs. A 50$`\AA `$ layer of AlAs is deposited, followed by a 4000$`\AA `$ layer of Al0.32Ga0.68As and capped by a heavily doped GaAs n+ layer, serving as a top gate. The AlAs layer at the interface with GaAs is intended to improve mobility. Sixteen Ni-Ge-Au contact pads are placed evenly along the edges of the mesa using standard optical lithography. One corner pad provides the contact to the top gate (see inset to Fig.1). Biasing the gate, we are able to continuously change the 2DEG density from $`5\times 10^9`$cm-2 to $`4.9\times 10^{11}`$cm-2. The peak mobility is $`1.1\times 10^7`$cm2/Vsec at $`2.3\times 10^{11}`$cm-2. The 2DEG density has a reproducible linear dependence on gate voltage. Only the lowest electronic subband is occupied for the whole density range. Fragile fractional quantum Hall states such as the 5/2 state and those around 3/2 are well developed for a wide range of densities, attesting to the high quality of the sample. The anisotropies at 9/2 and 11/2 are large, enabling us to identify the anisotropy axes unambiguously.

The electrical measurements are carried out in a dilution refrigerator with a base temperature of 20mK using standard low frequency (13-21Hz) lock-in techniques. An excitation current of 10nA is used to avoid electron heating. Rxx traces are taken at fixed perpendicular magnetic field while sweeping gate voltage to change density. The left and right panels of Fig. 1 show traces taken at 2.1 Tesla and 3.1 Tesla, respectively. At 2.1T, the 9/2 state occurs at $`2.3\times 10^{11}`$cm-2. At 3.1T, it occurs at $`3.4\times 10^{11}`$cm-2. The two central insets in Fig. 1 show the two orthogonal contact configurations used to measure the anisotropyCon . They are labelled $`110`$ and $`1\overline{1}0`$ according to the direction of the current flow. The $`1\overline{1}0`$ configuration is always represented by a solid line and the $`110`$ configuration by a dash dotted line. In the left panel, at 2.1T, the “easy” axis is along the $`110`$ direction, in agreement with all previous experiments. In the right panel, at 3.1T, the “easy” axis is *rotated* 90 and points now along the $`1\overline{1}0`$ direction. We conclude that the underlying stripe state has also rotated 90 and is now aligned along the $`1\overline{1}0`$ direction. This demonstrates that the pinning direction of the stripes is *not* unique. This observation is the most striking result of our experiments.

Our results indicate that the pinning mechanism of the stripes is more complex than previously assumed. A few conclusions can be drawn. First, contrary to earlier results, a $`1\overline{1}0`$ orientation of the stripes does exist in a perpendicular magnetic field. Second, assuming that the stripes are pinned by a unidirectional periodic potential, these stripes can be aligned either *parallel* or *perpendicular* to it. Although we do not know the mechanism responsible for the reorientation, we identify several consequences of changing density and discuss their effects. With increasing density, the electron wave function is squeezed and pressed harder against the interface. Therefore, the electrons experience the interface potential more strongly. On the other hand, with increasing density, screening also increases, effectively weakening the interface potential. In addition, the period of the stripes, calculated to be a few magnetic lengths at $`\nu =9/2`$ Koulakov et al. (1996), decreases with increasing density. Any of these factors could lead to the observed reorientation of the stripesAoyama et al. . Furthermore, we can not rule out the possibility of the coexistence of two orthogonal periodic potentials in our sample. Changing the period of the stripes could shift the relative effectiveness of these two potentials and may cause a reorientation.

To find the transition density for the reorientation of the stripes, we take density sweeps at fixed magnetic fields from 2.0T to 3.3T in steps of 0.1T and identify 2.7T and the corresponding density n$`=2.9\times 10^{11}`$cm-2 as the transition point at $`\nu =9/2`$. From 2.2T to 3.0T, we observe hysteresis between density up-sweeps and down-sweeps in a narrow range of filling factors ($`\mathrm{\Delta }\nu 0.2`$) around 9/2. As an example, Fig. 2a and b show the hysteresis at 2.6T and 40mK in both contact configurations. Sweep directions are indicated by arrows and the hysteretic area is hatched. The hysteresis decreases with increasing temperature and disappears altogether at 90mK while the anisotropy remains although it is weaker. These 90mK data unambiguously identify the orientation of the stripes in close vicinity to $`\nu =9/2`$. Fig. 2c shows the area of the hysteresis loop in the $`1\overline{1}0`$ configuration at 40mK as a function of magnetic field. It peaks near 2.5T and vanishes at 2.1T and 3.0T. In Fig. 2d, we plot the range of the hysteretic behavior in the $`\nu `$-B plane. This hysteresis “ellipse” is bounded by $`\nu =4.4`$ and $`\nu =4.6`$, which coincides with the filling factor range in which the stripe state is expected to exist theoreticallyFogler et al. (1996).

In general, the presence of hysteresis indicates a first order phase transition. As the system crosses the transition point, the initial phase continues to survive due to some pinning mechanism. In our sample, the phase transition corresponds to the reorientation of the stripes while impurities act as the pinning mechanism. In the $`\nu `$-B plane, the stripe state exists in a horizonal band between $`\nu `$=4.4 and $`\nu `$=4.6. The hysteresis “ellipse” is embedded within this band, indicating the transition region. To the left of the “ellipse”, the stripes have the $`110`$ orientation. To the right of the “ellipse”, the stripes have the $`1\overline{1}0`$ orientation. Therefore, the boundary separating the two phases must stretch through the whole bubble from the left side to the right side. We know from the 90mK data that the boundary intercepts the $`\nu =`$4.5 line at 2.7T. Furthermore, the orientation of the stripes along the border of the hysteresis “ellipse” can be deduced from data such as Fig. 2a and b since hysteresis always prolongs the initial phase. For example, at 2.6T, on density up-sweeps, the stripes are along the $`1\overline{1}0`$ direction, as seen from the minimum in the $`1\overline{1}0`$ configuration and the maximum in the $`110`$ configuration. This indicates that the stripes near $`\nu =4.4`$ have the $`1\overline{1}0`$ orientation. Conversely, density down-sweeps indicate the stripes near $`\nu =4.6`$ have the $`110`$ orientation. Similar identifications can be made for most other points along the border of the hysteresis “ellipse”. From these inputs, we conclude that the phase boundary must run from the lower left to the upper right of the “ellipse”, intercepting the $`\nu =4.5`$ line at 2.7T. Although the exact position of the phase boundary can not be derived from our data, a straight line as shown in Fig. 2d indicates its general position. Within the $`4.4<\nu <4.6`$ band, the stripes align along the $`110`$ direction to the left of the boundary. To the right of the boundary, the stripes align along the $`1\overline{1}0`$ direction. This analysis provides a phase diagram for the stripes around $`\nu =9/2`$. Hystereses with much smaller areas are also observed around $`\nu =11/2`$.

It remains unclear why in previous studies samples with densities inside the transition region of Fig. 2d neither show the $`1\overline{1}0`$ orientation nor the hysteresis Cooper et al. (2001); Willett et al. (2001b). It suggests the importance of the individual sample details. Along this line, the thin AlAs layer at the interface may play a role in establishing the specific transition density.

The application of an in-plane magnetic field has become a useful tool in investigating the native symmetry breaking potential since it gives us additional control over the orientation of the stripes. In previous fixed density samples, the stripes always aligned perpendicular to a sufficiently strong ($``$0.5T) in-plane field BipPan et al. (1999); Lilly et al. (1999b). This means that an in-plane field parallel to the Bip=0 orientation of the stripes reorients the stripes to the perpendicular alignment.

We apply an in-plane magnetic field Bip to the 2DEG by rotating the HIGFET *in situ* in our dilution refrigerator. Two cooldowns are required to align Bip along the $`110`$ or the $`1\overline{1}0`$ direction respectively. The rotator assembly is able to reach a base temperature of 70mK. Tilting angles were determined from a calibration of the rotator and from positions of well defined QHE and FQHE minima. The accuracy is better than 1. By sweeping density at a fixed angle, i.e. fixed perpendicular magnetic field, Bperp, and fixed Bip, we determine the orientation of the stripes at $`\nu =9/2`$. Fig. 3 shows the phase diagram in the Bperp-Bip plane constructed from such data. The left panel represents data taken with Bip along the $`110`$ direction. The right panel represents data taken with Bip along the $`1\overline{1}0`$ direction. A $`110`$ orientation of the stripes is represented by a hollow triangle and a $`1\overline{1}0`$ orientation by a solid square.

In order to assess the properties of the phase diagram, we first focus on the Bip=0 axis which simply reflects the reorientation of the stripes discussed earlier: below Bperp=2.7T (equivalent to n$`=2.9\times 10^{11}`$cm-2, right axis), the stripes have the $`110`$ orientation whereas above 2.7T, the $`1\overline{1}0`$ orientation prevails. As we increase Bip along either the $`110`$ or the $`1\overline{1}0`$ direction, the $`110`$ to $`1\overline{1}0`$ transition moves to higher densities. The density shift indicates a Bip induced energy gain of the $`110`$ phase. It can be viewed as a sum of two terms: an *isotropic* term and an *anisotropic* term. The isotropic term is independent of the direction of Bip. The anisotropic term reflects the difference in energy gain between a Bip parallel to the stripes and a Bip perpendicular to the stripes. Two conclusions can be drawn from the phase diagram. First, the roughly symmetric shape of the phase boundary indicates that in our sample, the dominant in-plane magnetic field effect is an isotropic energy gain of the $`110`$ phase. The anisotropic term only becomes significant for $`\left|B^{ip}\right|>0.4T`$. Second, in the high Bip regime, in the right panel, the transition to the $`1\overline{1}0`$ phase occurs at lower densities than in the left panel. From this, we conclude that stripes oriented *parallel* to Bip are energetically favored as compared to stripes oriented perpendicular to Bip. Both conclusions are in contrast to previous fixed density samples where the *anisotropic* term, which favors a *perpendicular* alignment of the stripes with respect to Bip, dominatesPan et al. (1999); Lilly et al. (1999b); Jungwirth et al. (1999); Stanescu et al. (2000); Two . Together with previous results, we must conclude that the reaction of the stripes to an in-plane magnetic field is non-universal and sensitively depends on the details of individual samples. The unique density tunability of a HIGFET is essential to such a systematic study. In fact, diverse behaviors are seen even within one specimen. As an inset to Fig. 3, we show the phase diagram of the 11/2 state, constructed during the same cooldowns as the 9/2 state using the same method. The Bip=0 transition also occurs at$``$2.7T. We can clearly see that the isotropic behavior observed at $`\nu `$=9/2 is missing. The differences between the behaviors of the 9/2 state and the 11/2 state strongly attest to the subtlety of the interactions involved. Examining other half-fillings, we note that the isotropic states of $`\nu =5/2`$ and $`\nu =7/2`$ become anisotropic in an in-plane field and the “hard” axis is always along the in-plane field direction. This is consistent with previous experimental findings and theoretical calculationsPan et al. (1999); Lilly et al. (1999b); Jungwirth et al. (1999).

In summary, we have observed a density induced reorientation of the stripe state at $`\nu =9/2`$. As the density of the 2DEG is raised above $`2.9\times 10^{11}`$cm-2, the stripes rotate from the $`110`$ orientation to the $`1\overline{1}0`$ orientation. Our results demonstrate that the pinning direction of the stripes is not unique.

###### Acknowledgements.

We thank A. Millis and S. Simon for helpful discussions and K. Baldwin and E. Chaban for experimental assistances. Financial support from the W. M. Keck Foundation is gratefully acknowledged.


format: text

text:
# Maximum mass and radius of strange stars in the linear approximation of the EOS

## 1 Introduction

One of the most important characteristics of compact relativistic astrophysical objects is their maximum allowed mass. The maximum mass is crucial for distinguishing between neutron stars and black holes in compact binaries and in determining the outcome of many astrophysical processes, including supernova collapse and the merger of binary neutron stars. The theoretical value of the maximum mass for white dwarfs and neutron stars was found by Chandrasekhar and Landau and is given by $`M_{\mathrm{max}}\left(\frac{\mathrm{}c}{G}m_B^{4/3}\right)^{3/2}`$ (Shapiro & Teukolsky (1983)), where $`m_B`$ is the mass of the baryons (in the case of white dwarfs, even pressure comes from electrons; most of the mass is in baryons). Thus, with the exception of composition-dependent numerical factors, the maximum mass of a degenerate star depends only on fundamental physical constants. The radius $`R_{\mathrm{max}}`$ of the degenerate star obeys the condition $`R_{\mathrm{max}}\frac{\mathrm{}}{mc}\left(\frac{\mathrm{}c}{Gm_B^2}\right)^{1/2}`$, with $`m`$ being the mass of either electron (white dwarfs) or neutron (neutron stars) (Shapiro & Teukolsky (1983)). White dwarfs are supported against of gravitational collapse by the degeneracy pressure of electrons whereas for neutron stars this pressure comes mainly from the nuclear force between nucleons (Shapiro & Teukolsky (1983)). For non-rotating neutron stars with the central pressure at their center tending to the limiting value $`\rho _cc^2`$, an upper bound of around $`3M_{}`$ has been found (Rhoades & Ruffini (1974)).

The quark structure of the nucleons, suggested by quantum cromodynamics, opens the possibility of a hadron-quark phase transition at high densities and/or temperatures, as suggested by Witten (Witten (1984)). In the theories of strong interaction, quark bag models suppose that breaking of physical vacuum takes place inside hadrons. As a result, vacuum energy densities inside and outside a hadron become essentially different and the vacuum pressure on the bag wall equilibrates the pressure of quarks, thus stabilizing the system. If the hypothesis of the quark matter is true, then some neutron stars could actually be strange stars, built entirely of strange matter (Alcock, Farhi & Olinto (1986), Haensel, Zdunik & Schaeffer (1986)). For a review of strange star properties, see (Cheng, Dai & Lu (1998)).

Most of the investigations of quark star properties have been done within the framework of the so-called MIT bag model. Assuming that interactions of quarks and gluons are sufficiently small, neglecting quark masses and supposing that quarks are confined to the bag volume (in the case of a bare strange star, the boundary of the bag coincides with the stellar surface), the energy density $`\rho c^2`$ and pressure $`p`$ of a quark-gluon plasma at temperature $`T`$ and chemical potential $`\mu _f`$ (the subscript $`f`$ denotes the various quark flavors $`u`$, $`d`$, $`s`$ etc.) are related, in the MIT bag model, by the equation of state (EOS) (Cheng, Dai & Lu (1998))

$$p=\frac{(\rho 4B)c^2}{3},$$
(1)
where $`B`$ is the difference between the energy density of the perturbative and non-perturbative QCD vacuum (the bag constant). Equation (1) is essentially the equation of state of a gas of massless particles with corrections due to the QCD trace anomaly and perturbative interactions. These corrections are always negative, reducing the energy density at given temperature by about a factor of two (Farhi & Jaffe (1984)). For quark stars obeying the bag model equation of state (1) the Chandrasekhar limit has been evaluated, from simple energy balance relations, in (Bannerjee, Ghosh & Raha (2000)). In addition to fundamental constants, the maximum mass also depends on the bag constant.

More sophisticated investigations of quark-gluon interactions have shown that Eq. (1) represents a limiting case of more general equations of state. For example, MIT bag models with massive strange quarks and lowest- order QCD interactions lead to some corrections terms in the equation of state of quark matter. Models incorporating restoration of chiral quark masses at high densities and giving absolutely stable strange matter can no longer be accurately described by using Eq. (1). On the other hand, in models in which quark interaction is described by an interquark potential originating from gluon exchange and by a density-dependent scalar potential which restores the chiral symmetry at high densities (Dey et al. (1998)), the equation of state $`P=P\left(\rho \right)`$ can be well approximated by a linear function in the energy density $`\rho `$ (Gondek-Rosinska et al. (2000)). It is interesting to note that Frieman and Olinto (Frieman & Olinto (1989)) and Haensel and Zdunik (Haensel & Zdunik (1989)) have already mentioned the approximation of the EOS by a linear function (see also Prakash et al. (Prakash, Baron & Prakash (1990)), Lattimer et al (Lattimer et al. (1990)). Recently Zdunik (Zdunik (2000)) has studied the linear approximation of the equation of state, obtaining all parameters of the EOS as polynomial functions of strange quark mass, the QCD coupling constant and bag constant. The scaling relations have been applied to the determination of the maximum frequency of a particle in a stable circular orbit around strange stars.

It is the purpose of this paper to obtain, by using a simple phenomenological approach (which is thermodynamical in its essence), the maximum mass and radius (the Chandrasekhar limits) for strange stars obeying a linear equation of state. Of course the maximum mass of compact astrophysical objects is a consequence of General Relativity and not of the character of motion of matter constituents. However, the formulae for maximum mass and radius, due to their simple analytical form, give a better insight into the underlying physics of quark stars, also allowing us to obtain some results which cannot be obtained by numerical methods. For example, from the obtained relations one can find the scaling relations for the maximum mass and radius of strange stars in a natural way.

The present paper is organized as follows. The maximum mass and radius of quark stars with a general linear equation of state is obtained in Section 2. In Section 3 we discuss our results conclude the paper.

## 2 Maximum mass and radius for strange stars in the linear approximation of the EOS

We assume that the strange star obeys an equation of state that can be obtained by interpolation with a linear function of density in the form:

$$p=a\left(\rho \rho _0\right)c^2,$$
(2)
where $`a`$ and $`\rho _0`$ are non-negative constants. $`\rho _0`$ is the energy density at zero pressure.

Such an equation of state has been proposed mainly to describe the strange matter built of $`u`$, $`d`$ and $`s`$ quarks (Gondek-Rosinska et al. (2000), Zdunik (2000)). The physical consistency of the model requires $`\rho _0>0`$.

The particle number density and the chemical potential corresponding to EOS (2) are given respectively by (Zdunik (2000))

$$n\left(p\right)=n_0\left(1+\frac{a+1}{a}\frac{p}{\rho _0c^2}\right)^{1/\left(a+1\right)},$$
(3)
$$\mu \left(p\right)=\mu _0\left(1+\frac{a+1}{a}\frac{p}{\rho _0c^2}\right)^{a/\left(a+1\right)},$$
(4)
where $`n_0`$ is the particle number density at zero pressure and $`\mu _0=\rho _0c^2/n_0`$.

The parameters $`a`$ and $`\rho _0`$ can be calculated, for realistic equations of state, by using a least squares fit method (Gondek-Rosinska et al. (2000), Zdunik (2000)). For the equations of state incorporating restoration of chiral quark masses at high densities proposed in Dey et al. (Dey et al. (1998)) one obtains the values $`a=0.463`$, $`\rho _0=1.15\times 10^{15}g/cm^3`$ and $`a=0.455`$, $`\rho _0=1.33\times 10^{15}g/cm^3`$, respectively (Gondek-Rosinska et al. (2000)). The standard bag model corresponds to $`a=0.333`$ and $`\rho _0=4\times 10^{14}g/cm^3`$ (Cheng, Dai & Lu (1998)).

From Eqs. (3)-(4) it follows that the particle number and chemical potential are related by the equation

$$\mu =\frac{\rho _0c^2}{n_0^{a+1}}n^a.$$
(5)
From the numerical studies of strange star models we know that the density profile of this type of astrophysical object is quite uniform (Glendenning (1996)). Therefore we can approximate $`nN/V`$, which leads to

$$\frac{\mu }{\mu _0}=\left(\frac{4\pi }{3}\right)^a\left(\frac{N}{n_0}\right)^aR^{3a},$$
(6)
where $`N`$ is the total number of particles in a star of radius $`R`$ and volume $`V`$.

With the use of Eqs. (2)-(6) one obtains the energy density of the star in the form

$$\rho =\frac{\rho _0}{a+1}\left(\frac{4\pi }{3}\right)^{a1}\left(\frac{N}{n_0}\right)^{a+1}R^{3\left(a+1\right)}+\frac{a}{a+1}\rho _0.$$
(7)
The total mass $`M`$ of the star is defined according to $`M=4\pi _0^R\rho r^2𝑑r\frac{4\pi }{3}\rho R^3`$ and is given by

$$M=\frac{\rho _0}{a+1}\left(\frac{4\pi }{3}\right)^a\left(\frac{N}{n_0}\right)^{a+1}R^{3a}+\frac{4\pi }{3}\frac{a}{a+1}\rho _0R^3,$$
(8)
where we assumed that the energy density is approximately constant inside the star.

Extremizing the mass with respect to the radius $`R`$ by means of $`M/R=0`$ gives the relation

$$\frac{\rho _0}{a+1}\left(\frac{4\pi }{3}\right)^a\left(\frac{N}{n_0}\right)^{a+1}R^{3a}=\frac{4\pi }{3}\frac{1}{a+1}\rho _0R^3.$$
(9)
Substituting Eq. (9) into Eq. (8) we obtain the maximum mass of the strange star in the linear approximation of the EOS:

$$M=\frac{4\pi }{3}\rho _0R^3.$$
(10)
This expression is very similar to the expression for the maximum mass of the quark star obtained assuming that the star is composed of three-flavour masslesss quarks, confined in a large bag (Bannerjee, Ghosh & Raha (2000), Cheng & Harko (2000)). From a physical point of view, Eq. (10) describes a uniform density zero pressure stellar type configuration.

The maximum equilibrium radius corresponds to a minimum total energy of the star (including the gravitational one), for any radius. For ordinary compact stars, the mass is entirely due to baryons, and the corresponding (Newtonian) gravitational potential energy is of the order $`E_G\alpha GM^2/R`$ ($`\alpha =3/5`$ for constant density Newtonian stars). For quark stars, assumed to be formed of massless quarks, the total mass can be calculated from the total (thermodynamic as well as confinment) energy in the star. One possibility for the estimation of the gravitational energy per fermion is to define an effective quark mass incorporating all the energy contributions (Bannerjee, Ghosh & Raha (2000)).

The gravitational energy per particle (the strange star is assumed to be formed from fermions) is

$$E_G=\frac{GMm_{eff}}{R},$$
(11)
where $`m_{eff}`$ is the effective mass of the particles inside the star, incorporating also effects such as quark confinment. For a star with $`N`$ particles one can write $`M=Nm_{eff}=\rho _0V`$, or $`m_{eff}=\rho _0/n`$. On the other hand one can assume $`\mu =\rho _0/2n`$ (Bannerjee, Ghosh & Raha (2000)), leading to $`m_{eff}=2\mu /c^2`$. Hence, with the use of Eqs. (6), (9) and (10) we can express the gravitational energy per particle as

$$E_G=2\left(\frac{4\pi }{3}\right)^2G\frac{\rho _0^2}{N}R^5.$$
(12)
The energy density per particle of the fermions follows from Eq. (7) and is given by:

$$E_F=\frac{4\pi }{3}\frac{1}{a+1}\frac{\rho _0}{N}R^3.$$
(13)
The total energy $`E`$ per particle is

$$E=\frac{4\pi }{3}\frac{1}{a+1}\frac{\rho _0}{N}R^32\left(\frac{4\pi }{3}\right)^2G\frac{\rho _0^2}{N}R^5.$$
(14)
Extremizing the total energy with respect to the radius (with the total particle number kept constant), $`\left(\frac{E}{R}\right)_{N=const.}=0`$, it follows that the maximum radius of the equilibrium configuration in the linear approximation of the EOS is given by:

$$R_{max}=R_0\frac{c}{\sqrt{\pi \left(a+1\right)G\rho _0}}.$$
(15)
The maximum mass of the star can be calculated from Eq. (10) and is:

$$M_{max}=\frac{4}{3}\frac{R_0^3}{\left(a+1\right)^{3/2}}\frac{c^3}{G}\frac{1}{\sqrt{\pi G\rho _0}}.$$
(16)
In Eqs. (15) and (16) $`R_0`$ is a numerical factor of the order $`R_00.474`$.

The maximum radius of the quark star given by Eq. (15) is the radius corresponding to the maximum mass. On the other hand for the existing models of strange stars, the configuration with maximum mass has a radius which is lower than the maximum radius. For example, for strange stars described by the bag model equation of state, the maximum radius is $`11.40km`$, while the radius corresponding to the maximum mass is $`10.93km`$, which is $`4\%`$ lower than the maximum radius. This difference is neglected in Eq. (15).

The maximum mass and radius of the star are strongly dependent on the numerical value of the coefficient $`R_0`$ and estimations based on other physical models could lead to different numerical estimates of the limiting values of the basic parameters of the static strange stars.

## 3 Discussions and final remarks

In the present paper we have shown that there is a maximum mass and radius (the Chandrasekhar limits) for quark stars whose equation of state can be approximated by a linear function of the density. We have also obtained the explicit expressions for $`M_{max}`$ and $`R_{max}`$.

With respect to the scaling of the parameter $`\rho _0`$ of the form $`\rho _0k\rho _0`$, the maximum mass and radius have the following scaling behaviors:

$$R_{max}k^{1/2}R_{\mathrm{max}},M_{\mathrm{max}}k^{1/2}M_{\mathrm{max}}.$$
(17)
For the maximum mass of the strange stars this scaling relation has also been found from the numerical study of the structure equations in the framework of the bag model (Witten (1984), Haensel, Zdunik & Schaeffer (1986)).

A rescaling of the parameter $`a`$ of the form $`a+1K\left(a+1\right)`$, with $`\rho _0`$ unscaled, leads to a transformation of the radius and mass of the form

$$R_{max}K^{1/2}R_{max},M_{max}K^{3/2}M_{max}.$$
(18)
A simultaneous rescaling of both $`a`$ and $`\rho _0`$, with $`a+1K\left(a+1\right),\rho _0k\rho _0`$ gives

$$R_{max}k^{1/2}K^{1/2}R_{max},M_{max}k^{1/2}K^{3/2}M_{max}.$$
(19)
The maximum mass and radius of strange stars with linear EOS is strongly dependent on the numerical value of $`\rho _0`$, the mass decreasing with increasing $`\rho _0`$. For $`\rho _0=4B`$, with the bag constant $`B=10^{14}g/cm^3`$ ($`56MeVfm^3`$) we obtain $`M_{max}=1.83M_{}`$, a value that must be compared to the value $`M_{max}=2M_{}`$ obtained by numerical methods (Witten (1984), Haensel, Zdunik & Schaeffer (1986)). The difference between the numerical and theoretical predictions is around $`10\%`$. For $`\rho _0=1.33\times 10^{15}gcm^3`$ the maximum mass of the star is about $`1M_{}`$.

Generally our formulae (15) and (16) underestimate the maximum values of the mass and radius because we have assumed that the density inside the star is uniform. It is obvious that near the surface the density is much lower than at the center of the compact object. Due to the approximations and simplifications used to derive the basic expressions, reflected mainly in the uncertainties in the exact value of the coefficient $`R_0`$, Eqs. (15) and (16) cannot provide high precison numerical values of the maximum mass and radius for linear EOS stars, which must be obtained by numerically integrating the gravitational field equations.

For the the linear EOS, $`M_{max}`$ and $`R_{max}`$ depend mainly on the fundamental constants $`c`$ and $`G`$ and on the zero pressure density $`\rho _0`$ (the bag constant). The Chandrasekhar expressions for the same physical parameters involve two more fundamental constants, $`\mathrm{}`$ and the mass of the electron or neutron.

For quark stars, usually one assume they are composed of a three-flavour system of massless quarks, confined in a large bag. Hence the mass of the quark cannot play any role in the mass formula. But the linear EOS with arbitrary $`a`$ can describe quark matter with non-zero quark masses (the mass of the strange quark $`m_s200MeV`$), forming a degenerate Fermi gas (Gondek-Rosinska et al. (2000), Zdunik (2000)). Therefore this system should also be described by the same formulae as white dwarfs or neutron stars, not only by Eqs. (15) - (16). Generally $`\rho _0`$ is a function of the mass of the strange quark, so this mass implicitly appears in the expression of the maximum mass and radius. But on the other hand we can assume that the Chandrasekhar limit also applies to quark stars with the baryon mass substituted by an effective quark mass $`m_{qeff}`$, representing the minimum mass of the quark bubbles composing the star. Hence we must have

$$\left(\frac{\mathrm{}c}{G}m_{qeff}^{4/3}\right)^{3/2}\frac{c^3}{G}\frac{1}{\sqrt{\pi G\rho _0}}.$$
(20)
Eq. (20) leads to the following expression of the effective mass of the ”elementary” quark bubble:

$$m_{qeff}\left(\frac{\mathrm{}\rho _0^{1/3}}{c}\right)^{3/4}.$$
(21)
The effective quark mass is determined only by elementary particle physics constants and is independent of $`G`$. From its construction $`m_{qeff}`$ should be relevant when the system is quantum mechanical and involves high velocities and energies. With respect to a scaling of the zero pressure density of the form $`\rho _0k\rho _0`$, the effective quark mass has the scaling behavior $`m_{qeff}m_{qeff}k^{1/4}`$, similar to the scaling of the strange quark mass (Zdunik (2000)).

For $`\rho _0=4\times 10^{14}gcm^3`$ we obtain $`m_{qeff}3.63\times 10^{25}g204MeV`$. For $`\rho _0=1.33\times 10^{15}gcm^3`$ Eq. (21) gives $`m_{qeff}4.9\times 10^{25}g275MeV`$. The mass given by Eq. (21) can be considered as the minimum mass of the stable quark bubble. It is of the same order of magnitude as the mass $`m_s`$ of the strange quark. Therefore the Chandrasekhar limit applies also for quark stars if we take $`m_{qeff}`$ for the mass of the elementary constituent of the star .

In the present paper we have considered the maximum mass and radius of strange stars in the linear approximation of the equation of state and the dependence of these quantities on the parameter $`a`$ has been found. We have also pointed out the existence of scaling relations for the maximum radius of strange stars, an aspect that has not been mentioned in previous investigations (Witten (1984), Haensel, Zdunik & Schaeffer (1986), Bannerjee, Ghosh & Raha (2000), Zdunik (2000)). Our formulae also lead to the transformation relations for the maximum mass and radius of strange stars with respect to separate and simultaneous scaling of the parameters $`a`$ and $`\rho _0`$. On the other hand the possibility of estimation of the mass of the strange quark from general astrophysical considerations can perhaps give a better understanding of the deep connection between micro- and macro-physics.

Acknowledgements. This work is partially supported by a RGC grant of Hong Kong Government and T.H. is supported by a studentship of the University of Hong Kong. The authors are very grateful to the anonymous referee whose comments helped to improve an earlier version of the manuscript.


text:
# Phase estimation as a quantum nondemolition measurement

(March 27, 2002)

## Abstract

The phase estimation algorithm, which is at the heart of a variety of quantum algorithms, including Shor’s factoring algorithm, allows a quantum computer to accurately determine an eigenvalue of an unitary operator. Quantum nondemolition measurements are a quantum mechanical procedure, used to overcome the standard quantum limit when measuring an observable. We show that the phase estimation algorithm, in both the discrete and continuous variable setting, can be viewed as a quantum nondemolition measurement.

There is currently a great deal of active research aimed at developing implementation schemes for quantum computation. For an overview of this field see Nielsen and Chuang Nielsen and Chuang (2000) or Preskill notes Preskill (1998) and references therein. The goal of quantum computation is to build a computer, which, using the laws of quantum mechanics, can out perform any classical computer. Using the laws of quantum mechanics to enhance a device’s performance is not new. In quantum optics, the procedure of quantum nondemolition (QND) measurement Caves et al. (1980); Milburn and Walls (1983); LaPorta et al. (1989); Holland et al. (1990); Poizat and Grangier (1993) aims to accurately measure an observable without perturbing it. This procedure is consistent with the uncertainty principle, as, ideally, noise is induced only in the conjugate of the observable being measured.

Much of the impetus for research in quantum computation stems from Shor’s algorithm Shor (1994), which allows a quantum computer to factor integers exponentially faster than can currently be done on a classical computer. The factoring of large integers allows the decryption of many public key encryption systems. It has been shown that Shor’s, and many other quantum algorithms, can be described in terms of the phase estimation algorithm Kitaev ; Cleve et al. (1998).

It is evident that both phase estimation and QND measurement have the same goal: they both aim to accurately measure the eigenvalue of an operator (observable). In this paper we show that they are essentially the same procedure. That is, every instance of the phase estimation algorithm can be considered a QND measurement.

The original motivation for the study of QND measurements was the need to accurately measure weak classical forces Caves et al. (1980). Such measurements will be necessary, for example, in the detection of gravitational radiation. Classically, it is possible to measure both the position and momentum of a particle with complete precision. However, quantum mechanics tells us that the product of the uncertainty in the position and momentum of a particle must be greater than Planck’s constant divided by $`4\pi `$,

$`\mathrm{\Delta }x\mathrm{\Delta }p`$ $``$ $`{\displaystyle \frac{\mathrm{}}{2}}.`$ (1)
The consequence of the uncertainty principle for an oscillator is the fact that any quantum mechanical state is characterized by an “error box” in the phase plane with area at least

$`\mathrm{\Delta }x\mathrm{\Delta }p/m\omega `$ $``$ $`\mathrm{}/2m\omega ,`$ (2)
where $`m`$ and $`\omega `$ are the mass and angular frequency respectively Caves et al. (1980). We can re-write Eq. (2) in terms of the quadrature variables, $`\widehat{X}`$ and $`\widehat{Y}`$,

$`\mathrm{\Delta }X\mathrm{\Delta }Y`$ $``$ $`\mathrm{}/2m\omega ,`$ (3)
where

$`\widehat{x}+i\widehat{p}/m\omega `$ $`=`$ $`(\widehat{X}+i\widehat{Y})e^{i\omega t}.`$ (4)
If we naïvely perform an *amplitude-and-phase* measurement of the quadrature variables, we obtain a round error box, with

$`\mathrm{\Delta }X=\mathrm{\Delta }Y`$ $`=`$ $`\sqrt{{\displaystyle \frac{\mathrm{}}{2m\omega }}}.`$ (5)
However, by performing a back-action-evading QND measurement of, say, $`\widehat{X}`$, the error box becomes a long thin ellipse, with

$`\mathrm{\Delta }X\mathrm{\Delta }Y`$ (6)
and $`\mathrm{\Delta }X`$ can be made as small as one wishes, in principle.

Although a precise definition for a QND measurement is hard to pin down from the literature, there are three basic criteria which a QND measurement must satisfy. Firstly, it is required that the QND variable $`\widehat{A}`$, which we wish to measure, commutes with the free Hamiltonian evolution of the system,

$$\begin{array}{ccc}& [\widehat{H}_F,\widehat{A}]=0,& \hfill \text{(QND 1)}\end{array}$$

so that if the system is measured in an eigenstate of $`\widehat{A}`$, it remains in this eigenstate for all subsequent times, (perhaps only in the interaction picture). The most important criterion which must be satisfied, is that the Hamiltonian, $`\widehat{H}_{SM}`$, coupling the meter variable $`\widehat{B}`$ to the system must commute with the QND variable $`\widehat{A}`$,

$$\begin{array}{ccc}& [\widehat{H}_{SM},\widehat{A}]=0.& \hfill \text{(QND 2)}\end{array}$$

This is known as the *back action evasion criterion*. Finally, it is also necessary that the meter variable $`\widehat{B}`$, does *not* commute with the interaction Hamiltonian,

$$\begin{array}{ccc}& [\widehat{H}_{SM},\widehat{B}]0.& \hfill \text{(QND 3)}\end{array}$$

If this requirement is not fullfilled, then measuring the meter variable $`\widehat{B}`$ will yield no information about the QND variable $`\widehat{A}`$.

To illustrate how a QND measurement procedure works, consider the example of an ideal quadrature phase measurement of a single mode field Walls and Milburn (1995). Suppose we have two degenerate modes $`a`$ and $`b`$, of an electromagnetic field. We wish to measure $`\widehat{X}_a`$ using the meter variable $`\widehat{X}_b`$. Now, if we assume that the amplitude quadratures of each mode can be coupled according to the interaction Hamiltonian,

$`\widehat{H}_{SM}`$ $`=`$ $`\mathrm{}\chi \widehat{X}_a\widehat{Y}_b,`$ (7)
where $`\chi `$ is the coupling strength, then clearly (QND 2) and (QND 3) are satisfied. (QND 1) is also satisfied by choosing a rotating frame of reference. In order to draw analogies with the phase estimation algorithm, we will divide the QND measurement into three steps. The first step involves initializing mode $`b`$ into the zero eigenstate of the meter variable, $`|0_b`$,

$`\widehat{X}_b|0_b`$ $`=`$ $`0.`$ (8)
Now, for simplicity, we will assume that mode $`a`$ is already in an eigenstate $`|x_a`$ of the QND variable,

$`\widehat{X}_a|x_a`$ $`=`$ $`x|x_a.`$ (9)
Thus, the initial state of the quantum system is

$`|x_a|0_b.`$ (10)
Re-writing the meter system in terms of the phase variable, $`\widehat{Y}_b`$, gives

$`{\displaystyle 𝑑y|x_a|y_b}.`$ (11)
Evolving according to the interaction Hamiltonian (Eq. (7)) results in the state

$`e^{i\chi \widehat{X}_a\widehat{Y}_b}{\displaystyle 𝑑y|x_a|y_b}`$ $`=`$ $`{\displaystyle 𝑑ye^{i\chi xy}|x_a|y_b}.`$ (12)
Finally re-writing mode $`b`$ in terms of the amplitude variable $`\widehat{X}_b`$, we have

$`|x_a|\chi x_b.`$ (13)
Therefore, by performing a quadrature measurement of the amplitude of mode $`b`$, we obtain information about the state of mode $`a`$.

In the above description we have used eigenstates of the quadrature variables, however, like position and momentum eigenstates, these are physically unrealisable. However, they can be approximated to arbitrary degree by using highly squeezed states.

The goal of the phase estimation algorithm is to obtain an eigenvalue of a unitary operator $`\widehat{U}`$. Although hybrid systems have been discussed in the literature Travaglione and Milburn (2001); Lloyd (2000), generally the phase estimation algorithm is described in terms of (discrete) qubit registers Kitaev ; Cleve et al. (1998). Discrete classical computation, is accomplished through the manipulation of binary digits (bits) using logic gates such as AND, OR, NOT and FANOUT, whereas discrete quantum computation uses arrays of two level quantum systems (qubit registers) which are manipulated via one and two qubit unitary operators. In analogy with the classical bit, it is conventional to use the symbols $`|0`$ and $`|1`$ to denote the two levels of the quantum system. Also in analogy with classical computation, we use the compact integer representation to represent the basis states of a qubit register. For example,

$`|1_3|1_2|0_1|1_0`$ $``$ $`|13.`$ (14)
For an $`m`$ qubit register, we also introduce the Hermitian operator,

$`\widehat{\mathrm{\Lambda }}`$ $`=`$ $`{\displaystyle \underset{j=0}{\overset{m1}{}}}I\mathrm{}I\underset{j^{\mathrm{th}}\mathrm{term}}{\underset{}{2^{j1}(IZ)}}I\mathrm{}I,`$ (15)
where $`I`$ is the identity operator, and $`Z`$ is the pauli operator, with matrix representation,

$`Z`$ $`=`$ $`\left[\begin{array}{cc}\hfill 1& \hfill 0\\ \hfill 0& \hfill 1\end{array}\right].`$ (18)
The computational basis states are eigenstates of the operator $`\widehat{\mathrm{\Lambda }}`$ with an eigenvalue corresponding to their integer representation,

$`\widehat{\mathrm{\Lambda }}|x`$ $`=`$ $`x|xx[0,1,\mathrm{},2^m1].`$ (19)
Another basis which is useful in understanding the phase estimation algorithm is the Fourier transformed (FT) basis. This is the basis which arises by performing the quantum Fourier tranformCoppersmith (1994), $`\widehat{Q}`$, on the computational basis,

$`|\overline{x}\widehat{Q}|xx[0,1,\mathrm{},2^m1].`$ (20)
The Hermitian operator associated with this basis is

$`\widehat{\mathrm{\Gamma }}`$ $`=`$ $`\widehat{Q}\widehat{\mathrm{\Lambda }}\widehat{Q}^1,`$ (21)
as

$`\widehat{\mathrm{\Gamma }}|\overline{x}`$ $`=`$ $`x|\overline{x}x[0,1,\mathrm{},2^m1].`$ (22)
Suppose that we have a register of qubits, the *target* register, which is in an eigenstate $`|\varphi _T`$, of some unitary operator $`\widehat{U}`$. We wish to determine the eigenvalue corresponding to $`|\varphi _T`$. As $`\widehat{U}`$ is unitary, we have

$`\widehat{U}|\varphi _T`$ $`=`$ $`e^{i\varphi }|\varphi _T,`$ (23)
where $`\varphi `$ is some number between 0 and $`2\pi `$. The phase estimation algorithm utilizes an auxillary register of qubits, which we shall call the *index* register. Fig. 1 depicts the quantum circuit diagram for the discrete phase estimation algorithm.

The algorithm begins by placing the index register in the zero of the FT basis,

$`|\mathrm{\Phi }_0`$ $`=`$ $`|\overline{0}_I|\varphi _T,`$ (24)
Assuming that each qubit in the index register begins the computation in the state $`|0`$, this step is accomplished by performing a Hadamard gate, $`\widehat{H}`$, on each qubit. The next step is to perform a series of controlled $`\widehat{U}^{2^k}`$ gates, which couple the target and index registers. For each $`k[0,1,\mathrm{},m1]`$, the gate $`\widehat{U}^{2^k}`$ is applied conditionally upon the state of the $`k`$-th qubit in the index register. The cummulative effect of applying these gates is the operation, $`\widehat{U}^{\widehat{\mathrm{\Lambda }}}`$. The state of the system after the application of $`\widehat{U}^{\widehat{\mathrm{\Lambda }}}`$ is

$`|\mathrm{\Phi }_1`$ $`=`$ $`\widehat{U}^{\widehat{\mathrm{\Lambda }}}|\mathrm{\Phi }_0`$ (25)
$`=`$ $`\widehat{U}^{\widehat{\mathrm{\Lambda }}}|\overline{0}_I|\varphi _T`$
$`=`$ $`\widehat{U}^{\widehat{\mathrm{\Lambda }}}{\displaystyle \frac{1}{\sqrt{M}}}{\displaystyle \underset{j=0}{\overset{M1}{}}}|j_I|\varphi _T`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{M}}}{\displaystyle \underset{j=0}{\overset{M1}{}}}|j_I\widehat{U}^j|\varphi _T`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{M}}}{\displaystyle \underset{j=0}{\overset{M1}{}}}|j_Ie^{i\varphi j}|\varphi _T,`$
where the last line makes use of Eq. (23). The final step in the algorithm involves measuring the index register in the FT basis, or equivalently, performing the operation $`\widehat{Q}^1`$, and measuring in the computational basis. This will, with high probability Travaglione and Milburn (2001), result in an $`m`$-bit approximation of $`\varphi `$.

Having re-cast the phase estimation algorithm in terms of the FT basis, it is now not difficult to see that the phase estimation algorithm is effectively a QND measurement. We are trying to find an eigenvalue of the unitary operator, $`\widehat{U}`$, which can of course be written as

$`\widehat{U}`$ $``$ $`e^{i\widehat{H}_U},`$ (26)
for some Hermitian operator, $`\widehat{H}_U`$. This operator is the QND variable that we are trying to measure, $`\widehat{\mathrm{\Gamma }}`$ is the meter variable, and $`\widehat{H}_U\widehat{\mathrm{\Lambda }}`$ is the interaction Hamiltonian, as

$`\widehat{U}^{\widehat{\mathrm{\Lambda }}}`$ $`=`$ $`e^{i\widehat{H}_U\widehat{\mathrm{\Lambda }}}.`$ (27)
In the quantum circuit model, it is assumed that the qubits undergo no free evolution, so (QND 1) is automatically satisfied. It is clear from Eq. (27) that the back action evasion criterion is satisfied (QND 2), and it is not difficult to see that the interaction hamiltonian does not commute with the meter variable (QND 3) as

$`[\widehat{\mathrm{\Lambda }},\widehat{\mathrm{\Gamma }}]`$ $``$ $`0.`$ (28)
To illustrate how the phase estimation algorithm acts as a QND measurement, we describe the simplest possible algorithm. Suppose we wish to determine the eigenvalue of the pauli $`X`$ operator, or NOT gate,

$`X`$ $`=`$ $`\left[\begin{array}{cc}\hfill 0& \hfill 1\\ \hfill 1& \hfill 0\end{array}\right].`$ (31)
The eigenstates of this operator (ignoring a normalization factor) are $`|0\pm |1`$. The quantum circuit diagram for this simple phase estimation algorithm is depicted in Fig. 2.

In this example the QND variable is

$`\widehat{H}_U`$ $`=`$ $`{\displaystyle \frac{\pi }{2}}(I_TX_T),`$ (32)
the meter variable is

$`\widehat{\mathrm{\Gamma }}=\widehat{Q}\widehat{\mathrm{\Lambda }}\widehat{Q}^1`$ $`=`$ $`{\displaystyle \frac{H_I(I_IZ_I)H_I}{2}}={\displaystyle \frac{1}{2}}\left[\begin{array}{cc}\hfill 1& \hfill 1\\ \hfill 1& \hfill 1\end{array}\right],`$ (35)
as the single qubit quantum Fourier transform is the self inverse Hadamard operator, $`H_I`$, and the interaction Hamiltonian is

$`\widehat{H}_U\widehat{\mathrm{\Lambda }}`$ $`=`$ $`{\displaystyle \frac{\pi }{4}}(I_TX_T)(I_IZ_I).`$ (36)
Obviously, all the QND criteria are satisfied, and the measurement proceeds by placing the meter system into the zero eigenstate of the meter variable, $`|0+|1`$, which is accomplished by the first Hadamard gate, evolving according to the interaction Hamiltonian, which is accomplished by the controlled NOT gate, and finally measuring the meter variable, by performing another Hadamard and then measuring in the computational basis.

We have just shown that the discrete phase estimation algorithm is equivalent to a QND measurement. Recently, a number of papers have been written exploring quantum computation using continuous variables Lloyd and Braunstein (1999); Gottesman et al. (2001); Ralph et al. (2001); Bartlett and Sanders (2001). We now show that continuous variable phase estimation is equivalent to continuous variable QND measurement.

As with all quantum algorithms, the phase estimation algorithm can be partitioned into three stages; the *initialization* stage, the *entangling* stage and the *measurement* stage. The initialization stage involves placing the index system into the zero eigenstate of the FT basis. The entangling stage involves the application of the operator

$`\widehat{U}^{\widehat{\mathrm{\Lambda }}}`$ $`=`$ $`e^{i\widehat{H}_U\widehat{\mathrm{\Lambda }}},`$ (37)
and the measurement stage involves measuring the index system in the conjugate of the computational basis. If we replace the discrete qubit registers with infinite level (continuous) quantum variables, and carry out each of these steps, we obtain something very similar to the ideal quadrature QND measurement discussed earlier. Suppose that the continuous variable operator, acting on the target system, whose eigenvalue we wish to determine is $`e^{i\widehat{A}_T}`$, and we choose the position of the continuous index system as the computational basis. Then following the steps of the phase estimation algorithm, we prepare our index system in the zero eigenstate of the conjugate basis, that is, the zero momentum eigenstate. This step corresponds to initialising our meter variable. Then we perform the $`\widehat{U}`$ operator, conditioned upon the computational basis, $`e^{i\widehat{x}_I\widehat{A}_T}`$, which corresponds to applying an interaction Hamiltonian which satisfies (QND 1) and (QND 2). Before finally performing a Fourier transform and measuring the position variable, which is equivalent to measuring the momentum of the index system (i.e. the meter variable). This is clearly a QND measurement of $`\widehat{A}`$.

Of course, position eigenstates are physically unrealisable, and any continuous variable implementation of the phase estimation algorithm would need to be done using approximate position eigenstates, such as squeezed states. Just as ideal QND measurements can only be approximated experimentally, ideal continuous variable phase estimation can also only be approximated.

The study of quantum computation brings together a number of disciplines including computer science, mathematics and quantum mechanics. We have taken the phase estimation algorithm, which has its roots in computer science, and shown that it is equivalent to the quantum mechanical procedure of QND measurement. We have also described a continuous variable version of the phase estimation algorithm. It is interesting to speculate on whether other algorithms such as Grover’s algorithm might have similarities with well studied quantum phenomenon.

###### Acknowledgements.

BCT and GJM thank S. Lloyd for helpful discussions.


text:
# Shear and Mixing in Oscillatory Doubly Diffusive Convection

## 1 Introduction

Doubly diffusive convection occurs when two material properties, such as heat and salt with differing diffusion rates, affect the fluid density and are stratified so that one is stabilizing and the other destabilizing (Turner, 1979). When the stably stratified property diffuses more rapidly, the onset of the convection is by way of a monotonic growth that gives rise to structures known as *salt fingers*. This kind of convection is an effective mixer, the more so because it induces large scale shears (Paparella & Spiegel, 1999, hereafter PS) in the same way that ordinary Rayleigh-Bénard convection does (Krishnamurti & Howard, 1981, Howard & Krishnamurti, 1986). On the other hand, when the destabilizing constituent is diffused more rapidly, the onset of convection is through growing oscillations. The mechanism of mixing in this case, that we shall refer to as oscillatory doubly diffusive, or ODD, convection is less well understood than in the salt finger case and this is the subject of the present work.

The question of mixing by ODD convection is not academic since this form of doubly diffusive convection does arise in natural circumstances. One example occurs in the cores of stars in late stages of their lives (Spiegel, 1969). In the version of this process called semiconvection by astronomers, the stabilizing constituent is a relatively heavy element, such as helium in hydrogen-rich stellar material. In many cases of interest, the resulting stabilizing gradient of molecular weight opposes an unstable entropy gradient. The problem that this raises is to determine the rate of mixing of the heavy element outward in the star. Though this process involves some complications peculiar to the astrophysical situation that we shall address in another place, we may usefully confine ourselves here to the exploration of the basic fluid dynamical aspects of mixing by ODD convection. Suffice it to say that if the mixing caused by semiconvection is effective, the convection can be self-regulating and this is an important ingredient in estimating the speed at which the star evolves.

Other natural occurrences of ODD convection are found in the Earth’s oceans, most notably, below the polar ice caps. There melting ice releases cold fresh water above warmer saltier water and produces the kind of situation we shall study here (Neal et al., 1969, Jacobs et al., 1981). Another natural occurrence of ODD convection is seen in *meddies*, vortices of warm and salty Mediterranean water common at mid depths in the eastern Atlantic (Ruddick, 1992). Intrusive staircases present an alternance of fingering and ODD stratification (Marmorino, 1991). Here too, there are additional complicating features that we shall not confront in studying the conventional plane-parallel situation stratified (initially) only in the vertical direction. This ideal situation is in fact the object of extensive experiments that have been carried out in binary mixtures (Kolodner et al., 1990). These have revealed much about the nature of weakly nonlinear behavior of ODD convection. Though the waves seen in these experiments will surely do some mixing, we want to argue here that, at higher Rayleigh numbers, we may expect to see induced shears like those seen by Krishnamurti & Howard (1981) in thermal convection. Those shears, in concert with the oscillatory convective motions, should lead to rapid mixing of constituents.

To bring out these results, we shall use the same modal procedure that we used in the monotonic case previously (PS). This permits surveys of the parameter space with good qualitative and semiquantitative results without the excessive demands on numerical resources that are necessary for full simulations with large aspect ratio. Modal truncations are known to yield results which depend on the number of modes retained. This is illustrated, for example, by Franceschini & Tebaldi (1985) for relatively high dimensional truncations of the Navier-Stokes equations. However, when the focus is more on physical mechanisms than on the details of the bifurcation diagram, even crude truncations have proven to be a useful guide in the exploration of the dynamics of the underlying PDE. Already in the simpler case of thermal convection, it was shown by Howard and Krishnamurti (1986) that a model consisting of a small number of ODEs captures some features of the large-scale shearing flow. In the particular case of ODD convection, a Galerkin approach was used by Veronis (1965), to show oscillatory solutions in a system of 5 ODEs (see also Da Costa et al., 1981). In a later study, Gough & Toomre (1982) used a truncation only in the horizontal direction, as had been used in the purely thermal problem (Gough et al., 1975), but they did not allow for a large-scale horizontal flow. Their results focus on the stability of an initial step-like density stratification and they omitted the horizontal shear of the Krishnamurti-Howard studies.

In our previous paper on the salt-finger case, we used modal truncation in the horizontal direction to argue for the importance of the horizontal flows for transport of heat and salt. Here we carry out an analogous study for the ODD case. In doing this, we would stress that a Galerkin study such as ours is an exploratory tool and is intended as a guide in rapidly identifying those features of the process that will then deserve a more realistic and resource-consuming study by means of direct numerical simulations.

Our plan is to present the equations in the next section in full and truncated form. We restrict our attention here to two-dimensional convection and refer to our previous paper (PS) for the derivation of the truncated equations, which at this stage are the same for both cases. In section 3 we explore the effect of a self-excited large-scale flow on ODD convection with parameters that allow for an analogy with the stellar case. Section 4 contains results for sea water parameters. We discuss and summarize our findings in section 5.

## 2 The equations

We consider the Boussinesq equations for double diffusion in two dimensions. The velocity field can be computed using a stream function $`\psi `$, for the horizontal and vertical velocities $`u=_z\psi `$ and $`w=_x\psi `$. The basic equations are:

$$\begin{array}{c}_t^2\psi J(\psi ,^2\psi )=\sigma R_T_xT+\sigma R_S_xS+\sigma ^4\psi \\ _tTJ(\psi ,T)=^2T\\ _tSJ(\psi ,S)=\frac{1}{\tau }^2S\end{array}$$
(1)
where the Rayleigh numbers are

$$R_T=\frac{g\alpha \mathrm{\Delta }Td^3}{\kappa _T\nu }\mathrm{and}R_S=\frac{g\beta \mathrm{\Delta }Sd^3}{\kappa _T\nu };$$
(2)
$`\alpha `$ and $`\beta `$ are respectively the thermal and saline expansion coefficients, $`g`$ is the acceleration of gravity, $`\nu `$ is the viscosity and $`\kappa _T`$ is the thermal diffusivity. We also use the Prandtl number $`\sigma =\nu /\kappa _T`$ and the Lewis number $`\tau =\kappa _T/\kappa _S`$, where $`\kappa _S`$ is the diffusivity of salt (or other solute). The equations have been made dimensionless with a unit of length, $`d`$, where the height of the fluid layer is $`\pi d`$, and a unit of time $`d^2/\kappa _T`$; the imposed temperature and salinity differences across the layer $`\mathrm{\Delta }T`$ and $`\mathrm{\Delta }S`$ serve as units for temperature and salt concentration. The signs of the Rayleigh numbers determine the particular kind of doubly-diffusive instability that triggers the motion. Here our interest is limited to the oscillatory case, which requires that the salinity difference across the fluid layer be stabilizing and that the temperature difference be destabilizing. In our notation, this is achieved when both Rayleigh numbers are negative.

In this study we consider only the case with an initially stable density stratification, that is $`R_T>R_S`$, where instability sets in through growing oscillations. (We recall also that oscillatory instability is possible for $`R_S>R_T`$, Baines & Gill, 1969.) As in PS, we use a modal version of these equations. We split the variables into a horizontally averaged part and a fluctuating part: $`\psi =\overline{\psi }(z,t)+\psi ^{}(x,z,t)`$; $`T=\overline{T}(z,t)+T^{}(x,z,t)`$; $`S=\overline{S}(z,t)+S^{}(x,z,t)`$. The fluctuating part is expanded as a Fourier sum along the $`x`$ direction, with amplitudes depending on $`(z,t)`$. The expansions are then plugged back into the equations (1), and truncated by retaining only the horizontally averaged quantities, and the amplitudes associated with a single horizontal wavenumber $`k_h`$.

Our model (modal) equations are these nine p.d.e.s:

$$\begin{array}{c}_tU=_z\left(\sigma _zU+\frac{1}{k_h}W^+_zW^{}\frac{1}{k_h}W^{}_zW^+\right)\\ _t\overline{T}=_z\left(_z\overline{T}W^+\mathrm{\Theta }^+W^{}\mathrm{\Theta }^{}\right)\\ _t\overline{S}=_z\left(\frac{1}{\tau }_z\overline{S}W^+\mathrm{\Phi }^+W^{}\mathrm{\Phi }^{}\right)\\ _t\mathrm{\Theta }^+=𝒟\mathrm{\Theta }^++k_hU\mathrm{\Theta }^{}W^+_z\overline{T}\\ _t\mathrm{\Theta }^{}=𝒟\mathrm{\Theta }^{}k_hU\mathrm{\Theta }^+W^{}_z\overline{T}\\ _t\mathrm{\Phi }^+=\frac{1}{\tau }𝒟\mathrm{\Phi }^++k_hU\mathrm{\Phi }^{}W^+_z\overline{S}\\ _t\mathrm{\Phi }^{}=\frac{1}{\tau }𝒟\mathrm{\Phi }^{}k_hU\mathrm{\Phi }^+W^{}_z\overline{S}\\ _t𝒟W^+=\sigma 𝒟^2W^+\sigma k_h^2\left(R_T\mathrm{\Theta }^+R_S\mathrm{\Phi }^+\right)k_hW^{}_{zz}U+k_hU𝒟W^{}\\ _t𝒟W^{}=\sigma 𝒟^2W^{}\sigma k_h^2\left(R_T\mathrm{\Theta }^{}R_S\mathrm{\Phi }^{}\right)+k_hW^+_{zz}Uk_hU𝒟W^+\end{array}$$
(3)
Here $`\mathrm{\Theta }^+`$, $`\mathrm{\Theta }^{}`$are the amplitudes associated with, respectively, the sine and the cosine modes of the temperature expansion. Analogously $`\mathrm{\Phi }^+`$, $`\mathrm{\Phi }^{}`$describe the salinity fluctuations. We use amplitudes $`W^+`$, $`W^{}`$ for the vertical velocity fluctuations, as well as the horizontal velocity $`U(z,t)`$, in place of the amplitudes of the modal expansion for $`\psi ^{}`$ and of $`\overline{\psi }`$. Finally, $`𝒟_{zz}k_h^2`$. The large scale flow represented by $`U`$ has often been neglected in previous studies of oscillatory convection. However, it is now experimentally well established that such a large scale, shearing flow does spontaneously appear in Rayleigh-Bénard convection (Siggia, 1994) where it changes the heat transport scaling law; we will show in the following the important role it may play in ODD convection.

## 3 ODD convection at low Prandtl and high Lewis number.

Conditions favorable to the onset of ODD convective instability are thought to occur in the core of a large variety of stars. In the simplest scenario, they are related to the existence of a central convection zone, in which hydrogen is slowly converted to helium. The mass fraction $`Y`$ of helium influences the differential buoyancy in the same way that the concentration of salt does in thermohaline convection. In modern evolutionary models, this central convection zone retreats, leaving behind a stabilizing radial helium gradient. In the case of massive stars, the entropy gradient in the same region is typically unstable and this leads to a form of ODD convection known in astrophysics as semiconvection because it is effective in transporting helium but not heat.

In stars, the thermal diffusivity is orders of magnitude greater than diffusivity of the helium. Typical values for the Prandtl and the Lewis numbers (Merryfield, 1995) are

$$\sigma 510^5,\tau 10^8.$$
(4)
While solving the equations for the stellar plasma in this parameter regime is a challenge for the numerical modeler, insight into the nature of the process may be gained by solving the modal equations (3) at low Prandtl number and high Lewis number. We would hope in this way to ascertain whether the overstable conditions can lead to a state of overturning convection. This is important in estimating the rate of mixing of the chemical species in the doubly-diffusive region. To get an insight on the nature of mixing in such conditions, we shall adopt here the values $`\sigma =10^2`$, and $`\tau =10^2`$ and present results for $`R_T=200`$, $`R_S=400`$. Those calculations have been performed using a resolution of 64 modes in $`z`$. We choose the horizontal wavenumber $`k_h=1`$, which is close to the wavenumber of maximum growth rate in the linear theory for that set of parameters. (For a discussion of the linear stability analysis of doubly diffusive convection appropriate for the boundary conditions that we use here, see Baines and Gill, 1969.) We have explored other parameter regimes than those described here and found a rich and complicated range of behaviors. Extremely long transients may precede the time-asymptotic states. We have tried to choose representative situations in the regime with low Prandtl number and high Lewis number. These display two basic dynamical states: time-dependence with intense mixing (quantified in section 3.2), and steady motion with no mixing. Both are present within the parameter set that we next discuss in detail.

### 3.1 The wave-shear mixing cycle

Equations (3) are solved with fixed temperature, fixed salinity and free slip boundary conditions, using the numerical code described in PS. Initially we set $`\overline{T}=\overline{S}=z`$; $`U=A\mathrm{cos}(z)`$; all the other variables are set equal to $`A\mathrm{sin}(z)`$.

When $`A`$ is small (we used $`A=10^4`$) the nonlinearities are initially negligible, and the initial perturbation evolves quickly into a good approximation to the fastest growing eigenfunction. As the nonlinear terms become important, the system settles to a steady solution, shown in figure 1. This solution is dominated by a large scale shearing velocity, sustained by a single convective cell, which, in our example, circulates in the clockwise direction. The reflected configuration with an anti-clockwise cell is also possible and is obtained by starting with a negative amplitude. The most remarkable feature of this solution is its steady, stable density inversion, associated with a saddle point in the stream function (cf. fig. 2). The formation of this structure, an incomplete wave roll-up, is preceded by a temporal development that begins in the linear stages of the growth with a pair of cells of opposite signs, as in ordinary convection. This motion induces a wave pattern in the temperature and salinity fields. When the amplitudes $`W^+`$ and $`W^{}`$ reach the order of unity, the onset of the large scale flow $`U`$ makes the convective cells more and more asymmetrical, advecting the temperature and salinity wave crests close to each other as they encroach on the weaker cell, and stretching them apart in the stronger one. Finally, the weaker cell disappears altogether, leaving a tongue of warm, salty fluid folded over a tongue of cold, fresh fluid, above a saddle point in the stream function, in a remarkable advective-diffusive equilibrium. Further simulations show that this is a stable solution.

We also performed a run where the absolute values of the Rayleigh numbers were slowly and continuously lowered at the rate of $`(1.510^6)`$ units per thermal time, while keeping their ratio constant. The steady solution survives up to $`R_T14.5`$ when it abruptly falls to the conductive state, supporting the supposition that it appears through a sub-critical bifurcation.

A completely different dynamics appears if one sets $`A=1`$ in the initial conditions. We then find a solution that proceeds in cycles of approximately $`t_c246`$ thermal times. In figure 3 we show the time evolution of $`U`$ at $`z=\pi `$, during one such cycle. For convenience we set $`t=0`$ just before the large-scale flow reaches its minimum amplitude. At this time, the stream function shows a large, anti-clockwise cell, placed beside a smaller region of weak clockwise circulation. Growing, wave-like patterns shape the temperature and the salinity fields as they are advected across the domain by the large-scale flow (fig. 4, first panel from top). In this case, however, the shear is a leading order effect, and slows the growth of the overstable oscillations.

Shortly thereafter, the asymmetry of the convective cells drives up the amplitude of the large-scale flow. The second panel from the top of figure 4 shows the fields at time $`t=24`$. The stream function is dominated by a single, anti-clockwise cell. All clockwise circulation is confined to two small regions close to the upper and the lower boundary. Gradients in the salinity field increase, and the wave pattern is skewed. On the other hand, because the cell’s turnover time is of the same order as the thermal diffusivity time, the temperature field is only slightly perturbed.

At a later time, ($`t=48`$, fig. 4, third panel from top), the large scale flow is close to its maximum intensity. All clockwise circulation has disappeared from the stream function. A sequence of wave roll-ups has made the salinity field homogeneous over some parts of the domain, and it has greatly increased the gradients elsewhere. Notice the similarity between the stream function at this time and the stream function of the steady solution. In both cases, the circulation has the same sign in the whole domain, and the convective cells are separated by saddle points. Here, however, the strength of convection is weaker than in the steady solution. The flow is too slow to advect the temperature field appreciably, and the salinity plumes are too heavy to be sustained by the flow in a steady equilibrium.

The sequence of wave roll-ups continues until the shear becomes too intense to allow for any further roll-up ( $`t=63`$, fig. 4, fourth panel from top). At that point, the high shear disrupts the convective cells and destroys the inhomogeneities in the temperature and salinity fields, bringing them close to the conductive solution.

For the remaining fraction of the whole cycle, the flow is dominated by an intense, $`x`$-independent shear which flows rightward in the upper half of the domain, and leftward in the lower half ($`t=163`$, fig. 4, fifth panel from top). As for the steady case, for each solution there is a symmetric one with the shear flowing in the opposite direction; the particular choice of initial conditions selects between the two possibilities. The temperature and the salinity fields closely resemble the conductive solution $`T=S=z`$. The amplitudes $`W^+`$, $`W^{}`$, $`\mathrm{\Theta }^+`$, $`\mathrm{\Theta }^{}`$, $`\mathrm{\Phi }^+`$ and $`\mathrm{\Phi }^{}`$ are very small (less than $`510^4`$ in our non-dimensional units) at any height $`z`$. The evolution of $`U`$ is essentially decoupled from the dynamics of the other variables, and it is reduced to a slow viscous decay. Finally, when the intensity of the shear has been sufficiently damped, the overstable oscillations restart their growth, and the cycle begins again.

The dynamics that we just described is summarized in figure 5, where we plot the value of the shearing velocity $`U`$ at the top boundary versus the value of $`W^+`$ in the middle of the slab.

With the parameters discussed here, the cycles are not perfectly periodic. In fact, longer simulations reveal that the cyclic behavior is only a long transient. After approximately 4500 thermal times the flow falls back to the steady solution. On the other hand, at lower Rayleigh numbers, such as $`R_T=100`$ and $`R_S=200`$, the wave-shear mixing cycle appears to be a perfectly periodic solution. Further exploration of the parameter space has shown that, either as a transient or as the asymptotic regime, the roll-up of overstable oscillations by a self-excited shear is a behavior characteristic of the low Prandtl number regime.

In order to bring out the key role of the large-scale flow, we have carried out a simulation where the variable $`U`$ was kept equal to zero at all times, which is equivalent to removing the first of the equations (3). Starting from initial conditions of small amplitude, the numerical solution undergoes an oscillatory instability in agreement with the predictions of the linear analysis. However, after about 100 thermal times the solution has converged to a steady state, shown in figure 6. In this solution, all the salinity gradients are concentrated at the boundaries. In the interior of the domain the salinity is almost constant. This brings the system to a state analogous to purely thermal convection without shear, with two symmetric cells in the stream function field, and two steady thermal plumes (one ascending, the other descending) in steady equilibrium between advection and diffusion.

### 3.2 A Lagrangian view of the wave-shear mixing.

The numerical simulations described above show that the fluid can alternate long periods of laminar motion with relatively shorter times of tumultuous overturning. To better investigate the details of these *mixing events* we have computed the trajectories of a sample of fluid particles. For this purpose, the variables $`U`$, $`W^+`$ and $`W^{}`$ are interpolated using cubic spines (taking care to use the correct boundary conditions in their spline representation), to recover the fluid velocity $`𝐮_i(x,z)`$ at the position of the $`i`$-th particle. The advection equation $`\dot{𝐱}_i=𝐮_i`$ is then integrated at each Eulerian time step with a second order Adams-Bashfort scheme to obtain the particles’ trajectories.

We release 2500 particles at time $`t=0`$ on a square array of 50 by 50 points, with side length of 0.5 non-dimensional units, and with the lower left corner at (0, $`\pi /3`$). As time goes on, the first part of each thermohaline oscillation stretches the ensemble of particles vertically, the large-scale flow slants it, and the second part of the oscillation folds it. The interplay of waves and shear produces a dynamics very similar to that of a Baker’s map, and it is pictured in the time sequence of figure 7 (for an introduction to the physics of mixing see Ottino (1989); a simple example of such mixing flow is given by Aref’s blinking vortex, Aref, 1984). The stretching and folding continues until the thermohaline oscillations are disrupted by the shearing velocity. By the time $`U`$ reaches its maximum amplitude, ($`t65`$), the distribution of tracers has become vertically homogeneous.

The trajectory of an individual Lagrangian tracer (cf. fig. 8 ) illustrates what happens on longer time scales. The mixing events are clearly singled out by the periods of wiggly motion, alternating with periods where the advection is due to only the large scale velocity $`U`$. Since the sign of the shear remains constant, when a Lagrangian tracer happens to be in the upper half of the computational domain, it travels rightward, while it travels leftward when in the lower half. Dispersion in the horizontal direction occurs because the $`z`$ coordinate of a particle after a mixing event has a sensitive dependence upon the particle’s position before the mixing event. Each particle switches in an unpredictable way between rightward and leftward motion, performing a random walk along the $`x`$-axis. Thus we expect the horizontal transport of particles to be described by a diffusion process. On a dimensional basis, we can express its diffusion coefficient as $`K=U_c^2t_c`$. We estimate $`U_c`$ as the time and space average of $`|U|`$ and $`t_c`$ as the duration of the viscous decay of the shearing velocity that separates a mixing event from the next. In our non-dimensional units we have $`U_c0.85`$, $`t_c180`$, which yields $`K130`$.

This qualitative description is confirmed by evaluating the single-particle dispersion $`\xi ^2`$

$$\xi ^2(t;𝐱_i(t_0))=\frac{1}{N}\underset{i=1}{\overset{N}{}}\left|𝐱_i(t)𝐱_i(t_0)\right|^2$$
(5)
where $`N`$ is the number of Lagrangian tracers and $`𝐱_i(t_o)`$ is the initial position of the $`i`$-th one. In the case of Brownian motion, $`\xi ^2`$ does not depend on the initial positions and it is related to the diffusion coefficient $`K`$ of the Brownian stochastic process by the formula $`\xi ^2(t)=2Kt`$ (Gardiner, 1996). In our case, since the Eulerian flow has very different characteristics along the $`x`$ and the $`z`$ direction, we apply the definition of single-particle dispersion separately to each direction.

The results are shown in figure 9. The overall trend of $`\xi ^2(t)`$ for horizontal displacements is well approximated by a linear fit, which gives an effective diffusion coefficient $`K=98.6\pm 1.0`$ in non-dimensional units, in good agreement with the above rough estimate. However, on short time scales, the behavior of $`\xi ^2(t)`$ is not perfectly linear. Initially, in the $`x`$-direction, dispersion grows slowly, because the tracers are all clustered together into a single, coherent convection cell. After the first vertical mixing event (marked by the arrow in figure 9), the horizontal dispersion grows faster with time. At this stage, the dispersion is approximately ballistic (that is $`\xi ^2t^2`$) until $`U`$ can be approximated by a constant. However, as the large scale flow $`U`$ decays, the dispersion slows down as well, until the next mixing event occurs.

Along the vertical direction, the single-particle dispersion grows from zero to its maximum value during the first occurrence of wave-rollups, which makes the tracers’ distribution homogeneous in the vertical. The vertical dispersion then remains approximately constant, with statistical fluctuations during the vertical mixing events.

To confirm that the flow of figure 4 produces authentic mixing, we look for evidence that particle paths are independent of each other by evaluating the pair dispersion $`\mathrm{\Xi }^2`$

$$\mathrm{\Xi }^2(t)=\frac{2}{N(N1)}\underset{i=1}{\overset{N}{}}\underset{j>i}{}\left|𝐱_i(t)𝐱_j(t)\right|^2.$$
(6)
If each particle follows an independent Brownian path, then $`\mathrm{\Xi }(t)=4Kt`$. On evaluating the diffusion coefficient along the $`x`$-direction by using pair dispersion, we find $`K=102.0\pm 1.0`$, in close agreement with the value obtained in the single particle analysis. This is strong evidence that the long-time transport properties of the flow in the horizontal direction are those of Brownian motion, and can thus be modeled by a suitable eddy diffusion coefficient.

We carried out a similar simulation for the case where the large-scale flow is not allowed. As expected when the Eulerian flow is steady, the tracers moved along periodic orbits. Even after 1000 thermal times the 2500 tracers, seeded at the same initial position as the tracers in the simulation described above, remain confined in a circular belt and do not disperse.

To complement the information given by the advection of individual Lagrangian tracers, and to give a further visual indication that no barriers to mixing exist in the flow of figure 4, we carried out some simulations in which we use the streamfunction generated by our model to advect a passive scalar field $`c`$, according to the equation

$$_tcJ(\psi ,c)=\frac{\kappa _c}{\kappa _T}^2c.$$
(7)
The small amount of diffusion $`\frac{\kappa _c}{\kappa _T}=0.001`$ is necessary in order to maintain the stability of the finite volume code used to integrate equation (7). Figure 10 shows a time sequence of the numerical solution, and illustrates the high degree of homogeneity reached by the passive scalar during a single wave-shear mixing cycle.

## 4 ODD convection with heat and salt.

Because conditions favoring the onset of ODD convection occur at several places in the Earth’s oceans we investigate the behavior of the model in the parameter range appropriate for sea water. For this, we adopt $`\sigma =10`$, and $`\tau =10^2`$, which are approximately the values of the Prandtl and the Lewis numbers of sea water. With these parameters, the ODD instability occurs when the density ratio $`R_\rho =R_T/R_S`$ approaches unity (Baines & Gill, 1969).

When $`R_T=200`$, $`R_\rho =0.99`$ the system converges to a limit cycle characterized by a pair of plumes which oscillate up and down as standing waves, with periods of $`0.61`$ thermal times. Figure 11 shows a snapshot of the stream function, temperature and salinity fields close to the extremum of one oscillation: the warm, salty plume on the right has just reached the maximum elevation, the cold, fresh plume on the left has reached its maximum depth, and they are both moving back toward the middle of the domain. The large scale flow is not excited, and initial perturbations on $`U`$ decay as $`e^{\sigma t}`$. No Lagrangian mixing occurs in this simple flow.

Searching for a self-excited large-scale flow, we increased the Rayleigh number up to $`R_T=1000`$, and we have set the density ratio to $`R_\rho =0.95`$. To accurately resolve the high salinity gradients that develop at the top and bottom boundaries, we increased the vertical resolution to 256 modes. The initial conditions are the same as in section 3.1. We have tried several values for the initial amplitude $`A`$ in the range $`[10^3,10]`$. The solution that we find exhibits extremely long cycles of approximately $`100`$ thermal times, independently of the value of $`A`$.

Initially this solution shows a pair of oscillating plumes analogous to those found at $`R_T=200`$. They move faster (the period of oscillation is 0.3 thermal times, but it shrinks as the amplitude of the oscillations grows), and lead to sharper gradients in the scalar fields. The large-scale flow is excited, but it is so weak that the associate skewness of the convective cells is imperceptible in the figures. Figure 12 shows the stream function, the temperature and the salinity fields for the simulation starting with $`A=0.01`$ at time $`t=4.3`$, when the salinity fluctuations reach their maximum amplitude. The oscillations, however, do not last indefinitely as in the previous case. They damp out in a few thermal times, and the motion in the fluid ceases almost completely.

By looking at the vertical profiles of temperature and salinity it is possible to find a hint of the mechanism that switches off the oscillations. Figure 13 shows the time evolution of $`_z\overline{T}`$ and $`_z\overline{S}`$ computed at $`z=\pi /2`$ (that is, in the middle of the slab). Initially $`_z\overline{T}=_z\overline{S}=1`$. As the doubly-diffusive oscillations grow and then damp out, the slope of the salinity profile becomes steeper than the slope of temperature, and one is left with $`_z\overline{S}<1`$, while the temperature slope promptly returns to the diffusive solution $`_z\overline{T}=1`$. This results in a decreased density ratio in the middle of the slab, which locally makes the fluid linearly stable. Figure 14 shows the horizontally averaged profiles of temperature and salinity at the time $`t=6`$. The density ratio at $`z=\pi /2`$ is now $`R_\rho =\frac{_z\overline{T}}{_z\overline{S}}R_{\rho _{\mathrm{initial}}}=0.831`$, which is less than the critical value $`R_{\rho _{\mathrm{critical}}}=0.9163\mathrm{}`$ required to have a linear instability for the chosen values of $`R_T`$, $`\sigma `$, and $`\tau `$. At this time all motion has become negligibly small. The only dynamic in the system is the slow diffusion of the salinity field, which relaxes to the conductive solution with a time scale which, in non dimensional units is $`O(\tau )`$. Finally, when $`t100`$, the growth of ODD instability resumes and the whole cycle starts over with an alternation of short-lived oscillatory convection and slow diffusive relaxation. This dynamics does not lead to mixing. A set of Lagrangian tracers seeded in the way described in section 3.2 remains closely packed during the whole duration of the oscillations.

While the salinity build-up stabilizes the bulk of the fluid, we observe in figure 14 that, near the boundaries, $`_z\overline{T}<_z\overline{S}`$, implying that the fluid there is unstably stratified. However, we do not observe direct convection being triggered at the boundaries. We assume that this is a consequence of the truncation of the Boussinesq equations (1): the range of horizontal wavenumbers which would be unstable is not represented in the truncated model (3). Numerical solutions of the untruncated equations are needed in order to show whether convection actually arises at the boundaries and in what ways it contributes to the establishment of the step-like profiles seen in laboratory experiments (Huppert & Linden, 1979).

In this simulation the large-scale flow $`U`$ does not seem to play a dominant role. The doubly-diffusive oscillations are suppressed by the salinity build-up in the center of the domain, and not by the shear induced by $`U`$, as in the low Prandtl number case. Indeed, a simulation with $`R_T=500`$ and $`R_\rho =0.95`$ shows the same qualitative dynamics as the simulations with $`R_T=1000`$, but there is no self-excitation of a large-scale flow.

However, large-scale flows are prevalent in geophysics where they need not be convectively driven. Such flows may play the kind of role in the convective process that the shears we have been discussing do. For this reason, we have investigated the interaction between ODD oscillations and an imposed shearing velocity, $`U=\stackrel{~}{U}+\omega _0z`$ in the equations (3), and solving for $`\stackrel{~}{U}`$. Here $`\omega _0`$ is a constant shear. In our spectral code we represent $`z`$ as a trigonometric sum in the interval $`[0,\pi ]`$, truncated at the Nyquist frequency. We keep $`R_T=1000`$, $`R_\rho =0.95`$, and we impose an external shear $`\omega _0=0.2`$. Once again, oscillating plumes develop as a result of the ODD instability. The presence of the external shear makes the convective cells sufficiently asymmetric that they now can contribute to the build-up of a large-scale flow, which can reach a peak intensity of $`30`$ non dimensional units close to the boundaries. A snapshot of the fields at time $`t=4.3`$ is shown in figure 15. As in the case with low Prandtl number, the large-scale flow induces wave roll-ups, which mix the fluid until the salinity build-up chokes the motion. The resulting Lagrangian chaos is illustrated in figure 16.

## 5 Discussion and conclusions

This exploratory study of sheared ODD convection has revealed some previously unknown phenomenology. The main finding is that the presence of a large-scale shearing velocity (be it self-excited or externally imposed) induces vigorous mixing. The mechanism is reminiscent of the mixing induced by the Baker’s map: ODD waves vertically stretch material lines, the shear slants them, so that they can subsequently be folded up. This engenders a rapid, effective homogenization of Lagrangian tracers along the vertical. Horizontal mixing is also enhanced, in a way which recalls the phenomenon of shear dispersion (see, e.g. Jones & Young, 1994). As first noticed by Taylor, in the presence of a shear, molecular diffusion is enhanced along the direction of the flow. In our model the wave-shear mixing plays the role of molecular diffusion.

Further work is needed to asses at which Rayleigh numbers the large-scale flow appears. Spontaneous generation of shear seems very elusive at high Prandtl number, as in the case of heat and salt in water. This is supported by the laboratory experiments of Stamp & Griffiths (1997) and Stamp et al. (1998). They observe the generation of large-scale flows in convecting regions above and below an interface having ODD-favorable stratification (therefore giving further indication that this is a common feature of convecting fluids). The flows have the same intensity and direction above and below the interface, so they do not cause any shear on the interface, but merely translate it, either leftward or rightward. In an annular geometry, in the reference frame comoving with the mean flow, they observe gross features in agreement with the results on step-like interfaces of Gough & Toomre (1982), who employ a truncated model similar to ours.

On the other hand, we find that large-scale shearing flow is easily generated at very low Prandtl number, which is of interest in astrophysical applications. Unfortunately, carrying out either laboratory experiments or reliable direct numerical simulations in this regime is challenging. Other mechanisms can lead to mixing, distinct from that illustrated in this work. In particular, the simulations of Merryfield (1995) and of Biello (2001) both show ODD oscillations breaking into small-scale turbulence, as in the earlier discussion of Stevenson (1979). A larger investigation is needed to decide which of these mechanisms may be the more significant in different contexts.

In all cases where we have found mixing-inducing motion, it lasts for a relatively short interval, and it is then followed by a longer interval of non-mixing dynamics. The physical cause of this intermittent behavior is different in the low Prandtl number regime and in the sea water regime. In the first case the stabilizing factor is the shearing velocity $`U`$; in the second case the fluid is made linearly stable by a salinity build-up away from the boundaries. The present study has suggested that the astrophysical situation is the one in which the convectively induced shears come into play. This points to the value of a study directed at low Prandtl number ODD convection and suggests that we next seek simplifications that may help us in that limit.

Acknowledgments We thank the GFD program for bringing us all together. E. A. Spiegel was supported by NSF Applied Mathematics under DMS-99-72022. S. Talon was supported by NSERC of Canada and by the Canada Research Chair in Stellar Astrophysics awarded to Gilles Fontaine. We are grateful to Barry Ruddick for his comments which contributed to improve this manuscript and to Bill Young who helped us clarifying the issue of what defines mixing.