text:
# Spin-Induced Disk Precession in the Supermassive Black Hole at the Galactic Center

## 1 Introduction

Sgr A\*’s radio emission is known to be variable, with a fluctuation amplitude that increases toward high frequencies (Zhao et al. 1993). At mm wavelengths, large-amplitude ($`100\%`$) variations have also been observed (Wright & Backer 1993; Miyazaki et al. 1999). Its time-averaged spectrum is roughly a power law below $`100`$ GHz, with a flux density $`S_\nu \nu ^\alpha `$, and $`\alpha 0.190.34`$. In the mm/sub-mm region, however, Sgr A\*’s radiative output is dominated by a “bump” extending above this power law (Zylka et al. 1992; Falcke et al. 1998).

Unlike the high level of polarization seen in this object at mm/sub-mm wavelengths (Aitken et al. 2000), Sgr A\* reveals a lack ($`<1\%`$) of linear polarization below $`112`$ GHz, though some circular polarization ($`1\%`$) has been detected (Bower et al. 1999; Bower et al. 2001). These prominent spectral and polarimetric differences (Melia, Bromley, & Liu 2001) between the cm and the mm/sub-mm bands suggest two different emission components in Sgr A\*. Because higher frequencies correspond to smaller spatial scales (Melia 1992; Melia et al. 1992; Narayan et al. 1995), the mm/sub-mm radiation is likely produced in the vicinity of the black hole. Earlier work (e.g., Melia 1992, 1994; Coker & Melia 1997) has indicated that Sgr A\* is accreting from the stellar winds surrounding the black hole and that the infalling gas circularizes at a radius of $`20800r_S`$. We recently (Melia, Liu, & Coker 2001; Bromley et al. 2001) demonstrated that the inner $`10r_S`$ of the resultant Keplerian structure can not only account for the mm/sub-mm properties via thermal synchrotron emission, but it may also produce Sgr A\*’s X-ray spectrum in the quiescent state (Baganoff et al. 2001) via Comptonization of the mm/sub-mm photons. On the other hand, the cm radio emission appears to be produced by non-thermal synchrotron emission (Liu & Melia 2001). Meanwhile, the recently detected X-ray flare (Baganoff et al. 2001) from Sgr A\* appears to be produced by thermal bremsstrahlung radiation resulting from a transient injection of additional mass flowing through the disk (Liu & Melia 2002).

The 106-day cycle seen at $`1.3`$ and $`2.0`$ cm (Zhao et al. 2001) is intriguing because high resolution VLA observations have already ruled out the possibility that such a period might be produced by an orbiting emitting object, which would have been distinguished spatially from the main source (Bower et al. 1998). This inference is supported by the observed lack of proper motion (Backer et al. 1999; Reid et al. 1999) of Sgr A\*, which precludes its possible association with rapidly moving components. In addition, any stellar origin for such a source falls short of the power required to account for the measured radio emission. It is reasonable to conclude, therefore, that the periodic variations in Sgr A\* are intrinsic to the source. In this Letter, we examine the possibility that the observed 106-day radio modulation is associated with a precession of the disk in a Kerr metric and explore the implications this may have on the nature of Sgr A\*.

## 2 Nature of the Long-Term Modulation

The characteristics of this 106-day cycle constrain the nature of its origin rather tightly (Zhao et al. 2001). First, the observed period is independent of the wavelength, and recent Submillimeter Array (SMA) and VLA observations of Sgr A\* show that the mm and cm variations are correlated (Zhao et al. 2001a). Thus, the fact that the emission at different frequencies is produced on different spatial scales (Melia et al. 1992), suggests that this period should be induced by a single process, which can cause correlated fluctuations across a broad range of wavelengths. Second, the period is four orders of magnitude longer than the dynamical time scale at the marginally stable orbit $`r_{ms}`$. So it can either be produced on a much larger spatial scale, or it may be associated with an intrinsic property of the supermassive black hole itself, such as a nonzero spin. However, recent VLA observations have indicated that the 2 cm emission is produced within a region no larger than $`140r_S`$ (Bower et al. 2002), for which the corresponding dynamical time scale is about one and a half days, which is much smaller than the observed period. Higher frequency emission is expected to be produced within still smaller regions, associated with even smaller time scales.

We may ask then, whether this modulation could be produced by a corrugation wave in an accretion disk, which is used to account for the quasi-periodic oscillations (QPOs) seen in low-mass X-ray binaries. These waves have periods that are much longer than the corresponding dynamical time scale (on a comparable spatial scale; Kota 1990), but the periods depend on radius and thus may not be able to account for the first feature described above. Moreover, the observed light curves show quite stable periodic fluctuations, rather than the uncorrelated segments constituting QPOs. We thus conclude that the period must be associated with a single process evolving on a small spatial scale (i.e., certainly less than $`100r_S`$ in size).

A schematic diagram of the assumed geometry for this system is shown in Figure 1. The oscillation frequency of a perturbation in the perpendicular direction to a circular orbit in the equatorial plane of a Kerr black hole has been derived by Kato (1990), and may be written in the form $`\nu _\theta ^2=\nu _\varphi ^2(14aM^{1/2}r^{3/2}+3a^2r^2),`$ where $`G=c=1`$, and $`\nu _\varphi =\pm M^{1/2}r^{3/2}\left[2\pi (1\pm aM^{1/2}r^{3/2})\right]^1`$ is the frequency of a circular orbit at radius $`r`$ measured by a static observer at infinity (Bardeen et al. 1972). For small values of $`a/M`$, the nodal precession frequency is then $`\nu _{p\theta }\nu _\varphi \nu _\theta aMr^3/\pi .`$ Because this frequency depends on $`r`$, the final configuration of a disk that is misaligned with the black hole spin at large radii depends on whether the viscous time scale $`t_{vis}`$ is shorter or longer than the precession period $`\nu _{p\theta }^{}{}_{}{}^{1}`$. For a thin cold disk, the differential Lense-Thirring precession will dominate the internal coupling of the plasma at small radii and will therefore lead to the so-called Bardeen-Petterson effect (Bardeen & Petterson, 1975), in which the inner region flattens toward the equatorial plane, producing a warped accretion pattern. As shown by Nelson & Papaloizou (2000), however, thicker disks with a mid-plane Mach number of $`5`$ or less can suppress warping effects due to the coupling provided by pressure gradients in the gas. The mid-plane Mach number in Sgr A\* is $`3`$ (Melia, Liu, & Coker 2001), so it appears that the disk in this system precesses more or less as a rigid body.

Since the surface density of the disk is roughly a constant (Melia, Liu, & Coker 2001), the total precessional torque acting on it is

$$T=\mathrm{sin}\xi _{r_i}^{r_o}2\pi \nu _{p\theta }v_kr2\pi r\mathrm{\Sigma }𝑑r=8\pi \mathrm{sin}\xi aM(GM)^{1/2}\mathrm{\Sigma }\left[1(r_i/r_o)^{1/2}\right]/r_i^{1/2},$$
(1)
where $`\xi `$ is the angle between the axis of the disk and the black hole spin, $`\mathrm{\Sigma }`$ is the surface density, $`v_k`$ is the Keplerian velocity and $`r_i`$ and $`r_o`$ are, respectively, the inner and outer radii. The disk’s total angular momentum is given by $`L=0.8\pi (GM)^{1/2}\mathrm{\Sigma }r_o^{5/2}\left[1(r_i/r_o)^{5/2}\right].`$ Thus, the precession period of the disk is

$$P=2\pi \mathrm{sin}\xi (L/T)=\frac{\pi r_o^{5/2}r_i^{1/2}\left[1(r_i/r_o)^{5/2}\right]}{5aM\left[1(r_i/r_o)^{1/2}\right]}.$$
(2)
In these calculations we have neglected the higher order general relativistic effects, which will introduce corrections of order $`r_S/r`$. The gas will fall into the black hole very quickly once it crosses $`r_{ms}`$. The inner boundary $`r_i`$ of the accretion disk may therefore be taken to be $`r_{ms}`$, which is $`3r_S`$ for a black hole with small spin. Corrections introduced by the higher order general relativistic effects should then be no bigger than $`30\%`$.

The disk’s outer radius, on the other hand, is not well constrained, even though several observational lines of evidence point to a compact, hot, magnetized, Keplerian region as the dominant radiator (Melia, Liu, & Coker 2000; 2001; Liu & Melia 2001). Earlier hydrodynamical simulations, however, do suggest that the accretion flow in Sgr A\* circularizes at a radius of $`20800r_S`$ (Coker & Melia 1997). In the following discussion, we will take $`r_o=30r_S`$ as a fiducial outer boundary. It is straightforward to see that an identification of this precession period with the detected radio modulation then implies a black hole spin

$$a/M=0.088(r_i/3r_S)^{1/2}(r_o/30r_S)^{5/2}\frac{1(10r_i/r_o)^{5/2}}{1(10r_i/r_o)^{1/2}}.$$
(3)
(Note that in these units, $`r_S=2M`$.)

There are several questions that arise from this discussion. The first concerns whether or not the accreted angular momentum vector can disrupt the small inner disk as it precesses around the black hole’s spin axis, since $`t_{vis}P`$ (Liu, & Melia 2002). It should be emphasized, however, that the viscous time scale is the time required for the gas to flow through the disk, while the time required to realign the disk via accretion is determined by the disk’s total angular momentum vector and the angular momentum flux through it. This depends primarily on whether or not angular momentum is transported efficiently through the disk. The infalling gas carries angular momentum inward. However, a magnetic stress applied to the flow by the accretor can introduce an outward angular mementum flux that can effectively cancel much of the inward flux. Recent work on processes below $`r_{ms}`$ has shown that the torque at the inner edge of the disk is not necessarily zero, and that angular momentum can therefore be extracted from the magnetized gas below $`r_{ms}`$ (see, e.g., Krolik 1999; Gammie, 1999). Consequently, the angular momentum flux through the disk may be sufficiently small (perhaps even zero) that its effect on the inner structure may be negligible. However, the disk can be disrupted by significant fluctuations in the angular momentum induced by variations in the large scale accretion flow. As we alluded to earlier, this type of event is expected to occur roughly every few hundred years (Coker & Melia 1997).

Second, how does a precessing disk affect the cm emission, given that earlier work has already established the implausibility of an accretion flow producing this radiation via a thermal synchrotron process (Liu & Melia 2001)? One possibility is that the $`120`$ cm spectrum of Sgr A\* is produced by shock-accelerated non-thermal particles in the region outside $`30r_S`$, where the infalling gas is circularizing. In that case, the radio modulation could be the result of optical depth effects as the flattened distribution of particles precesses about the black hole’s spin axis. One should note, however, that due to the weak internal coupling in the quasi-spherical region (Liu, & Melia 2001), the precessing disk should not affect this large scale accretion. A more interesting possibility is suggested by the fact that the integrated $`120`$ cm luminosity of Sgr A\* would be comparable to the power extracted from its spin energy via a Blandford-Znajek type of electromagnetic process if $`a/M0.088`$. The 106-day modulation would then presumably arise when the precessing disk periodically shadows the non-thermal particles flooding the region surrounding the black hole as they escape from their creation site near the event horizon.

For this mechanism to work, however, the vacuum must break down as $`r`$ approaches $`r_S`$, but one can see right away that if $`a/M`$ is indeed as low as the value inferred from the radio luminosity of Sgr A\*, this break down cannot happen with $`\gamma `$-rays produced via curvature radiation (Blandford & Znajek 1977). Instead, the prevalence of mm/sub-mm photons offers the possibility that inverse Compton scatterings may produce a sufficient number of energetic $`\gamma `$-rays that can materialize to form the required sea of electron-positron pairs. To quantify this thought, we note that Sgr A\*’s sub-mm spectrum turns over at $`\nu _{mm}236`$ GHz with a flux density of $`2.8`$ Jy (Falcke et al. 1998). The corresponding radiation energy density near the event horizon is therefore $`1.02`$ erg cm-3 (assuming a distance to the Galactic Center of $`8.5`$ kpc) which is larger than the critical radiation energy density of $`0.57`$ ergs cm-3 required to initiate the cascade of electron-positron pair creation (see Eq. 2.9 of Blandford & Znajek 1977). To break down the vacuum, this photon number density must be matched to an adequately intense magnetic field: $`B29(0.088M/a)(2.6\times 10^6M_{}/M)(236\mathrm{G}\mathrm{H}\mathrm{z}/\nu _{mm})\text{gauss}.`$ By comparison, the magnetic field in the accretion disk is inferred to be about $`20`$ gauss (Melia, Liu, & Coker 2001). However, Krolik (1999) has argued that below $`r_{ms}`$, the magnetic field energy density will increase inward and should saturate to the local rest mass energy density of the infalling gas. Gammie (1999) has built a specific model where the magnetic field energy density increases significantly inside $`r_{ms}`$. Based on these estimates, we would therefore expect the magnetic field to be even larger than $`29`$ gauss near the event horizon. Under these conditions, the power extracted from the spinning black hole would scale in the following fashion: $`L1.1\times 10^{34}(B/29\mathrm{gauss})^2(a/0.088M)^2(M/2.6\times 10^6M_{})^2\text{ergs s}^1,`$ a value that is based on the integrated $`120`$ cm radio emissivity of Sgr A\*.

Unfortunately, the actual process by which the energy is transferred from the black hole to the particles is not known with sufficient certainty for us to predict with confidence the characteristics of their radiative output, other than their maximum available power. Nonetheless, by assuming an equipartition between their energy density and that of the magnetic field $`B`$, one can show that the angular size $`\theta `$ and the value of $`B`$ for this non-thermal source may be inferred from its observed turnover frequency $`\nu _{mm}`$ and the corresponding flux density $`F_{mm}`$ (Marscher 1983): $`\theta =\left[8\pi \times 10^{10}n(\alpha )/b(\alpha )^2(2\alpha 1)(\gamma _1m_ec^2)^{2\alpha 1}D_{\mathrm{Gpc}}\right]^{1/(4\alpha +15)}\nu _{mm}^1F_{mm}^{(2\alpha +7)/(4\alpha +15)}`$ mas, $`B=10^5b(\alpha )\left[8\pi \times 10^{10}n(\alpha )/b(\alpha )^2(2\alpha 1)(\gamma _1m_ec^2)^{2\alpha 1}D_{\mathrm{Gpc}}\right]^{4/(4\alpha +15)}\nu _{mm}F_{mm}^{2/(4\alpha +15)}`$ gauss, where $`\nu _{mm}`$ and $`F_{mm}`$ are in units of GHz and Jy respectively, $`\alpha `$ is the spectral index of the power-law synchrotron emission in the optically thin region, $`n(\alpha )`$ and $`b(\alpha )`$ are parameters given in Marscher (1983), $`\gamma _1m_ec^2`$ is the low-energy cutoff of the electron distribution, and $`D_{\mathrm{Gpc}}`$ is the distance to the source in Gpc. The results are not very sensitive to the spectral index $`\alpha `$. For Sgr A\*, adopting an $`\alpha `$ of $`1`$, we find that the $`22`$ cm radio emission ($`0.53`$ Jy) is produced within a region whose projected size is $`4.7`$ mas, corresponding to $`760r_S`$. Its magnetic field intensity is $`0.32`$ gauss. At $`1.3`$ cm (with a flux of $`1.1`$ Jy), where the 106-day period was detected, the corresponding source size is $`65r_S`$. We therefore conclude that a precessing disk with a diameter of $`60r_S`$ can shadow a significant fraction of the emitting region at $`1.3`$ cm, and thereby produce a noticeable modulation at this wavelength, but not at $`22`$ cm. This expected behavior is consistent with what was observed, in that the amplitude of the modulation was greatest at $`1.3`$ cm, dropping monotonically toward longer wavelengths, and finally falling below statistically significant levels beyond $`6`$ cm.

## 3 Discussion and Conclusions

If this picture is correct, there are several immediate consequences for the nature of Sgr A\*’s spectrum at mm/sub-mm and X-ray wavelengths. If we assume that the non-thermal emission from this process extends into the sub-mm region with a spectral index of $`0.2`$ and it turns over at $`236`$ GHz, then the flux density at the turnover is $`1.5`$ Jy. The calculated source size and magnetic field intensity at that frequency are then $`7.4r_S`$ and $`49`$ gauss, respectively, which are also consistent with the magnetic field inferred from our accretion model. But the fact that the observed $`236`$ GHz flux density ($`2.8`$ Jy) is larger than that produced by the non-thermal particles then points to the need for a hot accretion disk and its contribution to the mm/sub-mm bump.

Thinking about the spectrum produced by the non-thermal particles in the region beyond the observed break at $`\nu _{mm}236`$ GHz, we note that extending the best fit X-ray spectral index (Baganoff et al. 2001) of $`1.2`$ backward toward the mm/sub-mm region predicts a flux density of $`1.4`$ Jy at $`\nu _{mm}`$, which is very close to the value one obtains by extrapolating the non-thermal radio emission toward this frequency (Liu & Melia 2001). Although this may itself be a coincidence, the convergence of Sgr A\*’s low- and high-energy spectral components at $`236`$ GHz suggests that this may therefore be the turnover frequency that separates the optically thick and thin non-thermal emission components in Sgr A\*. The Chandra-detected X-rays in the quiescent state would then be due to optically-thin non-thermal synchrotron processes. In this instance, however, $`a/M`$ would need to be somewhat larger than $`0.088`$ in order to accommodate the implied non-thermal particle luminosity of $`2.9\times 10^{34}`$ ergs s-1.

There is clearly much to be done with this picture. Not coincidentally, the physics of accretion below the marginally stable orbit has attracted some attention in recent years (Krolik 1999; Gammie 1999; Agol & Krolik 2000). It appears that a nonzero stress at this radius can change the properties of the accretion disk significantly, and strong dissipation may be possible inside this orbit. Since the sharp increase in radial velocity toward the event horizon leads to a significantly lower number density, a continued equipartition between the particles and the magnetic field may produce very high-energy particles via magnetic dissipation. Under some circumstances, this combination of low particle density and magnetic reconnection may permit some particles to escape and form a directed, relativistic outflow, a scenario that provides an alternative to the trapped non-thermal particle model we have been describing here. Further study of these various possibilities is warranted.

We note, finally, that the shadowing effect from the precessing disk should be even more pronounced at $`43`$ GHz than at $`22`$ GHz, based on our estimate for the relative sizes of the emitting regions, suggesting that an observation of a 106-day period at this higher frequency should reveal an even stronger modulation than has been seen thus far at the longer wavelengths.

We are very grateful to the anonymous referee, whose comments have led to a significant improvement of the manuscript. This research was partially supported by NASA under grants NAG5-8239 and NAG5-9205, and has made use of NASA’s Astrophysics Data System Abstract Service. FM is grateful to the University of Melbourne for its support (through a Miegunyah Fellowship).


format: text

text:
# Random inelasticity and velocity fluctuations in a driven granular gas

## 1 Introduction

Whereas equilibrium statistical mechanics has reached a rather mature phase, the understanding of non-equilibrium processes is far from complete. In particular, granular (and thus inelastic) gases Jaeger  driven into a non-equilibrium steady state by a suitable injection of energy define a stimulating research field where theoretical predictions can be confronted against model experiments, with the aim to understand the possible deviations from equilibrium behavior. A good probe to quantify these deviations is the velocity distribution of the grains, $`P(v)`$, which has focused sustained attention recently, and has been shown to exhibit pronounced differences from Maxwell-Boltzmann statistics Olafsen ; Losert ; Rouyer ; Kudrolli ; Aranson . Several authors reported a stretched exponential law \[on the whole range of velocities available, which covers an accuracy of $`5`$ to $`6`$ orders of magnitude for $`P(v)`$\]

$$P(v)\mathrm{exp}[(v/v_0)^\nu ],$$
(1)
with an exponent $`\nu `$ close to $`3/2`$ Losert ; Rouyer ; Aranson  (here $`v_0`$ is the “thermal” r.m.s. velocity). This behavior was observed for the horizontal velocity components of a vertically vibrated 2D system of steel beads in a wide range of driving frequencies and densities Rouyer , but also in a three dimensional electrostatically driven granular gas Aranson .

At this point, some questions naturally arise, that will be addressed below: (i) is it possible to find a consistent model where this velocity distribution would emerge? (ii) what physical ingredients are required?

One possible approach consists in performing “realistic” molecular dynamics simulations. The model of inelastic hard spheres (IHS) with binary momentum conserving collisions, and a “reasonable” restitution coefficient Jaeger  provides the simplest candidate. The energy loss in a collision is proportional to the inelasticity parameter $`1\alpha _0^2`$ where $`\alpha _0`$ is the coefficient of normal restitution ($`0<\alpha _01`$), which in the simplest and efficient approximation is a constant independent on the relative velocity of colliding partners. Such an approach has been presented in vibrated ; Brey , and allows to reproduce the experimental velocity statistics with a good accuracy. The possible lack of universality has also been addressed in Brey .

Another route, which contrary to the previous numerical one has the merit to allow an analytical derivation of $`\nu `$ in some cases twan , consists in formulating an effective modeling of the energy injection, considering idealized homogeneous systems of inelastic hard spheres (given the experiments reported in Rouyer , the assumption of homogeneity is well founded, see below). From this point of view, a simple and popular model consists in IHS with constant inelasticity, with a homogeneous forcing described by a “stochastic thermostat” Williams ; Puglisi ; twan ; Pre ; Montanero ; Moon ; PRE2 ; Garzo . This model has attracted attention, in particular because it has been shown analytically twan  that $`P(v)`$ exhibits a high energy tail of the form of Eq. (1) with $`\nu =3/2`$, independent on dimension and restitution coefficient, in apparent agreement with the experiments. This result holds at the mean-field level for an homogeneous system. The above model, where an external white-noise driving force acts on the particles and thus injects energy through random “kicks” between the collisions, is therefore considered to provide a relevant theoretical framework to quantify the non Gaussian character of velocity distributions.

However, we shall see below that this uniformly heated model –in its simplest version– is unable to reproduce the experimental data: if simulated in dimension 2 or higher with an experimentally relevant value of the restitution coefficient ($`\alpha _0`$ between $`0.7`$ and $`1`$), the obtained distribution $`P(v)`$ is indistinguishable from a Gaussian within the experimental accuracy 111 if this model is simulated in dimension $`1`$ with $`\alpha _0[0.7;1]`$, the obtained non Gaussian behaviour at thermal velocities corresponds to $`\nu >2`$ instead of $`\nu 3/2`$ in the experiments.; in fact, the range of velocities for which the high energy behavior $`\mathrm{exp}[(v/v_0)^{3/2}]`$ may be observed is even beyond reach of precise numerical procedures and corresponds to a regime where $`P(v)`$ is practically vanishing \[lower than $`10^6P(0)`$\]. The predictions of this model are consequently incompatible with the experimental velocity distributions, that show important non Gaussian features already at thermal velocity scale Losert ; Rouyer ; Aranson . We emphasize however that generalizations of the aforementioned heated model have been proposed Cafiero\_prl ; Cafiero . With a convenient choice of the extra parameters introduced, one may obtain a velocity distribution close to that measured in the experiments (we will come back to this point in section 4). However, in such approaches, energy is injected in between the collisions whereas in the experiments we are interested in, the transfer of horizontal momentum takes place at every inter-particle collision (see section 3). We will therefore focus on this feature for the 2D experiment reported in Rouyer  and investigate in details the collision dynamics in the horizontal direction (section 2). The vibrated system under study there shows important density and granular temperature gradients, especially close to the boundaries which inject energy, but since the shaking is violent, there is a region where both gradients are very small simultaneously. The velocity acquisition in Rouyer  has been restricted to this region, where the system, although open, may be considered as homogeneous. In this article, we thus consider the following question : remaining at the level of a homogeneous system, what ingredients are required for a self-consistent effective description of the horizontal degrees of freedom, that exhibit the stretched exponential law (1) ?

## 2 Effective restitution coefficients

In order to characterize the collision process in the horizontally projected system, we have measured directly the effective 1D restitution coefficient from the experimental data provided by K. Feitosa and N. Menon, for a gas of stainless steel spheres (the system investigated in Rouyer ) but also for glass, brass and aluminum beads Feitosa , which allow to sample a wide range of nominal inelasticities. Let us recall briefly the experimental set-up. The balls (diameter: $`d=1.600\pm 0.002mm`$) are confined to a vertical, rectangular cage (32 $`d`$ high x 48 $`d`$ wide x 1.1 $`d`$ thick) sandwiched between two parallel plates of Plexiglas. The cage is vibrated vertically at a frequency of $`60Hz`$ and amplitudes up to $`2.4d`$, producing maximum accelerations, $`\mathrm{\Gamma }`$, and velocities, $`v_0`$, of $`56g`$ and $`1.45m/s`$ respectively. The motions of the balls are recorded with a high-speed camera which allows a location of each ball with a precision of 0.03 $`d`$. The results we discuss here are taken in a rectangular (10 $`d`$ x 20 $`d`$) window around the geometrical center of the cell, where, as mentioned above, density and granular temperature are almost homogeneous Rouyer . Moreover, the measured velocity distributions do not vary with height nor with the phase of the vibration cycle. The experimental data can thus be considered as obtained in the bulk of a two-dimensional homogeneous (but open) system, reasonably far from the boundaries.

The horizontal component of relative velocities are computed before ($`g_x`$) and after ($`g_x^{}`$) each collision, from which we deduce the effective restitution coefficient

$$\alpha _{1d}=\frac{|g_x^{}|}{|g_x|}.$$
(2)
Fig. 1 displays the histogram $`\mu _{exp}(\alpha _{1d})`$ obtained from the experimental data for different materials. At large $`\alpha _{1d}`$, a power-law tail is evidenced. Note that values $`\alpha _{1d}>1`$ are expected, due to the transfer from vertical to horizontal translational kinetic energy Barrat\_Trizac .

The strong correlations between relative horizontal velocity $`g_x`$ and $`\alpha _{1d}`$ are clearly seen in the scatter-plot (inset of Fig. 1), with a very sharp cutoff above the second bisector $`\alpha _{1d}1/g_x`$. This cut-off follows from the definition of $`\alpha _{1d}`$: since the post-collision velocity is finite, large values of $`\alpha _{1d}`$ may only be obtained for small values of $`g_x`$, in which case Eq. (2) implies that the maximum $`\alpha _{1d}`$ is of order $`1/g_x`$.

These features are qualitatively the same for all materials investigated. Moreover, three different densities have been investigated for steel beads, and the same distributions have been obtained Feitosa . Similar distributions have also been measured in a different experimental set-up with rolling beads blair . The relatively small number of collisions investigated does not however allow us to get accurate histograms for the joint distributions $`\mu (\alpha _{1d},g_x)`$ nor the conditional $`\mu (\alpha _{1d}|g_x)`$. The correlations evidenced in Fig. 1 nevertheless play a crucial role, as will be shown below.

More insight into the conditional distributions $`\mu (\alpha _{1d}|g_x)`$ has been obtained in molecular dynamics of two-dimensional IHS driven by vibrating walls in Ref. vibrated . Histograms of $`\mu (\alpha _{1d})`$ similar to the experimental results have been obtained. It turns out that the global $`\mu (\alpha _{1d})`$ is almost insensitive to the details of the system (density, velocity of the vibrating walls…), while this dependence exists for $`\mu (\alpha _{1d}|g_x)`$ (and also for $`P(v)`$). The following characteristics of $`\mu (\alpha _{1d}|g_x)`$ have been obtained: at constant $`g_x`$, $`\mu (\alpha _{1d}|g_x)`$ is almost constant for $`\alpha _{1d}[0,\alpha _0]`$, has a small peak at $`\alpha _0`$, and decreases as $`\mathrm{exp}(A(g_x)\alpha _{1d}^2)`$ for $`\alpha _{1d}>\alpha _0`$, with $`A(g_x)g_x^2`$.

## 3 Random Restitution Coefficient model

We have considered the possibility to mimic the experimental distributions by including randomness in the restitution coefficient, following the approach of Barrat\_Trizac  where a random restitution coefficient (RRC) model was introduced to account for the fact that, in vertically shaken granular gases, the energy is transferred to the vertical degrees of freedom by the moving piston, and then to the horizontal ones through grain/grain collisions only. The heating of horizontal degrees of freedom thus occurs through the inter-particle collisions, and not in between as in the “stochastic thermostat” approach. Moreover, a globally dissipative collision may correspond to an energy gain for the horizontal components of the velocities. This leads to the study of horizontally projected collisions with an effective restitution coefficient that can be either smaller or larger than 1, as the experimental data of Fig. 1 indeed show. In our case, the RRC model is therefore an effective approach in $`1`$ dimension (since the original collisions are two-dimensional), in which IHS undergo binary, momentum-conserving collisions with a restitution coefficient $`\alpha _{1d}`$ drawn randomly at each collision from a given distribution $`\mu (\alpha _{1d})`$ which should mimic $`\mu _{exp}(\alpha _{1d})`$. Even if the original collisions are not random, the projected ones may be considered as such.

Note that the energy injection is here given solely by the values of $`\alpha _{1d}`$ larger than $`1`$. The simulations are performed at the mean-field level of the homogeneous non-linear Boltzmann equation for point particles, solved by the Direct Simulation Monte Carlo method (DSMC) Bird . The velocity statistics is computed in the non-equilibrium steady state which is reached after a transient.

We have first considered distributions with a large tail in order to reproduce the experimental $`\mu (\alpha _{1d})`$, but without any correlations with the relative velocities of the colliding particles; in this case, it turns out that $`P(v)`$ has a power-law decay at large $`v`$, in marked contrast with experimental results 222It is noteworthy that a large $`v`$ power law is actually obtained for any distribution $`\mu (\alpha _{1d})`$, see Barrat\_Trizac .. This approach consequently needs to be refined and the next crucial step is to take into account the correlations between $`\alpha _{1d}`$ and $`g_x`$, with the insight given by the experimental scatter-plot (inset of Fig. 1) and by the molecular dynamics results vibrated .

We have thus used distributions decaying as $`\mu (\alpha _{1d}|\stackrel{~}{g}_x)\mathrm{exp}[(\alpha _{1d}\stackrel{~}{g}_x)^2/R]`$ at large $`\alpha _{1d}`$, where $`\stackrel{~}{g}_x=g_x/v_0`$ is the rescaled velocity defined from the total kinetic energy of the system ($`v_0^2=v^2`$), and the parameter $`R`$ can be varied with values of order 1. The function $`\mu (\alpha _{1d}|g_x)`$ is the only input needed to simulate the RRC model. For the consistency of the approach, the distribution $`\mu (\alpha _{1d})`$ measured in the simulation needs to be close to its experimental counterpart. This comparison is displayed in Fig. 1 and justifies a posteriori the choice made for $`\mu (\alpha _{1d}|g_x)`$. Both experimental and numerical distributions $`\mu (\alpha _{1d})`$ display a power law tail of the form $`\alpha _{1d}^n`$ with $`n3`$.

The velocity distribution obtained from the RRC model is compared to the experimental measure in Fig. 2. The agreement is satisfactory over the whole range of velocities; in particular, the RRC distributions is compatible with the stretched exponential behavior reported in Rouyer , with an exponent $`\nu `$ close to 1.5. Note that, since no precise experimental data is available for the distributions of restitution parameters conditioned by relative precollisional velocity, the parameters $`R`$ are tunable \[as long as the global $`\mu (\alpha _{1d})`$ coincides with the experimental one\], and the value giving the best agreement has been chosen ($`2R4`$). The agreement remains satisfactory upon changing $`R`$, provided that the resulting large $`\alpha `$ cutoff remains sharp (i.e. $`R`$ should not be too large).

The RRC model therefore provides a self-consistent framework which allows to reproduce the experimental $`P(v)`$ if implemented with the correct distribution of effective coefficients. Moreover, the velocity statistics depends on the distribution of effective restitution coefficients: a broader $`\mu (\alpha _{1d}|g_x)`$ leads to a broader $`P(v)`$, consistently with the numerical study of vibrated  which showed both broader $`P(v)`$ and $`\mu (\alpha _{1d}|g_x)`$ as e.g. the density is increased. Both distributions $`P`$ and $`\mu `$ are equally sensitive to a possible non universality (dependence on material properties). As a consequence, an accurate experimental measure of $`\mu `$ appears as complementary to the direct computation of $`P(v)`$, in order to assess the experimentally difficult question of the velocity statistics universality.

## 4 Stochastic thermostat model

For comparison, we have also considered heating of inelastic hard discs (2D) through the “stochastic thermostat”, in the framework of the non-linear homogeneous Boltzmann equation, where the large velocity tail has been shown to behave like $`\mathrm{exp}(v^{3/2})`$ twan : this result has subsequently often been considered as an agreement with the experimental results. In this model, the energy injection is achieved through a random force $`\eta (t)`$ acting on each particle

$$\frac{d𝐯}{dt}=𝐅+\eta (t),\eta _i(t)\eta _j(t^{})=2D\delta (tt^{})\delta _{ij}$$
(3)
where $`D`$ is the amplitude of the injected power and $`𝐅`$ the systematic force due to inelastic collisions. The variance of $`\eta `$ determines the granular temperature in the non equilibrium steady state, but has no influence on the form of the rescaled distribution function $`P(c_x)`$.

With the accuracy of Fig. 3 (the current experimental resolution), the corresponding numerical velocity distributions are then found indistinguishable from a Gaussian, for physically relevant inelasticities in the range $`0.7\alpha _01`$. Departure from Maxwell-Boltzmann behavior becomes manifest below $`\alpha _0=0.6`$ (squares in Fig. 3), which is unphysically low, but the velocity distribution is still incompatible with its experimental counterpart. We also investigated the possibility to describe the effective horizontal dynamics with the 1D stochastic thermostat: for $`0.7<\alpha _0<1`$ the velocity distributions are incompatible with the experimental $`P(v)`$ displayed in Figs. 2, with opposite non Maxwellian features (underpopulated both at vanishing and high energies BBRTW ). More precisely, within the stochastic thermostat approach the kurtosis $`v_x^4/v_x^2^23`$ of the distribution is negative for $`\alpha _0>1/\sqrt{2}0.71`$ twan , irrespective of dimension, which corresponds to an underpopulated low velocity behavior at variance with the experimental data shown in Fig. 2.

We have also considered the stochastic thermostat model for two-dimensional IHS with both tangential $`\alpha _t`$ and normal $`\alpha _n(=\alpha _0)`$ restitution coefficients Luding-McNamara . No numerical studies of $`P(v)`$ can indeed be found in the literature in this case, although an investigation into the non-equipartition between translational and rotational kinetic energies has been performed in Luding-McNamara . The resulting velocity distributions $`P(v)`$ remain very close to a Gaussian for $`\alpha _n0.7`$ and arbitrary $`\alpha _t`$ (where $`1\alpha _t1`$), as for smooth spheres (corresponding to $`\alpha _t=1`$). However, such a two parameter model may be too schematic compared to the experiments Foerster  and we have also considered a more realistic approach with Coulomb friction along the simplifications discussed in Walton : a friction coefficient $`\mu `$, is introduced in addition to ($`\alpha _t,\alpha _n`$Luding\_pre . This does not change significantly $`P(c_x)`$ (see the squares and dashed line in Fig. 3). This seems to discard the relevance of such an approach for the comparison of the velocity distributions with experimental data.

The above analysis shows that the stochastic thermostat in its original formulation (including some possible extensions) does not provide a relevant model of energy injection as far as the velocity distribution is concerned, although it may be useful to investigate other features such as kinetic energy non-equipartition in granular mixtures Granular . However, a variant of this model may improve the picture. In particular, a multiplicative driving \[corresponding to a velocity dependent amplitude $`D|𝐯|^{2\delta }`$ in Eq. (3)\] has been studied in Refs. Cafiero\_prl ; Cafiero . We have performed DSMC simulations for this model, choosing the value of $`\delta `$ that, for a given reasonable inelasticity parameter ($`\alpha =0.9`$), gives the best agreement with the experimental $`P(v)`$. We obtained $`\delta 0.6`$, which leads to the distribution shown by the triangles in Fig. 3. The agreement with the experimental data is satisfactory, and displays a similar accuracy as obtained within the RRC model. Consequently, models with homogeneous energy injection may also describe quite accurately the experimental $`P(v)`$, with the problem of predicting the values of the various parameters involved.

## 5 Conclusion

In conclusion, the Random Restitution Coefficient model with point particles captures the essential features responsible for the observed non Gaussian character of the velocity distribution $`P(v)`$ in vibrated granular gases experiments Rouyer , and represents therefore a self-consistent framework. The conditional distribution $`\mu (\alpha _{1d}|g_x)`$ defining the appropriate collision rule has been shown to encode the relevant dynamic information and provides an alternative route to characterize the non equilibrium steady state, complementary to the direct measure of $`P(v)`$. An interesting point would be to obtain an analytical prediction for $`\mu (\alpha _{1d}|g)`$. It is noteworthy that our approach is mean-field (Boltzmann) like, the only correlations considered being in the collision law. Its self-consistency, which was not an obvious point a priori, has been established by comparison with experiments.

While we have restricted our analysis to the two-dimensional case, such investigations can be extended to three-dimensional systems, for which however experimental measures of effective restitution coefficients seem more difficult. As implemented here, without an analytical knowledge of $`\mu (\alpha _{1d}|g_x)`$, the RRC model is not predictive since an experimental input is required to obtain the correct velocity statistics. Our results however suggest to assess experimentally the question of the universality of $`P(v)`$ from the direct measure of the distribution of effective restitution coefficients. These characteristics are indeed linked within the RRC model: the exponent $`\nu `$, and the whole shape of $`P(v)`$, depend on the functional form of $`\mu (\alpha _{1d}|g_x)`$. At this point, the fact that with a rather poor accuracy, similar $`P`$ and $`\mu `$ have been obtained for the case investigated in Rouyer ; Feitosa  simply confirms the RRC picture, and calls for experiments with widely different collisional properties, such as hollow spheres.

Acknowledgments: We are grateful to K. Feitosa and N. Menon for generous provision of unpublished experimental data and interesting correspondence.


text:
# Preparing encoded states in an oscillator

## I Introduction

It appears, in principle, that the laws of quantum mechanics allow certain mathematical problems to be solved more rapidly than can be done using a classical computer Nielsen and Chuang (2000); Preskill (1998). However, in order to accomplish this task, the state of a quantum system must maintain coherence, despite unwanted interactions with the environment. There have been a number of proposed mechanisms for protecting quantum information during a computation Shor (1995); Steane (1996); Knill and Laflamme (1996); Chau (1997); Rains (1997); Gottesman (1996); Calderbank et al. (1997). Recently, it has been shown Gottesman et al. (2001) that a $`d`$-dimensional quantum system (here we only consider $`d=2`$) can be embedded in an infinite-dimensional Hilbert space, such that a universal set of fault-tolerant quantum gates can be implemented using linear optical operations, squeezing, homodyne detection, and photon counting. The qubits are embedded in the continuous system in a manner which protects the quantum information against small shifts in the canonical quantum variables, $`q`$ and $`p`$. Ideally, the encoded states are an infinite sum of delta functions in both $`q`$ and $`p`$. Of course, such states are non-normalizable, and unphysical. Hence they must be approximated. It has been proposed Gottesman et al. (2001) that these approximate encoded states could be generated by a procedure involving a non-linear interaction Hamiltonian of the form,

$`H^{}`$ $``$ $`qb^{}b,`$ (1)
where $`q`$ is the position operator of one variable, and $`b`$ ($`b^{}`$) is the annihilation (creation) operator of a second variable. Unfortunately, interactions of the form given in Eq. (1) have proven very difficult to implement. They generally require the radiation pressure of photons to move a macroscopic object (a mirror) Giovannetti et al. (2000).

Here we show that approximate encoded states can be generated by coupling the continuous variable to a *single* qubit, and performing a sequence of operations similar to a quantum random walk algorithm Travaglione and Milburn (2002).

In Sec. II, we briefly review the continuous variable encoding scheme proposed in Gottesman et al. Gottesman et al. (2001). In Sec. III we show how approximate encoded states can be non-deterministically generated by coupling the continuous variable to a qubit. We then discuss in Sec. IV how error recovery can be performing by deterministically preparing ancilla variables. Finally, in Sec. V we discuss how an ion trap quantum computer could be used to generate approximate encoded states, and therefore provide an important proof of principle.

## II Encoding a qubit in an oscillator

Quantum computation is generally formulated in terms of interacting two level quantum systems, or qubits. The choice of two level quantum systems is partially because it is easy to draw analogies with the classical bit, but also because a two level system is the simplest non-trivial system; and increasing the number of levels only increases the computation efficiency by a constant of proportionality.

However, with the goal of building a quantum computer in mind, two level quantum systems are by no means the most natural choice. Most physical systems, even in their most elemental form, are represented by many more than two levels. Indeed, many quantum systems are naturally described by a continuous variable (infinite dimensional Hilbert space). Such continuous quantum systems have been well studied, and proposals have been made for performing analog quantum computation using such systems Braunstein (1998); Lloyd and Slotine (1998); Lloyd and Braunstein (1999).

### II.1 Ideal Encoded States

Gottesman et al. Gottesman et al. (2001) discuss how to embed a qubit in a continuous quantum system, so that the extra degrees of freedom within the system can be used to correct errors which arise from unwanted interactions with the environment. Ideally, an encoded zero state, $`|\overline{0}`$, will be represented in position space by the wave function

$`q|\overline{0}`$ $`=`$ $`{\displaystyle \underset{s=\mathrm{}}{\overset{\mathrm{}}{}}}\delta (q2\alpha s)={\displaystyle \underset{s=\mathrm{}}{\overset{\mathrm{}}{}}}e^{i\pi sq/\alpha },`$ (2)
and thus in momentum space, it has the wave function

$`p|\overline{0}`$ $`=`$ $`{\displaystyle \underset{s=\mathrm{}}{\overset{\mathrm{}}{}}}\delta (p{\displaystyle \frac{\pi s}{\alpha }})={\displaystyle \underset{s=\mathrm{}}{\overset{\mathrm{}}{}}}e^{i2sp\alpha }.`$ (3)
Whilst the encoded one state, $`|\overline{1}`$, is represented in position and momentum space by the wave functions,

$`q|\overline{1}`$ $`=`$ $`{\displaystyle \underset{s=\mathrm{}}{\overset{\mathrm{}}{}}}\delta (q2\alpha (s{\displaystyle \frac{1}{2}}))={\displaystyle \underset{s=\mathrm{}}{\overset{\mathrm{}}{}}}e^{i\pi (\frac{sq}{\alpha }1)}`$ (4)
$`p|\overline{1}`$ $`=`$ $`{\displaystyle \underset{s=\mathrm{}}{\overset{\mathrm{}}{}}}(1)^s\delta (p{\displaystyle \frac{\pi s}{\alpha }})={\displaystyle \underset{s=\mathrm{}}{\overset{\mathrm{}}{}}}e^{i(2s1)p\alpha }.`$ (5)
The wave functions for the encoded zero state are depicted in Fig. 1(a), whilst Fig. 1(b) depicts the wave functions for the encoded one state. Clearly the zero and one encoded states are orthogonal,

$`\overline{0}|\overline{1}`$ $`=`$ $`0.`$ (6)
### II.2 Error Recovery

For the details of how quantum computation is performed with these encoded states we direct the reader to Gottesman et al. Gottesman et al. (2001). Here we review the error recovery procedure, which protects these encoded states against shifts in position, $`q`$, and momentum, $`p`$, of size

$`|\mathrm{\Delta }q|<{\displaystyle \frac{\alpha }{2}}`$ $`\mathrm{and}`$ $`|\mathrm{\Delta }p|<{\displaystyle \frac{\pi }{2\alpha }}.`$ (7)
Suppose we have an encoded qubit in some arbitrary superposition of zero and one,

$`|\psi _e`$ $`=`$ $`c_0|\overline{0}+c_1|\overline{1}.`$ (8)
Assume an error occurs to the state $`|\psi _e`$, such that the wave function is shifted in the position variable by some amount $`ϵ<\alpha /2`$. We wish to correct this error without destroying the state. This can be accomplished by using an ancilla variable, prepared in the state

$`|\varphi _a`$ $`=`$ $`(|\overline{0}+|\overline{1})/\sqrt{2},`$ (9)
and an interaction Hamiltonian of the form

$`H_1`$ $`=`$ $`q_ep_a,`$ (10)
where the subscript $`e`$ denotes the encoded qubit variable, and the subscript $`a`$ denotes the ancilla variable. After the two systems have interacted, we can measure the $`q`$ variable of the ancilla system, which will allow us to determine the value of $`ϵ`$. This error can then be corrected by applying an appropriate displacement operation to the encoded qubit system. Likewise, a shift of $`ϵ<\pi /2\alpha `$ in the momentum variable can be corrected using an ancilla system prepared in the $`|\overline{0}`$ state, and evolving according to the interaction Hamiltonian,

$`H_2`$ $`=`$ $`p_eq_a.`$ (11)
## III Preparing encoded states using a qubit

Once prepared, it is hoped that the error recovery procedure will be able to maintain the encoded states. However, preparation of the encoded states is not trivial. As has already been stated, we can only prepare approximate encoded states. In this section we show how approximate encoded states can be prepared with the aid of a single ancilla qubit. Our preparation scheme is non-deterministic, in that a valid approximate encoded state will only be prepared with some probability less than one, however, we will know when our preparation procedure has worked.

We shall denote approximate encoded zero and one states with the symbols $`|\stackrel{~}{0}`$ and $`|\stackrel{~}{1}`$. As in Gottesman et al. (2001), we begin the preparation procedure with the quantum system in the ground state of the oscillator, $`|0`$, and apply squeezing in the $`q`$ quadrature. This creates the state

$`q|s`$ $`=`$ $`g(q,\mathrm{\Delta }),`$ (12)
where

$`g(q,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{e^{q^2/2\mathrm{\Delta }^2}}{\sqrt{\mathrm{\Delta }\sqrt{\pi }}}},`$ (13)
and $`\mathrm{\Delta }`$ is the width of the Gaussian and a measure of the degree of squeezing. $`\mathrm{\Delta }=1`$ corresponds to the oscillator ground state, and $`\mathrm{\Delta }<1`$ indicates a squeezed state. Using an ancilla qubit, initially in the zero state, $`|0`$, the approximate encoded one state, $`|\stackrel{~}{1}_1`$ is then created by applying the sequence of operators,

$`\widehat{H}e^{i\alpha p_e\sigma _z}\widehat{H},`$ (14)
where $`\sigma _z`$ is the Pauli $`z`$ matrix,

$`\sigma _z`$ $`=`$ $`\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),`$ (17)
applied to the qubit, and $`\widehat{H}`$ is the Hadamard gate,

$`\widehat{H}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right),`$ (20)
applied to the qubit. Measuring the qubit in the zero state, which will occur with probability 1/2, results in the continuous variable being left in the state,

$`q|\stackrel{~}{1}_1`$ $`=`$ $`{\displaystyle \frac{N}{\sqrt{2}}}\left(g(q\alpha ,\mathrm{\Delta })+g(q+\alpha ,\mathrm{\Delta })\right),`$ (21)
where $`N`$ is a normalization factor, which is approximately equal to one, if $`\mathrm{\Delta }/\alpha `$ is small compared to one. If the qubit is measured in the one state, the encoded variable is discarded and we try again. To create improved approximate encoded states, we iterate the following procedure:
Given $`|\stackrel{~}{1}_{n1}`$, and a qubit in the state $`|0`$.

* Apply the operators:

$`\widehat{H}e^{i2^{n1}\alpha p_e\sigma _z}\widehat{H}`$ (22)
* Measure the qubit.
* If the qubit is found in the state $`|0`$, then we have created $`|\stackrel{~}{1}_n`$.
* Else discard and start again.

Thus, with probability $`1/2^n`$, we create the approximate encoded state

$`q|\stackrel{~}{1}_n`$ $`=`$ $`{\displaystyle \frac{N}{\sqrt{2^n}}}{\displaystyle \underset{s=1}{\overset{2^n}{}}}g(q+\alpha (1+2^n2s),\mathrm{\Delta }).`$ (23)
In momentum space the approximate encoded state has wave function

$`p|\stackrel{~}{1}_n`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Delta }}{2^n\sqrt{\pi }}}\right)^{1/2}Ne^{(p\mathrm{\Delta })^2/2}{\displaystyle \frac{\mathrm{sin}\alpha 2^np}{\mathrm{sin}\alpha p}}.`$ (24)
Fig. 2 depicts the approximate encoded state $`|\stackrel{~}{1}_3`$, with $`\mathrm{\Delta }=0.15`$ and $`\alpha =\sqrt{\pi /2}`$. This state will be generated with probability 1/8.

The approximate encoded zero state $`|\stackrel{~}{0}_n`$ is created by displacing the state $`|\stackrel{~}{1}_n`$ by an amount $`\alpha `$ in the position variable. Thus

$`q|\stackrel{~}{0}_n`$ $`=`$ $`{\displaystyle \underset{s=1}{\overset{2^n}{}}}g(q+\alpha (2^n2s),\mathrm{\Delta }),`$ (25)
and

$`p|\stackrel{~}{0}_n`$ $`=`$ $`e^{i\alpha p}p|\stackrel{~}{1}_n.`$ (26)
Because of the the $`2^n`$ term in Eq. (22), the average energy of the approximate encoded states will increase exponentially with $`n`$, however, as we see in the following section, the probability of error decreases exponentially with $`n`$.

It is perhaps also worth noting that alternative approximate encoded states, where the sign changes occur in position space rather than momentum space can be created by discarding the states when a $`|0`$ is measured instead of a $`|1`$.

### III.1 Fidelity of approximate encoded states

As in Gottesman et al. (2001), the approximate encoded states $`|\stackrel{~}{0}`$ and $`|\stackrel{~}{1}`$ will have negligible overlap if $`\mathrm{\Delta }`$ is small compared to $`\alpha `$. In position space, the probability of mistaking an approximate encoded zero, $`|\stackrel{~}{0}`$, for an approximate encoded one, $`|\stackrel{~}{1}`$ is simply the probability of measuring the zero state nearer to an odd multiple of $`\alpha `$ than an even multiple. The error probability will be bounded by the sum of each of the Gaussians’ tails,

Error Prob $`<`$ $`2^n2{\displaystyle _{\alpha /2}^{\mathrm{}}}𝑑q\left|{\displaystyle \frac{g(q,\mathrm{\Delta })}{\sqrt{2^n}}}\right|^2.`$ (27)
Thus the error probability is independent of $`n`$, and using the asymptotic expansion of the error function,

$`{\displaystyle _x^{\mathrm{}}}𝑑te^{t^2}`$ $`=`$ $`\left({\displaystyle \frac{1}{2x}}\right)e^{x^2}\left[1O(1/x^2)\right],`$ (28)
it is not hard to show that error probability will be bounded by

Error Prob $`<`$ $`{\displaystyle \frac{4\mathrm{\Delta }}{\sqrt{\pi }\alpha }}e^{\frac{1}{8}\left(\frac{\alpha }{\mathrm{\Delta }}\right)^2}.`$ (29)
Therefore the likelihood of error becomes exponentially small for small $`\mathrm{\Delta }/\alpha `$.

In momentum space, we wish to determine the probability of finding $`(|\stackrel{~}{0}|\stackrel{~}{1})/2`$ closer to an even multiple of $`\pi /\alpha `$ than an odd multiple. Assuming $`N1`$, using Eqs. (24) and (26), we calculate the area under periodic part of the probability function,

$`{\displaystyle \frac{|p|\stackrel{~}{0}_n1|\stackrel{~}{1}_n|^2}{2}}`$ (30)
about each even multiple of $`\pi /\alpha `$, divide this by the width $`2\pi /\alpha `$, and multiple by the area of the Gaussian envelope,

$`{\displaystyle 𝑑pe^{(p\mathrm{\Delta })^2}}.`$ (31)
This gives a bound on the error probability of

Error Prob $`<`$ $`{\displaystyle \frac{1}{\pi 2^{n+1}}},`$ (32)
which becomes exponentially small with $`n`$.

## IV Deterministic error recovery

For robust quantum computation, it is necessary that our encoded states are comb-like in both the position and momentum quadratures. However, this is not necessary for the ancilla systems used in error recovery. To correct an error in position, it is only necessary that the ancilla system is comb-like in position, and to correct an error in momentum it is only necessary that the ancilla system is comb-like in momentum.

This allows us to deterministically prepare ancilla systems for error recovery. The ancilla system states can be prepared using the procedure described in Sec. III, except that we continue with the preparation procedure for $`n`$ iterations, irrespective of whether the qubit is measured in the $`|0`$ or $`|1`$ state. Thus, after three iterations, if the sequence of qubit measurements were say, $`|1,|0`$ and $`|1`$, then we would be left with the state, $`|a`$, depicted in Fig. 3. This state is no longer comb-like in momentum space, but it is still comb-like in position space. Thus, it could be used to perform position error recovery.

Likewise, ancilla variables appropriate for momentum quadrature error recovery can be prepared by first squeezing the vacuum in the momentum quadrature, and replacing the qubit - continuous system interaction operator with

$`e^{i2^{n1}\alpha q_e\sigma _z}.`$ (33)
## V Implementing in an ion trap

There are several physical systems which enable a coupling between a continuous quantum system and a discrete quantum system, such as a cavity QED system or an ion trap. Here we discuss the possibility of creating approximate encoded states in an ion trap.

Though scalable continuous variable quantum computation using ion traps seems unlikely, the ion trap provides a good test bed for such first steps as creating approximate encoded states, as the processes of decoherence within the ion trap are well understood.

Consider a single $`{}_{}{}^{9}\mathrm{Be}_{}^{+}`$ ion, confined in a coaxial-resonator radio frequency (RF)-ion trap, as described in Monroe et al. (1996), and references therein. The continuous quantum system is the vibrational mode of the ion, and the two-level discrete system is the ground and first excited electronic levels of the ion.

First it would be necessary to laser-cool the ion to the motional and electronic ground state, as described in Monroe et al. (1995). Ideally, we would then need to squeeze the vibrational mode of the ion. This could prove a difficult task. However, it is possible to create the sequence of operations described in Eq. (14). The Hadamard operation is accomplished by a $`\pi /2`$-pulse, creating an equal superposition of the ground and excited electronic states. A displacement beam is then applied which excites the motion correlated to the excited state. A $`\pi `$-pulse is then applied to exchange the internal states, and the displacement beam is applied again. Finally another $`\pi /2`$-pulse is applied, executing the second Hadamard gate. The electronic level of the ion is then measured using another laser pulse, tuned to a transition between the first excited level and a higher level. If fluorescence is observed, the ion has been measured in the $`|1`$ state. The absence of fluorescence indicates that the ion is in the ground state. In addition to the operations which we wish to implement, the ion trap system will undergo free evolution, so it will be necessary to couple the qubit, and measure only once every period of oscillation. In order to verify that the desired approximate encoded state had been created it would then be necessary to carry out state tomography on the system.

## VI Conclusions

For a quantum computer to become a reality, the daunting task of providing adequate error correction needs to be fulfilled. At this point in time, it is unclear which, if any, implementation scheme for quantum computation will become viable. As the quantum mechanical oscillator is so prevalent in the study of quantum mechanics, it appears to be a natural test bed for quantum computation. Here we have shown how a continuous quantum system can be coupled to a discrete two level quantum system in a manner which allows the continuous quantum system to encode qubit. The ion trap provides a convenient setting for this encoding scheme as it contains the required discrete and continuous quantum variables.

###### Acknowledgements.

BCT would like to thank R. Polkinghorne, G. Kociuba and P. Cochrane for helpful discussions.


text:
# AGING, DOUBLE HELIX AND SMALL WORLD PROPERTY IN GENETIC ALGORITHMS

## I Introduction

The carrier of genetic information in the Nature, the DNA molecule, has very peculiar structure. This is a spiral-staircase shaped object, whose steps are made from pairs selected from among only four other kinds of molecules. Those are: cytosine (C), guanine (G), adenine (A) and thymine (T). There are only two correct kinds of those pairs, not $`\left(\genfrac{}{}{0pt}{}{4}{2}\right)=6`$, as one might naively suspect. Namely, the cytosine always forms a pair with guanine (C $``$ G), and adenine “likes” thymine (G $``$ T). So, if the sequence in a part of DNA molecule is something like

$$\mathrm{}ATTGTCCA\mathrm{}$$

then the other half of DNA’s double helix has the structure

$$\mathrm{}TAACAGGT\mathrm{}$$

Such a redundancy has not yet been exploited by genetic algorithms, at least the present author failed to find in available literature anything similar in conjunction with this class of optimization algorithms (or any other class …).

In this paper we will show, that the idea of doubled (or duplicated) genetic information may be very useful in fine-tuning of genetic algorithms and in formulating the stopping criteria, which are quite general, independent of the problem under study.

## II The role of mutations

Many elements of existing genetic algorithms are straightforward computer implementations of the natural evolutionary phenomena. We have an evolving population, usually fixed in size, in which individuals mate, have offspring and are mutated. All those processes are driven by one or more random number generators, that is by purely stochastic forces, and, additionally, by the Darwinian rule of survival of the fittest.

It is well known that the crossover processes alone cannot guarantee finding the optimal region in a search space, at least when the population is small. Without mutations, the premature convergence is then almost certain, thus revealing the inability of the algorithm to find the desired solution. The mutations cause rapid changes of location in the search space thus making possible to reach an explore the regions, which are at all not accessible without them.

The majority of researchers and practitioners prefer rather low rate of mutations. This is because the very frequent mutations turn the genetic algorithm into generic Monte Carlo procedure, which is to be avoided, since we hope that some kind of intelligence will lead us much faster to the desired solution than completely blind trials.

## III The fate of mutated genes

Suppose, for simplicity, that the smallest part of the chromosome, called here gene, consists of only a single bit. Let this bit has initial value of $`0`$. The mutation process is then nothing else than negating (flipping) this bit. The question we want to examine now is:

> What will be the state of this bit after $`k`$ generations?

Let the probability of bit flip per generation is equal to $`p`$. We can write

$$\begin{array}{ccccc}\hfill p_{k+1}(1)& =& p_k(0)p\hfill & +& p_k(1)(1p)\hfill \\ \hfill p_{k+1}(0)& =& p_k(0)(1p)\hfill & +& p_k(1)p\hfill \end{array}$$
(1)
where $`k`$ numerates consecutive generations and the number in parenthesis means the state of the bit in question. Taking into account that for any $`k`$

$$p_k(0)+p_k(1)=1$$
(2)
we can rewrite the first row of relation (1) as

$`p_{k+1}(1)`$ $`=`$ $`\left[1p_k(1)\right]p+p_k(1)(1p)`$ (3)
$`=`$ $`p_k(1)\left(12p\right)+p`$ (4)
or, in shorter form

$$p_{k+1}=(12p)p_k+p$$
(5)
where $`p_j`$ is the probability that the examined bit is in the state $`1`$ after $`j`$ generations.

The recurrence formula (5) seems very similar to the construction of a simple random number generator, so one might expect that it produces chaotic sequences of numbers. Surprisingly (?) this is not the case, as we shall see.

Consider the behavior of $`p_k`$ when $`k\mathrm{}`$. It is easy to see, that for $`p=0`$ (no mutations) the sequence $`p_k`$ is constant (and thus convergent); it remains at whatever initial state – in our case $`0`$. The other limiting case, $`p=1`$, is also easy: now the recurrence relation (5) simplifies to

$$p_{k+1}=1p_k$$
(6)
and our bit permanently oscillates between two possible states $`01`$ — there is no convergence. For all other cases, i.e. $`0<p<1`$, the sequence (5) is always convergent and we have

$$\underset{k\mathrm{}}{lim}p_k=\frac{1}{2}$$
(7)
We skip the easy proof. It is interesting, however, that the convergence has oscillatory character when $`p>\frac{1}{2}`$, while for $`p<\frac{1}{2}`$ the sequence $`p_k`$ grows monotonically. There is no chaotic behavior. The case $`p=\frac{1}{2}`$ is special again: it produces constant sequence $`p_k=\frac{1}{2}`$.

## IV Aging of the population

It is widely believed, that the mutations are essentially wrong thing for any individual. They usually damage the genetic material, sometimes to the extent which makes further existence of the individual impossible. It only rarely happens that the mutated individual is better fitted to its environment than the average one in a given population. On the other hand not all mutations are lethal. In the well known Penna ,  model of biological evolution, for example, the subsequent mutations are accumulated in every chromosome. The older the chromosome, the more mutations it carries and finally dies either after single lethal mutation or due to the accumulation of many less damaging defects (we do not discuss here the so called Verhulst factor describing the decrease of the population size due to overcrowding and thus the limited accesss to necessary resources by any single individual). If no replication occurs, then the entire population is eventually extinct. By contrast, the “small” changes in the genotypes are silently transmitted to the offspring chromosomes.

In the Nature, however, an error correcting capability is at work. The double helix structure of the DNA strand limits the proliferation of defective genes. The two parts of DNA (and RNA as well) are not independent, but complementary. After unzipping, just before the replication, only one strand of DNA contains the distorted genetic information, while the other sequence is correct. Therefore only $`50`$% of offspring DNA helices will be damaged.

## V Implementation in genetic algorithms

The idea is to use in genetic algorithms the chromosomes, which are “doubled”, i.e. which consist of regular, well known part (“visible”) and the other (“invisible”) parallel structure, of identical length, which is initially the exact negation of the first part, as below:

| 0 | 0 | 1 | 0 | 0 | 1 | 1 | 0 |  | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 1 | — visible part |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | 1 | 0 | 1 | 1 | 0 | 0 | 1 |  | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | — invisible part |

(8)
We do not mark any gene boundaries here. A single gene may consist of one or more bits; the lengths of consecutive genes need not to be equal to each other. The chromosome becomes mutated when at least one of its bits is flipped.

The genetic material, representing trial solution of the original problem, is placed in the visible part. Every mutation is reflected in the state of this very part, the invisible part being completely insensitive to mutations. Quite contrary, the crossover operation is performed on both parts, mixing independently two visible and two invisible parts of the involved chromosomes, at the same crossover point(s). It is easy to see, that the crossover operation alone, no matter one- or multipoint, will never destroy the complementarity of the visible and invisible part of any chromosome. Of course, the fitness factor will always be computed for the visible part only. Now the process of aging of the population, generation after generation, may be tracked by looking at the invisible part of every chromosome. We can easily count all the bits in the entire population, which were changed by mutations — it is enough to compare the visible and invisible parts of all chromosomes and count the bits, which are identical in both parts. Of course, the even number of mutations applied to the same bit will go unnoticed. This is sometimes called backward evolution. From the earlier considerations, and assuming that the mutation probability $`p`$ per bit and per generation is small (more precisely $`p<\frac{1}{2}`$, as is always the case), we can conclude that the fraction of (effectively) mutated bits in the entire population will be an increasing function of time. The stochastic limit for this fraction is equal to $`\frac{1}{2}`$. The value of $`p=\frac{1}{2}`$ uniquely separates the genetic algorithms from Monte Carlo type of optimization.

## VI Discussion

For $`p\frac{1}{2}`$ the fraction of mutated bits initially increases linearly with the number of generations. This observation does not take into account directly neither the structure of the individual chromosome (number of genes or bits it consists of) nor the number of chromosomes in the population. Such information is hidden in the value of parameter $`p`$ — probability of a particular bit being flipped during a single cycle of simulation (single epoch). Thus the the number of generations to reach $`50`$% of flipped bits may be roughly estimated as $`N_{gen}=\frac{1}{2p}`$. In practice, as the performed simulations show, the fraction of effectively mutated bits reaches the value of $`\frac{1}{2}`$ — for the first time — after roughly $`5`$ to $`6`$ times greater number of epochs. The exact analytical formula for $`p_k`$ as a function of $`k`$ is neither easy to obtain nor very informative. For any $`k`$, $`p_k`$ is a polynomial of order $`k`$ in variable $`p`$. The graph of $`p_k`$ vs. $`k`$, however, shows striking similarity to the graph of the expression:

$$\frac{1}{2}\left[1\mathrm{exp}\left(2kp\right)\right]$$
(9)
This is only an observation based on several evolutionary processes, guessed rather than formally derived. Anyway, utilizing this approximation, we can say that after $`\frac{1}{2p}`$ epochs the population is mature and gradually looses its innovative forces, and after

$$N_{gen}=\frac{3}{p}$$
(10)
epochs it becomes old and practically useless. It is time then to switch to other locally searching routine to improve further the optimal set of unknown parameters, if appropriate. The formula (10) is the main result of this paper. It sets the upper limit for the number of necessary generations.

In conclusion, during calculations we should monitor the behavior of the fraction of mutated bits. When this variable reaches the value of $`\frac{1}{2}`$, then we may say that every second bit in the population was flipped at least once. This is the other way of saying that the search space has been explored quite thoroughly and no significant improvement of the fitness should be expected.

Why do we make such a claim?

If the fitness of all chromosomes was the same, i.e. there was no evolutionary pressure of any kind, then after $`\frac{3}{p}`$ generations the search space would be covered quite uniformly, although irregularly, with trial points. Further sampling would only make finer the already existing irregular grid of trial points. When there is a preference for better fitted chromosomes, then after $`\frac{1}{2p}`$ generations all interesting regions, or at least their neighborhoods, should have already been found. Observing the values of the fitness function alone, either as an average for the population or for the best individual only, may be very misleading.

Considering the efficiency of the algorithm, understood as the number of generations necessary to find the optimal fitness, we can see that it is inversely proportional to $`p`$. For every practical purpose the condition

$$pN_{bits}=1$$
(11)
(average number of bits flipped in the whole population during one epoch is equal to one) sets the lowest sensible limit for $`p`$. $`N_{bits}`$ is the total number of bits in the population (in its visible part). Of course, higher values of $`p`$ will speed up the evolutionary search.

Is the prescription $`N_{gen}\frac{1}{2p}`$ perfect?

Certainly not. It gives the reasonable estimate of the number of generations needed in fairly regular cases. The genotype space with only one well fitted chromosome and all other equally bad is an evident exception. What can be done in such cases is to increase the number of bits assigned to every continuous unknown (kind of oversampling), for the price of increased computational effort, of course, since this approach is equivalent to the use of finer grid of points in the search space. To purely combinatorial problems (only integer and/or logical unknowns) our estimate should be of similar value.

As an example let us take the data from . Here the population consisted of 30 chromosomes, $`20`$ bits long each. and $`p`$ was set equal to $`0.03`$. According to (11), $`p`$ should be set as at least $`1/(30\times 20)0.0017`$ — while the quoted mutation rate was $`18`$ times higher than that number. The population should become mature after some $`1/(2\times 0.03)17`$ generations and the evolutionary search should be terminated after $`100`$ generations at the latest. In fact, the author of  reports that satisfactory results were achieved after $`10`$$`19`$ generations (every run was limited to $`20`$ generations). The search space had only $`2^{20}=1048576`$ points, so this example may be considered small and not representative. Nevertheless, no more than only $`600`$ evaluations of the fitness factor were enough to reach valuable conclusions.

## VII Why is genetic algorithm so effective anyway?

We must be aware, that the search space, no matter how large, is always discrete and finite for this class of optimization algorithms. It may be considered as a random graph, in which every chromosome is a vertex (we are not interested in the edges). This graph is by no means random, but clearly exhibits the small world property, i.e. the average (Hamming) distance between its vertices scales as the logarithm of their number. Indeed, $`N`$-bit chromosomes are points in $`2^N`$ element universe. The maximum distance between any two chromosomes is equal to $`N`$, so the average distance must never exceed this number. Incidentally $`\mathrm{log}_22^N=N`$. This fact alone explains why genetic algorithms are fairly insensitive with respect to the number of unknowns. On the other hand, speaking of the neighboring points in the search space makes sense, especially when we think of nearest neighbors. It is therefore intuitively appealing that we should be able to find a (hopefully small) subset of points in the search space with the property that any point is distant from this set no more than $`1`$ unit – something similar to the backbone or spanning tree known from the graph theory. In case $`N=3`$ it is easy to check that exactly two points, namely $`\left(000\right)`$ and $`\left(111\right)`$, are enough to form such a subset with requested property. Evaluating fitness for each member of this subset should be nearly equivalent to the exhaustive (“brute force”) search, since then, for arbitrarily chosen point in the search space either this point itself or one of its nearest neighbors was visited and evaluated during evolutionary process. Unfortunately, today we don’t know how to find such a set in general case; we don’t even know what is its cardinality – hopefully much lower than that of original search space. It is certain, however, that the solution of this problem need not to be unique.111Consider the simple case of $`N=2`$ -bit chromosomes. The universe consists of just $`4`$ elements: $`\left(00\right)`$, $`\left(01\right)`$, $`\left(10\right)`$ and $`\left(11\right)`$. There are two subsets with desired property: $`\{\left(01\right),\left(10\right)\}`$ and $`\{\left(00\right),\left(11\right)\}`$. Perhaps the genetic algorithm is the best available tool for approximate construction of such subsets?

## VIII Comments on performance

Using “doubled” chromosomes we successfully mimic the double helix structure of DNA. The cost is moderate: we need to double the storage for the population. Counting the flipped bits is performed only once per generation, so this cost should be negligible in comparison with evaluations of fitness function. Instead of the true implementation of the “double helix”, one can use the simplified formula (10) as a stopping criterion. Direct observation of the fraction of effectively mutated bits will signal the end of calculations usually much earlier.

## IX Acknowledgment

This work was done as part of author’s statutory activity in Institute of Physics, Polish Academy of Sciences.