text:
# Dynamical properties of vibrofluidized granular mixtures

## I Introduction

The present keen interest in the dynamical properties of granular materials is motivated both by the challenge of understanding the complex processes involved and by the important practical applications in engineering, industry and technology general . These materials are peculiar in many respects and display several intriguing phenomena such as clustering Zanetti , shear instability Brey ; Baldassa  and lack of energy equipartition, which make their behavior different from ordinary molecular fluids. The dissipation of kinetic energy during the inelastic collisions makes them special. The main motivation of the present paper stems from two recent experiments Feitosa ,Parker  which demonstrated that when a mixture constituted by two different species of grains is vibrated, each component attains its own “granular temperature”, i.e. the average kinetic energy per particle does not take on the universal value $`fKT`$, where $`f`$ is the number of degrees of freedom and $`K`$ a constant, as it occurs in molecular gases. On the contrary, one observes that the ratio $`T_1/T_2`$ varies with concentration, inelasticity parameters, particle sizes, masses and driving mechanism. Even in the absence of energy injection, the inelastic gas cools, but one observes that the temperature ratio asymptotically remains constant. On the other hand, it is understood that while the only relevant hydrodynamic field is the global temperature $`T=xT_1+(1x)T_2`$, transport properties depend on that ratio dufty2 .

These are of course manifestations of the non-equilibrium nature of Granular systems, which can only be maintained stationary by a continuous energy feeding to compensate the energy losses due to the inelastic collisions and to friction.

A theoretical understanding of such a behavior of granular mixtures has been achieved in the case of homogeneous driving mechanisms by means of a combination of models and approximations including the pseudo-Maxwell inelastic gas and the Inelastic Hard Sphere model treated by means of the Boltzmann-Enskog equation. Both models have been studied analytically and numerically in the free cooling dufty ; ioepuglio  and in the driven case ioepuglio2 ; Pagnani ; Barrat ; barrat2 . Apart from the studies of ref Pagnani ; barrat2  none of these investigations considered the role played by the gravitational field, by the strongly inhomogeneous boundary conditions employed in the experiments of refs. Feitosa ; Parker , by the roughness of the grains and by their rotational degrees of freedom. Moreover, there is still an open debate about the “best” energy feeding mechanism. Whereas theoreticians seem to favour a uniform thermal gaussian bath, because it lends itself to a great deal of analytical work, a numerical computer experiment can test directly driving mechanisms which are closer to those employed in a laboratory. The structure of the paper is the following: in section 2 we illustrate briefly the model, leaving the technical details to the appendix; in section 3 we discuss the results for the geometry of ref. Feitosa . In sec. 4 we consider a different aspect ratio, namely a longer box, where gravity plays a more relevant role. Finally in sec. 5 we present our conclusions.

## II Model system

We decided to remain as close as possible to a commonly employed experimental set-up, by constraining the grains to move on a vertical rectangular domain of dimensions $`L_x\times L_y`$. The gravitational force acts along the negative $`y`$ direction. The grains are assumed to be spherical and free to rotate about an axis normal to the $`xy`$ plane. They receive energy by colliding with the horizontal walls, harmonically vibrating at frequency $`\nu `$. The side walls instead are immobile and were chosen either smooth or rough according to the numerical experiment. When side walls are considered to be rough, they are assumed to have the same friction coefficient $`\mu `$ as the particles.

The collisional model adopted in the present paper corresponds to the one proposed by Walton Walton . It conserves both the linear and the angular momentum of a colliding pair, but allows energy to be dissipated by means of a normal restitution coefficient and a friction coefficient $`\mu `$. The collision rule (given in detail in the appendix) takes into account a reduction of normal relative velocity of the two particles ($`V_n`$), a reduction of total tangential relative velocity ($`V_r`$) and an exchange of energy between those two degrees of freedom. The reduction of normal relative velocity is modeled by means of a non constant restitution coefficients $`\alpha _{ij}[0,1]`$, whose dependence by the relative velocity is of the form:

$$\alpha _{ij}(V_n)=\{\begin{array}{ccc}1(1r_{ij})\left(\frac{|V_n|}{v_0}\right)^{\frac{3}{4}}\hfill & \text{ for }\hfill & V_n<v_0\hfill \\ r_{ij}\hfill & \text{ for }\hfill & V_n>v_0\hfill \end{array}$$

where $`i`$ and $`j`$ are the numbers indicating the species of the colliding particles, $`r_{ij}`$ are constants related to the three types of colliding pairs, $`v_0\sqrt{gd}`$, where $`d`$ is the average diameter of the particles and $`g`$ the gravitational acceleration Bizon .

Simulated collisions are of two types: with sliding or sticking point of contact. When the following condition is satisfied (high relative tangential velocity), the collision happens in a sliding fashion, otherwise it is sticking:

$$\frac{|V_r|}{V_n}\frac{l+1}{l}\mu (1+\alpha _{ij})$$
(1)
where $`l`$ is the dimensionless moment of inertia (equal to $`2/5`$ for spheres), while $`\mu `$ is a static friction coefficient characterizing the surface roughness of particles, assumed equal to the dynamical friction coefficient. The full dynamics consists of interparticle collisions, and wall-particle collisions. The trajectories between collision events are parabolic arcs due to the presence of the gravitational field.

An efficient Event Driven (ED) simulation code was employed to evolve the system Cattuto .

The two species were chosen to be spheres of equal diameters $`d=0.16cm`$, and unequal masses $`m_1=1.5810^2g`$ and $`m_2=5.2110^3g`$, respectively. The driving frequency was set to $`50`$ Hz, the vibration amplitude $`A=3.5`$ diameters so that the corresponding dimensionless acceleration $`\mathrm{\Gamma }=A\omega ^2/g=56`$.

All the averages quantities reported in the following have been obtained by employing equally spaced data points separated by time intervals $`\mathrm{\Delta }t=10^1`$ s in order to assure statistical independence of measures. We performed a number $`n=1.510^5`$ of vibration cycles.

Barrat and Trizac barrat2  have recently considered one of the systems studied in the present paper. However, our treatment presents some differences:

* The bottom and top walls of barrat2  move in a sawtooth manner with a negligible excursion so that their positions are considered fixed, while our walls move sinusoidally with a non-negligible amplitude.
* The walls are not smooth in our treatment, but have a friction coefficient $`\mu >0`$.
* The collisional models of our treatment and that of Barrat-Trizac are different (they do not distinguish between sticking and sliding collisions).
* We take into account gravity.

## III Results for a short system

The box dimensions were $`L_x=48d`$ and $`L_y=32d`$. Simulation runs were carried out using $`N_1=150`$ grains of each species. The static friction coefficient has been always chosen as $`\mu =0.1`$. The stationary state is determined by the balance between the energy input provided by the vibrating walls at a frequency of $`50`$ Hz and the dissipation due to inelastic collisions. The typical collision frequencies are of the order of $`\nu _1580`$ Hz and $`\nu _2850`$ Hz for the heavy and light balls, respectively. A typical microscopic configuration of the system is shown in fig.1.

### III.1 Temperature profiles

The restitution coefficients were first set to $`r_{11}=0.93`$ for 1-1 collisions, $`r_{22}=0.7`$ for 2-2, and $`r_{12}=0.8`$ for 1-2 collisions. This means that the more massive particles are also the more elastic ones. In figure 2 we show the partial translational temperature profiles for the two species, and observe that close to the vertical boundaries the two temperatures are essentially determined by the energy injected by the vibrating walls. Indeed, interparticle collisions are rare within this region, and play no signifcant role because of the low local density (see figure 4). In addition, grains 1 and 2 impinging with the same speed on the mobile wall bounce with the same velocity ($`VA\omega `$), hence the local value of the temperature ratio, $`\gamma =T_2/T_1`$ near the vibrating walls, turns out to be approximately $`\gamma m_1/m_2`$, as shown in figure 2.

On the other hand, the temperature drops as the distance from the walls increases, while the ratio $`\gamma `$ grows up to a plateau value, indicating that collisions tend to cool the mixture and render the two partial temperatures closer. Figure 2 clearly displays the breakdown of the kinetic energy equipartition already noticed in previous experimental and theoretical studies.

We also measured the rotational temperature profiles, shown in figure 2. We observe for rotational temperatures the same kind of equipartition breakdown that holds for translational degrees of freedom. Moreover, the absolute values of rotational and translational temperatures are quite different, as already reported by Luding Luding  for a one-component, vibrated granulate. On the other hand, the ratio of the two rotational temperature profiles seems to be quite close to that of the translational temperature profiles.

In figure 3, we changed the restitution coefficients and set $`r_{11}=0.7`$, $`r_{22}=0.93`$ for 2-2, and $`r_{12}=0.81`$. Therefore now the more massive particles are the more inelastic. We observe that, whereas the temperature profiles near the vibrating planes are nearly unchanged, because collisions are rare, the value of the temperature of the heavier species is lower and the temperature ratio is closer to $`1`$. In this case the larger inelasticity of the heavier particles competes with the mass asymmetry which instead tends to make $`T_2/T_1`$ smaller.

### III.2 Area fraction profiles

Due to the large value of the parameter $`\mathrm{\Gamma }`$, the partial density profiles tend to be rather symmetric, with a maximum near the center. We also notice small differences in the density profiles, which reveal that the heavier species has higher concentration close to the center of the cell, while the lighter species is more spread.

Comparing the present results for the area fraction profiles with those recently obtained by employing a Direct Simulation Monte Carlo technique Pagnani , we notice that the agreement is only qualitative, whereas the temperature ratios are in significantly better agreement.

We verified that removing tangential friction (coupling the translational and rotational degrees of freedom) does not change significantly the above scenario. In fact, the translational temperature and density profiles of the cases $`\mu =0`$ and $`\mu =0.1`$ show only small quantitative differences, more pronounced close to the horizontal walls.

### III.3 Transversal profiles

We also studied the density and temperature profiles along the horizontal direction. To the best of out knowledge, no such measure has been reported in experimental works. In figures 5 and 6 we observe that temperature profiles vanish close to side walls, while density profiles display their maxima in the same region. In order to gain further insight, we analyzed a sequence of snaphsots of the dynamics, and observed that the system bears a denser cloud of grains in the vicinity of one of the side walls. Such a configuration was maintained over an interval of time much longer than the vibration period $`2\pi /\omega `$.

The cluster eventually “evaporates” to form again close to a randomly selected side wall. Over several periods of oscillation of the cell, we noticed an effective horizontal simmetry breaking, i.e. the number of particles in the right-hand and left-hand sides of the cell were rather different. The system is unstable with respect to horizontal density fluctuations and clusterizes spontaneously, until the vibrating bases wash out the cluster.

Moreover, we also observed some spontaneous tendency of the system to segregate the species, a fact which becomes more apparent at high densities. Both phenomena have their origin in the inelasticity. A possible qualitative explanation of the observed dynamics is as follows: particles with smaller restitution coefficient tend to group together, since the more energy they dissipate through 2-2 collisions, the denser the segregated domain becomes. Particles with higher restitution coefficient bounce for a longer time after a collision, and dilate more quickly.

### III.4 Velocity Probability Distributions

Figure 7 shows the distributions functions for the translational velocities along the horizontal direction. For sake of comparison, all the velocities were rescaled by their mean square values. We notice that the two transversal velocity distribution functions deviate from a Gaussian and are fitted by $`f(c)=A/(exp(\beta c^\alpha )+exp(\beta c^\alpha ))`$, with: $`A=1.068`$, $`\beta =1.74`$, $`\alpha =0.918`$ .

The exponential tails appear to have a smaller slope than the theoretically predicted 3/2 value, for a uniform system stochastically driven vannoije . We notice the phenomenon already reported in Pagnani : the rescaled velocity distribution of lighter species display higher tails.

Figure 7 also shows the angular velocity distributions (rescaled by their mean square value, see caption).

Our simulations show that such distributions have pronounced non-Gaussian tails. The angular velocity distributions can be described by the same scaling function we used for translational velocity distributions.

### III.5 Collision time distribution

Recently Blair and Kudrolli utilizing high speed digital photography measured the collision statistics of grains bound to move on an inclined plane blair\_poster . They determined the distribution of path lengths, $`P(l)`$ and showed that it deviates from the theoretical prediction for elastic hard spheres. In particular, $`P(l)`$ shows a peak in the small $`l`$ region, not present in elastic systems. In order to assess the existence of such a behavior in the IHS model we performed similar measurements. Figure 8 displays the results for the collision times of the two species. One clearly sees that the probability density that a particle suffers a collision in a short interval is enhanced with respect to the elastic case. The physical reason for that lies in the existence of strong correlations which lead to the presence of clusters where the path lengths are shorter than in a uniform system. For the sake of comparison we plotted in figure 8 also the corresponding distribution of an elastic gas having the same density and granular temperature of our inelastic mixture. The effect of the shorter average collision time is clearly visible. We also notice that the ratio of the collision frequencies is approximately equal to the the ratio of the average velocities, as expected from elementary kinetic arguments.

## IV Long system

In the case of a box of dimensions $`L_x=48d`$ and $`L_y=64d`$ with $`300`$ grains of each species the effect of gravity is much more evident, resulting in a stronger inhomogeneity and asymmetry of the system. In fact, the density and temperatures profiles are not symmetric with respect to the vertical direction. A configuration of such system is shown in fig. 9. Most of the particles remain suspended above the bottom wall and are hit by those which are between the bottom wall and and the bulk. Very few particles reach the upper vibrating wall, so that the granular temperature of the system is much lower in the top than in the bottom. The partial temperature profiles and density profiles are shown in fig. 10. One sees that near the lower vibrating wall the temperature profiles are similar to those of the shorter system, whereas in the bulk they are appreciably different. The average ratio is lower than in the $`L_y=32`$ d case. An interesting feature present in fig.10 is the presence of a region of increasing temperatures. Such a phenomenon has been predicted by the hydrodynamic theory of Brey et al. Brey . The theory predicts a temperature varying as $`Tz^{3/2}`$ in the region above the minimum. Such a prediction is verified by our simulation results.

Physically, the increase is due to the competition between the energy input from the walls and the small dissipation due to the reduced number of collisions associated with the low density region. This can be appreciated in the inset of figure 10 where the quantity $`\xi (z)=n(z)T(z)^{3/2}`$ is shown: $`\xi (z)`$ is proportional to the energy dissipation due to inelastic collisions among particles, being the collision rate $`nT^{1/2}`$ and the average dissipated energy in a single collision $`T`$. The graph of $`\xi (z)`$ reveals that the energy dissipation is much more relevant in the bottom region than in the middle and upper regions.

Correspondingly the velocity PDF shows an interesting behavior (see figure 11). In fact, the shape of the rescaled PDF becomes narrower with the height. In the central region the measured PDF resembles the exponential PDF measured in the short system (figure 7), whereas in the bottom region the PDF has more extended tails. Such a lack of universality in the velocity PDf was noted experimentally by Blair and Kudrolli blair . It was also noticed in puglisi1 : there it was shown that the tails of the PDF became broader when the dissipation rate was increased. Here the mechanism is similar: the broader PDFs are those measured in the bottom region, where the dissipation rate is higher.

## V Conclusion

Summarizing, we studied a systems of inelastic particles in vertically vibrated containers and subject to the gravitational field. The system was numerically investigated by using an Event-Driven dynamics, which affords the exploration of a wide range of dynamical parameters. We examined two similar setups which differ only for the aspect ratio of the container and for the number of particles.

In the shorter system we determined the partial granular temperatures, their ratio, and the area fraction profiles along the vertical direction. Those measures are in qualitative agreement with the experimental results, but a quantitative comparison requires a more detailed knowledge of the experimental parameters. Physically the lack of energy equipartition between the two species has two causes: a) the vibrating walls feed energy proportionally to the mass of each species; b) the dissipation is in general different for the two components. The two effects may conspire (for example if the heavier species is the more elastic) to give a small value of the temperature ratio or they may cancel and give a temperature ratio close to one (for example when the heavier species is the more inelastic). Moreover, we observed that density profiles are non-uniform along the horizontal direction, as well, indicating that the particles tend to clusterize in the vicinity of side walls. We also measured the velocity distributions, verifying that they can be collapsed onto each other, under proper rescaling, and the scaling function has a stretched exponential behavior. Finally the distribution of flight times between succesive collisions has been measured and compared to that of an elastic system: the inelasticity has the effect of enhancing the statistics of very short times.

In the longer system we have again obtained the temperature and area fraction profiles, observing a stronger inhomogeneity and asymmetry. In this case the area fraction is much larger near the bottom wall, reaching higher values than in the previous experiment. The middle-upper region seems to be qualitatively well described by recent hydrodynamic theories developed for one component systems. Here the velocity PDFs display a non-universal behavior, with broader tails in the more dense (and dissipative) regions.

## VI Appendix

In order to make it simpler for the reader to interpret the present model, we present an appendix in which we explicitly state the collision rulesFarkas .

The colliding particles are characterized by radii $`R_1`$ and $`R_2`$, positions $`𝐫_1`$ and $`𝐫_2`$, translational velocities $`𝐯_1`$ and $`𝐯_2`$ and rotational (angular) velocities $`𝝎_1`$ and $`𝝎_2`$ (we assume that if $`𝝎`$ is parallel and in the direction of the $`z`$ axis, than the rotation is anticlockwise if seen from above the $`xy`$ plane). We introduce the normal unitary vector joining the centers of the particles $`𝐧=(𝐫_2𝐫_1)/|𝐫_2𝐫_1|`$ and the tangential unitary vector $`𝐭`$ obtained rotating $`𝐧`$ by an anti-clockwise angle $`\pi /2`$. Then we introduce the relative velocity $`𝐠`$, the velocity of the center of mass $`𝐕`$ and the velocities of the particles in the center of mass frame $`𝜻_1`$ and $`𝜻_2`$:

$`𝐠`$ $`=𝐯_1𝐯_2`$ (2a)
$`𝐕`$ $`={\displaystyle \frac{m_1𝐯_1+m_2𝐕_2}{m_1+m2}}`$ (2b)
$`𝜻_1`$ $`={\displaystyle \frac{m_{eff}}{m_1}}𝐠`$ (2c)
$`𝜻_2`$ $`={\displaystyle \frac{m_{eff}}{m_2}}𝐠`$ (2d)
where $`m_{eff}=m_1m_2/(m_1+m_2)`$.

Then we decompose the relative velocity $`𝐠`$ on the orthonormal basis given $`𝐧`$ and $`𝐭`$, as well as the velocities of the particles in the center of mass frame, i.e.:

$`g_n`$ $`=(𝐯_\mathrm{𝟏}𝐯_\mathrm{𝟐})𝐧`$ (3a)
$`g_t`$ $`=(𝐯_\mathrm{𝟏}𝐯_\mathrm{𝟐})𝐭`$ (3b)
$`\zeta _{n\beta }`$ $`={\displaystyle \frac{m_{eff}}{m_\beta }}g_n`$ (3c)
$`\zeta _{t\beta }`$ $`={\displaystyle \frac{m_{eff}}{m_\beta }}g_t`$ (3d)
with $`\beta =1,2`$ the index of the particle.

We finally introduce $`g_c`$ as the relative circular velocity at the point of contact and $`g_r`$ as the total tangential relative velocity (circular and translational) at the point of contact:

$`g_c`$ $`=R_1\omega _1+R_2\omega _2`$ (4a)
$`g_r`$ $`=g_c+g_t`$ (4b)
To characterize the collision rules we use a model that take into account a reduction of normal relative velocity ($`g_n`$), a reduction of total tangential relative velocity ($`g_r`$) and an exchange of energy between those two degrees of freedom. The reduction of normal relative velocity is modeled as usual by means of a restitution coefficient $`\alpha [0,1]`$:

$$g_n^{}=\alpha g_n.$$
(5)
We assume a dependence of $`\alpha `$ by the relative velocity of the form:

$$\alpha (g_n)=\{\begin{array}{ccc}1(1\alpha _0)\frac{|g_n|}{v_0}^{\frac{3}{4}}\hfill & \text{ per }\hfill & g_n<v_0\hfill \\ \alpha _0\hfill & \text{ per }\hfill & g_n>v_0\hfill \end{array}$$

where $`\alpha _0`$ is a constant which stands for $`r_{ij}`$, $`v_0=\sqrt{gD}`$, with $`D`$ the average diameter of the particles (and $`g`$ is the gravity acceleration). From equations (5) we obtain the update of normal velocities in the center of mass frame:

$`\zeta _{n1}^{}`$ $`=\alpha {\displaystyle \frac{m_{eff}}{m_1}}g_n`$ (6a)
$`\zeta _{n2}^{}`$ $`=\alpha {\displaystyle \frac{m_{eff}}{m_2}}g_n`$
What lacks now is an expression for the tangential and angular velocities after the collisions. We distinguish between two possible cases: sliding or sticking collisions. The condition that allows to determine if a collision is sticking or sliding is the following:

$$\frac{|g_r|}{g_n}\frac{l+1}{l}\mu (1+\alpha )$$
(7)
where $`l`$ is the adimensionalized inertia moment and is equal to $`1/2`$ or $`2/5`$ if the particle is a disk or a sphere respectively, while $`\mu `$ is the static friction coefficient of the surface of the particles which in the following will be assumed to be equal to the dynamic friction coefficient.

In the sliding case we use the following rules to update the tangential components of the velocities of the particles in the center of mass frame:

$`\zeta _{t1}^{}`$ $`=\zeta _{t1}\mu (1+\alpha ){\displaystyle \frac{m_{eff}}{m_1}}g_nsign(g_r)`$ (8a)
$`\zeta _{t2}^{}`$ $`=\zeta _{t2}+\mu (1+\alpha ){\displaystyle \frac{m_{eff}}{m_2}}g_nsign(g_r)`$ (8b)
$`R_1\omega _1^{}`$ $`=R_1\omega _1{\displaystyle \frac{\mu (1+\alpha )}{l}}{\displaystyle \frac{m_{eff}}{m_1}}g_nsign(g_r)`$ (8c)
$`R_2\omega _2^{}`$ $`=R_2\omega _2{\displaystyle \frac{\mu (1+\alpha )}{l}}{\displaystyle \frac{m_{eff}}{m_2}}g_nsign(g_r)`$ (8d)
In the case of a sticking collision, instead, the update rules are obtained considering that:

$$\zeta _{t1}^{}\zeta _{t2}^{}+R_1\omega _1^{}+R_2\omega _2^{}=0$$
(9)
from which, after calculations, one gets:

$`\zeta _{t1}^{}`$ $`={\displaystyle \frac{1}{l+1}}\zeta _{t1}{\displaystyle \frac{l}{l+1}}{\displaystyle \frac{m_{eff}}{m_1}}(R_1\omega _1+R_2\omega _2)`$ (10a)
$`\zeta _{t2}^{}`$ $`={\displaystyle \frac{1}{l+1}}\zeta _{t2}+{\displaystyle \frac{l}{l+1}}{\displaystyle \frac{m_{eff}}{m_2}}(R_1\omega _1+R_2\omega _2)`$ (10b)
$`R_1\omega _1^{}`$ $`=R_1\omega 1\left[{\displaystyle \frac{l}{l+1}}+{\displaystyle \frac{m_{eff}}{(l+1)m_2}}\right]`$ (10c)
$`R_2\omega _2{\displaystyle \frac{m_{eff}}{(l+1)m_1}}{\displaystyle \frac{1}{l+1}}\zeta _{t1}`$ (10d)
$`R_2\omega _2^{}`$ $`=R_2\omega 2\left[{\displaystyle \frac{l}{l+1}}+{\displaystyle \frac{m_{eff}}{(l+1)m_1}}\right]`$ (10e)
$`R_1\omega _1{\displaystyle \frac{m_{eff}}{(l+1)m_2}}+{\displaystyle \frac{1}{l+1}}\zeta _{t2}`$ (10f)
The velocity of particles in the absolute frame are finally obtained in the two cases (considering that the center of mass is not perturbed by the collision) by the equation:

$$𝐯_\beta ^{}=𝐕+\zeta _{n\beta }^{}𝐧+\zeta _{t\beta }^{}𝐭$$
(11)
(for particle of index $`\beta `$) leading to the following global collision rule for translational velocities:

$`𝐯_1^{}`$ $`=𝐯_1(1+\alpha ){\displaystyle \frac{m_2}{m_1+m_2}}[(𝐯_1𝐯_2)𝐧]𝐧\{\begin{array}{ccc}sign(g_r)\mu (1+\alpha )\frac{m_2}{m_1+m_2}[(𝐯_1𝐯_2)𝐧]𝐭\text{ (sliding)}\hfill & & \\ \frac{m_1lm_2}{(l+1)(m_1+m_2)}[(𝐯_1𝐯_2)𝐭]𝐭\frac{l}{l+1}\frac{m_2}{m_1+m_2}(R_1\omega _1+R_2\omega _2)\text{ (stick)}\hfill & & \end{array}`$ (12c)
$`𝐯_2^{}`$ $`=𝐯_2+(1+\alpha ){\displaystyle \frac{m_1}{m_1+m_2}}[(𝐯_1𝐯_2)𝐧]𝐧+\{\begin{array}{ccc}sign(g_r)\mu (1+\alpha )\frac{m_1}{m_1+m_2}[(𝐯_1𝐯_2)𝐧]𝐭\text{ (sliding)}\hfill & & \\ \frac{m_2lm_1}{(l+1)(m_1+m_2)}[(𝐯_1𝐯_2)𝐭]𝐭+\frac{l}{l+1}\frac{m_1}{m_1+m_2}(R_1\omega _1+R_2\omega _2)\text{ (stick)}\hfill & & \end{array}`$ (12f)
while for rotational velocities:

$`R_1\omega _1^{}`$ $`=\{\begin{array}{ccc}R_1\omega _1\frac{\mu (1+\alpha )}{l}\frac{m_2}{m_1+m_2}g_nsign(g_r)\text{ (stick)}\hfill & & \\ R_1\omega 1\left[\frac{l}{l+1}+\frac{m_1}{(l+1)(m_1+m_2)}\right]R_2\omega _2\frac{m_2}{(l+1)(m_1+m_2)}\frac{1}{l+1}\zeta _{t1}\text{ (slide)}\hfill & & \end{array}`$ (13c)
$`R_2\omega _2^{}`$ $`=\{\begin{array}{ccc}R_2\omega _2\frac{\mu (1+\alpha )}{l}\frac{m_1}{m_1+m_2}g_nsign(g_r)\text{ (stick)}\hfill & & \\ R_2\omega _2\left[\frac{l}{l+1}+\frac{m_2}{(l+1)(m_1+m_2)}\right]R_1\omega _1\frac{m_1}{(l+1)(m_1+m_2)}+\frac{1}{l+1}\zeta _{t1}\text{ (slide)}\hfill & & \end{array}`$ (13f)
###### Acknowledgements.

This work was supported by Ministero dell’Istruzione, dell’Università e della Ricerca, Cofin 2001 Prot. 2001023848, by INFM and INFM Center for Statistical Mechanics and Complexity (SMC).


format: text

text:
# A Two-Player Game of Life

## 1 Introduction

Conway’s game of life  is the best known two-dimensional cellular automaton , due to the complex behaviour it generates from a simple set of rules. The game of life, introduced to the world at large by Martin Gardner in 1970 , has provided throughout the years challenging problems to many enthusiasts. A large number of these are documented in Stephen Silver’s comprehensive life lexicon . Detailed accounts of the game of life have been given by Poundstone in , which uses the game of life as an illustration of complexity, and by Sigmund in , which uses the game of life in the context of artificial life and the ability of a cellular automaton to self-replicate. An interesting extension of the game of life to three dimensions was suggested by Bays .

The game of life is played on a square lattice with interactions to nearest and to next-nearest neighbours, where each cell can be either empty or occupied by a token and is surrounded by eight neighbouring cells. The evolution of the game is governed by the following simple rules. If a cell is empty it gives birth to a token if exactly three of its neighbours are occupied, and if it is occupied it survives if either two or three of its neighbours are occupied. In all other cases either the cell remains empty or it dies, i.e. it becomes empty. The game evolves by repeated applications of these rules to produce further configurations. The single player of the game decides what the initial configuration of the lattice will be and then watches the game evolve. One of the questions that researchers have investigated is the asymptotic behaviour of this evolutionary process. It transpires that when the initial configuration is random and its density is high enough, then the game eventually stabilises to a density of about 0.0287 ; see also .

One aspect that is missing from Conway’s game of life is the competitiveness element of two-player games, as Gardner noted “Attempts have also been made to invent competitive games based on “Life”, for two or more players, but so far without memorable results.” \[11, p. 255\]. This is the challenge that we take up in this paper.

An interesting version of the game of life for two players is known as black and white or immigration , where cells are either white or black and when a birth occurs the colour of the token is decided according to the majority of neighbouring cells. Although this variation is interesting in its own right, survival does not involve the two colours and remains non-competitive as in the single player version. We propose a new two-player version of the game of life where both birth and survival are competitive, and provide a preliminary analysis of its behaviour.

The rest of the paper is organised as follows. In Section 2 we present the rule set for our two-player version of the game of life and in Section 3 we provide a mean-field analysis of the game. In Section 4 we present results of simulations to ascertain the asymptotic density of the game and, finally, in Section 5 we give our concluding remarks and provide a web link to an implementation of the game.

## 2 Rules of the Game

In the two-player variation of the game of life, which we call p2life, the players, white and black, are competing for space. Conway’s game of life is considered to be “interesting” since its simple set of rules lead to complex and unpredictable behaviour. (The notion that from simple rules complex behaviour can emerge, which may help us understand the diversity of natural phenomena, is discussed in great detail by Wolfram in  but can already be learned from the 130 year old thesis of van der Waals on liquid-vapour equilibria.) P2life maintains the “interesting” behaviour of the game of life by preserving the essence of Conway’s game and adding to its rules a competitive element to decide who will give birth and who will survive.

The rules of p2life, from white’s point of view (the rules from black’s point of view are symmetric), are as follows:

If a cell is empty, then we consider two cases:

1. The cell has exactly three white neighbours and the number of black neighbours is different from three. In this case a white token is born in the cell.
2. The cell has exactly three white and three black neighbours. In this case an unbiased coin determines whether a white or black token is born in the cell.

If a cell is occupied by a white token, then we consider two cases:

1. If the difference between the number of white and black neighbours is two or three, then the white token survives.
2. If the difference between the number of white and black neighbours is one and the number of white neighbours is at least two, then the white token survives.

It is clear that if there is only one colour on the lattice then p2life reduces to the standard one-player version of the game. We note that having non-symmetric rules for white and black would allow us to investigate how different behaviour sets interact, but finding such a set of rules which would be “interesting” is an open problem. During the process of deciding the rule set for p2life we considered several variations, which we now briefly discuss:

1. In the first part of the birth rule, insisting that the number of black tokens should be less than three, which is similar to birth in the black and white variant. This decreases the birth rate but still seems to be “interesting”.
2. Omitting the second part of the birth rule, i.e. the random choice when there are exactly three white and three black neighbours. Again, this decreases the birth rate but still seems to be “interesting”.
3. Omitting the second part of the survival rule, i.e. survival when the difference between the number of white and black tokens is one and the number of white tokens is at least two. This decreases the survival rate and the conflict between the two players but, once again seems “interesting”
4. In the second part of the survival rule omitting the side condition that the number of white neighbours must be at least two, or in the first part of the survival rule omitting the condition that the difference between white and black is at most three. Due to the increased survival rates, these modifications lead to steady growth with spatial boundaries between the two players. It would be interesting to compare these rules to Schelling’s models of segregation  or the Ising model .

We illustrate two configurations which lead to interesting confrontations between white and black. The configuration shown on the left-hand side of Table 1 leads to a black block and two white gliders as shown on the right-hand side of the table. In the one-player game of life this initial configuration annihilates itself into empty space. On the other hand, the configuration shown on the left-hand side of Table 2 leads to two black blocks as shown on the right-hand side of the table. Thus black wins from this position. In the one-player game of life this initial configuration leads to six blinkers.

## 3 Mean-Field Theory for P2life

To compute the initial configuration given a specified initial density of $`p`$, we use the following procedure for each cell in the lattice. Firstly we determine whether it is occupied or empty according to the given probability $`p`$ and then, if it should be occupied we toss an unbiased coin to decide whether it will be occupied by a white or black token. In a similar fashion to the game of life , we now use mean-field theory to determine the density of the tokens after applying the rules to the initial configuration.

It can be verified by the rules of p2life that the density $`p^{}`$, after a single application of the rules to a square lattice having initial density $`p`$, is given the mean-field equation,

$$p^{}=21p^363p^4+\frac{105}{2}p^5+35p^6\frac{1505}{16}p^7+\frac{1057}{16}p^8\frac{553}{32}p^9.$$
(1)
From (1) we can compute the maximum density of $`p^{}`$, which is $`0.3895`$ at an initial density of $`p=0.6206`$. This can also be seen from the plot of the mean-field equation shown in Figure 1. As we would expect, simulations show a very close match with this plot, although obviously due to correlations we cannot use the mean-field equation to predict the long-term behaviour of the game. An interesting point to note about the plot is that when $`p=1`$, $`p^{}=0.2188`$, which is different from the original game of life (and the black and white variant) where $`p^{}=0`$. The reason is that although dense regions of a single colour die off immediately, mixed spaces between the two colours allow for survival of tokens. We also note that from Figure 1 we can verify that the only fixed-point of the mean-field equation for p2life is zero. To improve the prediction power for repeated applications of the rules we could extend the mean-field analysis using the local structure theory of Gutowitz and Victor .

## 4 Asymptotic Density for P2life

We have investigated the properties of p2life with particular emphasis on the estimation of its asymptotic density via simulations using matlab. This approach offers flexibility and reasonable performance (0.7 million updates per second) for a particularly computationally intensive task. We have performed simulations on both periodic (toroidal) and cutoff boundary conditions on square lattices of sizes from $`10\times 10`$ to $`1000\times 1000`$. In particular, we have investigated the dependence of the asymptotic density $`p_{\mathrm{}}`$ on the initial density $`p`$ of a random, uniformly and independently distributed, initial configuration. Each configuration is iterated until the game configuration reaches a stable or oscillatory state. In general, it is straightforward to detect when such conditions occur, with the notable exception of the case of periodic boundary conditions when in the final configuration there exists a glider which travels around the torus without colliding with other tokens in occupied cells.

The results for a square lattice of size $`100\times 100`$ are shown in Figure 2. Each experimentally computed point plotted in Figure 2 is the average over one hundred independently selected initial configurations. It is immediately evident that in direct contrast to the one-player version of the game of life, the asymptotic density of p2life is non-zero for non-zero initial density. In the one-player version of the game, if the initial density increases above approximately $`0.85`$ the asymptotic density is zero due to the annihilation of all tokens during the first iteration due to overcrowding. In fact, for $`p=1`$ we estimate that p2life has asymptotic density $`p_{\mathrm{}}=0.0362`$. This is due to the fact that in p2life the survival rules allow members of both token populations to carry over to the next iteration with $`p^{}=0.2188`$, which is consistent with the value obtained via the mean-field approach.

We make several other interesting observations on the asymptotic behaviour of p2life. Firstly, we have estimated that the maximum asymptotic density of $`p_{\mathrm{}}=0.0362`$ is reached at about $`p=0.3`$ and remains fairly constant until $`p=1`$. This is different to the behaviour of the one-player version of the game, where a fairly constant asymptotic density $`p_{\mathrm{}}=0.0287`$ is reached at about $`p=0.15`$ and remains at that level until approximately $`p=0.70`$ . Secondly, when periodic boundary conditions are used, the asymptotic density $`p_{\mathrm{}}`$ increases by approximately 5% to $`0.0381`$ compared to the situation of cutoff boundary conditions. Thirdly, it appears that the size of the lattice does not seem to affect the asymptotic density, although, as in the one-player version there may be small finite-size effects . Finally, we estimated the ratio of the loser population over the winner population at the final state. A histogram of the results aggregated over 400 runs for a square lattice of size $`100\times 100`$, where 100 runs were carried out for initial densities of 0.25, 0.5, 0.75 and 1.00, is shown in Figure 3. We observe that over 69% of the runs resulted in the ratio being over 0.5, i.e. the loser having more that one third of the final population.

## 5 Concluding Remarks

Our main contribution is to have shown that, by injecting competitive elements into Conway’s game of life, “Life” can be “interesting” when more than one player participates in the game. An applet demonstrating p2life can be accessed at

http://hades.dcs.bbk.ac.uk/users/gr/p2life/p2life.asp.

A problem that we are now investigating is how to convert p2life into a “real” game, i.e. where players are allowed to make moves between generations, which change the configuration of the tokens, and devise strategies to overpower or live side-by-side with their opponent.


text:
# Complex Ginzburg–Landau equation for self–oscillatory systems under amplitude modulated forcing in a 1:1 resonance

## I Introduction

The periodic forcing of spatially extended self-oscillatory systems is a classical method to excite the formation of spatial patterns in such systems. This kind of forcing admits a universal description when the system is operated near the oscillation threshold, and forcing acts on a $`n:m`$ resonance, defined by the relation $`\omega _e=\left(n/m\right)\left(\omega _0+\nu \right)`$ between the external forcing frequency $`\omega _e`$ and the natural frequency of oscillations $`\omega _0`$, where $`n/m`$ is an irreducible integer fraction and $`\nu `$ is a small mistuning. In such case the slowly varying complex amplitude of the oscillations $`u`$ verifies the following generalized complex Ginzburg-Landau equation (CGLE)

$$_tu=a_1u+a_2^2u+a_3\left|u\right|^2u+a_4\overline{u}^{n1},$$
(1)
where $`\overline{u}`$ stands for the complex conjugate of $`u`$, $`a_{i=1,2,3}`$ are complex coefficients and $`a_4`$ is proportional to the $`m`$th power of the forcing amplitude. Eq. (1) is valid in principle for perfectly periodic forcings and has been introduced making use of symmetry arguments . That elegant reasoning is based on the existence of a discrete time invariance when a perfectly periodic forcing is applied to a system. However if, e.g., the amplitude of forcing is modulated in time in an arbitrary way all temporal symmetries are broken and Eq. (1) is not rigorously justified. In this work we prove that Eq. (1) does provide a universal description of the close to threshold dynamics of self-oscillatory extended systems forced in a 1:1 resonance ($`n=1`$ and $`m=1`$ in Eq. (1)). The analysis is based on the technique of multiple scales and generalizes the concept of resonant forcing as it considers (almost) periodic forcings that are nonuniform across the system and/or modulated in time. Our derivation hence allows a rigorous simplified study of both noisy forcings as well as multi-frequency forcings within a 1:1 resonance. Generalizations to other resonances ($`n1`$ or $`m1`$), although cumbersome, can be made straightforwardly following the lines of the present derivation.

## II Model

We consider a generic two-dimensional system described by $`N`$ real dynamical variables $`\left\{U_i(𝐫,t)\right\}_{i=1}^N`$ whose time evolution is governed by the following set of real equations written in vector form,

$$_t𝐔(𝐫,t)=𝐟(\mu ,\alpha ;𝐔,^2𝐔),$$
(2)
where $`𝐟`$ is a sufficiently differentiable function of its arguments. We assume that the dependence of $`𝐟`$ on spatial derivatives of $`𝐔`$ is through its Laplacian $`^2=_x^2+_y^2`$, $`𝐫=(x,y)`$. This is the simplest dependence on derivatives in rotationally invariant systems at the time it corresponds to physical systems of most relevance like, e.g., reaction-diffusion and nonlinear optical systems. $`\mu `$ is the bifurcation parameter and $`\alpha (𝐫,t)`$ is the forcing parameter, which is allowed to vary on time and space. Physically $`\alpha `$ may represent either an independent parameter, or the modulated part of any other parameter.

We assume that in the absence of forcing ($`\alpha =0`$) Eq. (2) supports a steady, spatially homogeneous state $`𝐔=𝐔_s\left(\mu \right)`$ $`\left(_t𝐔_s=_x𝐔_s=_y𝐔_s=0\right)`$, which looses stability at $`\mu =\mu _0`$ giving rise to a self-oscillatory, spatially homogeneous state. In other words, we assume that the reference state $`𝐔_s`$ suffers a homogeneous Hopf bifurcation at $`\mu =\mu _0`$. We wish to study the small amplitude solutions that form in the system close to the bifurcation when the arbitrary parameter $`\alpha `$ is modulated in time with a frequency close to that of the free oscillations.

For the sake of convenience we introduce a new vector

$$𝐮(𝐫,t)=𝐔(𝐫,t)𝐔_s,$$
(3)
which measures the deviation of the system from the reference state, in terms of which we rewrite Eq. (2) as a Taylor series,

$`_t𝐮`$ $`=`$ $`𝐅(\mu ,\alpha )+𝒥(\mu ,\alpha )𝐮+𝒟(\mu ,\alpha )^2𝐮`$ (5)
$`+𝐊(\mu ,\alpha ;𝐮,𝐮)+𝐋(\mu ,\alpha ;𝐮,𝐮,𝐮)+\mathrm{h}.\mathrm{o}.\mathrm{t}.,`$
where $`\mathrm{h}.\mathrm{o}.\mathrm{t}.`$ denotes terms of higher order than $`3`$ in $`𝐮`$ or than $`1`$ in $`^2𝐮`$. These $`\mathrm{h}.\mathrm{o}.\mathrm{t}.`$ have no influence near the bifurcation $`\left(\mu \mu _0\right)`$ since, as we show below, they are $`𝒪\left(\left|\mu \mu _0\right|^2\right)`$ or smaller and only terms up to $`𝒪\left(\left|\mu \mu _0\right|^{3/2}\right)`$ contribute to the leading order dynamics of the system whenever the bifurcation is supercritical, which is the case we assume.

The different elements of the expansion (5) are defined as

$`𝐅(\mu ,\alpha )`$ $`=`$ $`𝐟_s,`$
$`𝒥_{ij}(\mu ,\alpha )`$ $`=`$ $`\left[f_i/U_j\right]_s,𝒟_{ij}(\mu ,\alpha )=\left[f_i/^2U_j\right]_s,`$
$`𝐊(\mu ,\alpha ;𝐚,𝐛)`$ $`=`$ $`\frac{1}{2!}_{i,j=1}^N\left[^2𝐟/U_iU_j\right]_sa_ib_j,`$
$`𝐋(\mu ,\alpha ;𝐚,𝐛,𝐜)`$ $`=`$ $`\frac{1}{3!}_{i,j,k=1}^N\left[^3𝐟/U_iU_jU_k\right]_sa_ib_jc_k,`$
where $`𝐚,𝐛,𝐜`$ are arbitrary vectors and the subscript $`s`$ denotes $`𝐔=𝐔_s\left(\mu \right)`$. Vector $`𝐅(\mu ,\alpha )`$ is subjected to the condition

$$𝐅(\mu ,0)=0,$$
(6)
since in the absence of forcing ($`\alpha =0`$) the reference state $`𝐮=0`$ is a steady state of Eq. (5) by hypothesis. $`𝒥`$ and $`𝒟`$ are matrices, and vector $`𝐊`$ ($`𝐋`$) is a symmetric and bilinear (symmetric and trilinear) function of its two (three) last arguments.

## III The Hopf bifurcation

In the absence of forcing the stability of the reference state against small perturbations $`\delta 𝐮`$ is governed by the following equation

$$_t\delta 𝐮=𝒥(\mu ,0)\delta 𝐮+𝒟(\mu ,0)^2\delta 𝐮,$$
(7)
obtained upon linearizing Eq. (5) for $`\alpha =0`$ with respect to $`\delta 𝐮`$. The general solution to Eq. (7) is a superposition of plane waves of the form

$$\delta 𝐮(𝐫,t)=_j𝐰_j\mathrm{exp}\left(\mathrm{\Lambda }_jt\right)\mathrm{exp}\left(i𝐤_j𝐫\right),$$
(8)
with

$`\mathrm{\Lambda }_j𝐰_j`$ $`=`$ $`(\mu ,k_j^2)𝐰_j,`$ (9)
$`(\mu ,k^2)`$ $`=`$ $`𝒥(\mu ,0)k^2𝒟(\mu ,0),`$ (10)
($`k^2=𝐤𝐤`$) hence eigenvalues and eigenvectors of matrix $``$ depend on $`𝐤`$ only through $`k^2`$ as $``$ does. As we are assuming that the reference state looses stability at $`\mu =\mu _0`$ via a homogeneous Hopf bifurcation in the absence of forcing, matrix $`(\mu ,k^2)`$ must have a pair of complex-conjugate eigenvalues $`\{\mathrm{\Lambda }_1,\mathrm{\Lambda }_2\}=\{\lambda (\mu ,k^2),\overline{\lambda }(\mu ,k^2)\}`$ (the overbar denotes complex conjugation) governing the instability, i.e.:

1. Close to the bifurcation $`Re\mathrm{\Lambda }_{i3}<0`$ whilst $`Re\lambda `$ can become positive for some $`k`$’s,
2. At the bifurcation $`Re\lambda `$ is maximum and null at $`k=0`$ (the perturbation with largest growth rate is spatially homogeneous):

$$Re\lambda _0=0,\left(_kRe\lambda \right)_0=0,\left(_k^2Re\lambda \right)_0<0.$$
(11)
where, here and in the following,

subscript $`0`$ affecting functions denotes $`\left\{\mu =\mu _0,\alpha =0,k=0\right\}\begin{array}{c}.\end{array}`$

1. The instability is oscillatory, i.e.,

$$Im\lambda _0=\omega _00.$$
(12)

Finally, the preceding properties imply that:

1. All eigenvalues of $`_0=𝒥_0`$, see Eq. (10), have negative real part but $`\{\lambda _0,\overline{\lambda }_0\}=\{i\omega _0,i\omega _0\}`$.

For the sake of later use we introduce the right and left eigenvectors of $`𝒥_0`$ associated with eigenvalues $`\{i\omega _0,i\omega _0\}`$,

$$\begin{array}{cc}𝒥_0𝐡=i\omega _0𝐡,& 𝒥_0\overline{𝐡}=i\omega _0\overline{𝐡},\\ 𝐡^{}𝒥_0=i\omega _0𝐡^{},& \overline{𝐡^{}}𝒥_0=i\omega _0\overline{𝐡^{}},\end{array}$$
(13)
where the short-hand notation $`𝐡=𝐰_1\left(\mu =\mu _0,k^2=0\right),\overline{𝐡}=𝐰_2\left(\mu =\mu _0,k^2=0\right)`$ has been introduced. These vectors verify the following orthonormality relations:

$`𝐡^{}\overline{𝐡}=0,𝐡^{}𝐡=1\text{.}`$

## IV Scales

We are interested in determining the small amplitude solutions that from in the system close to the Hopf bifurcation, which we define by

$$\mu =\mu _0+\epsilon ^2\mu _2,$$
(14)
where $`\epsilon `$ is a smallness parameter $`\left(0<\epsilon 1\right)`$. The study is based on the widely used technique of multiple scales. These spatial and temporal scales appear naturally close to the bifurcation and are those on which the asymptotic dynamics of the system evolves. As is well known, in a homogeneous Hopf bifurcation these slow scales are given by

$$T=\epsilon ^2t,X=\epsilon x,Y=\epsilon y.$$
(15)
These scales follow from the behaviour of $`\lambda `$ close to the bifurcation, Eq. (14), for values of $`k`$ close to the most unstable mode $`k=0`$:

$$\lambda (\mu _0+\epsilon ^2\mu _2,k^2)=\lambda _0+\left(_\mu \lambda \right)_0\epsilon ^2\mu _2+\frac{1}{2}\left(_k^2\lambda \right)_0k^2+\mathrm{max}\{𝒪\left(\epsilon ^4\right),𝒪\left(k^2\right)\},$$
(16)
where the term $`\left(_k\lambda \right)_0k`$ has not been included since $`\left(_k\lambda \right)_0=0`$ as $`\lambda `$ is an even function of $`k`$. From this equation we obtain, making use of Eq. (11),

$`Re\lambda (\mu _0+\epsilon ^2\mu _2,k^2)=\left(_\mu Re\lambda \right)_0\epsilon ^2\mu _2{\displaystyle \frac{1}{2}}\left|_k^2Re\lambda \right|_0k^2+\mathrm{max}\{𝒪\left(\epsilon ^4\right),𝒪\left(k^2\right)\},`$

what indicates that the only modes which can experience linear growth verify $`k=𝒪\left(\epsilon \right)`$, hence the asymptotic dynamics of the system exhibits spatial variations on a scale $`x,yk^1\epsilon ^1`$and the slow spatial scales $`\left(X=\epsilon x,Y=\epsilon y\right)`$ follow. Thus, putting $`k=\epsilon k_1`$,

$`Re\lambda (\mu _0+\epsilon ^2\mu _2,\epsilon ^2k_1^2)=\epsilon ^2\left[\left(_\mu Re\lambda \right)_0\mu _2{\displaystyle \frac{1}{2}}\left|_k^2Re\lambda \right|_0k_1^2\right]+𝒪\left(\epsilon ^4\right),`$

which shows that the growth of perturbations occurs on a scale $`t\left(Re\lambda \right)^1\epsilon ^2`$ and the slow timescale $`T=\epsilon ^2t`$ follows. On the other hand, making use of Eq. (12) and putting $`k=\epsilon k_1`$ again, Eq. (16) yields

$`Im\lambda (\mu _0+\epsilon ^2\mu _2,\epsilon ^2k_1^2)=\omega _0+\epsilon ^2\left[\left(_\mu Im\lambda \right)_0\mu _2+{\displaystyle \frac{1}{2}}\left(_k^2Im\lambda \right)_0k_1^2\right]+𝒪\left(\epsilon ^4\right),`$

whose first term, $`\omega _0=𝒪\left(\epsilon ^0\right)`$, indicates that the original timescale $`t`$ must be retained, whilst the rest of terms do not introduce other relevant timescales.

As for the external forcing we assume that its form is consistent with the previous scales. In particular we assume that it is weak and of the form

$`\alpha (𝐫,t)`$ $`=`$ $`\epsilon ^3\alpha _3(X,Y,T,t),`$ (17)
$`\alpha _3(X,Y,T,t)`$ $`=`$ $`{\displaystyle _{p=1}^{\mathrm{}}}\left[\alpha _{3,p}(X,Y,T)\mathrm{exp}\left(ip\omega _0t\right)+\overline{\alpha }_{3,p}(X,Y,T)\mathrm{exp}\left(ip\omega _0t\right)\right],`$ (18)
Thus we are considering a resonant forcing with slowly varying amplitude. Notice that, according to Eqs. (17,18), $`\alpha (𝐫,t+\frac{2\pi }{\omega _0})=\alpha (𝐫,t)+𝒪\left(\epsilon ^5\right)`$ and hence the considered forcing is almost periodic of fundamental frequency $`\omega _0`$.

Under these conditions a multiple scale analysis is possible and we look for asymptotic solutions to Eq. (5) in the form

$$𝐮(𝐫,t)=_{m=1}^{\mathrm{}}\epsilon ^m𝐮_m(X,Y,T).$$
(19)
We finally introduce Eqs. (14) and (1719) into Eq. (5) making use of the following chain rules for differentiation

$`_t𝐮`$ $`=`$ $`{\displaystyle _{m=1}^{\mathrm{}}}\epsilon ^m\left(_t𝐮_m+\epsilon ^2_T𝐮_m\right),`$
$`^2𝐮`$ $`=`$ $`{\displaystyle _{m=1}^{\mathrm{}}}\epsilon ^{m+2}\left(_X^2+_Y^2\right)𝐮_m,`$
and solve at increasing orders in $`\epsilon `$.

## V The complex Ginzburg–Landau equation

The general form of Eq. (5) at any order $`\epsilon ^m`$ is found to be

$$𝔍\left(𝐮_m\right)=𝐠_m(X,Y,T,t),$$
(20)
where

$$𝔍\left(𝐮\right)_t𝐮𝒥_0𝐮,$$
(21)
and $`𝐠_m`$ does not depend on $`𝐮_m`$ (but on $`𝐮_{n<m}`$). Clearly, as $`𝔍\left(\mathrm{exp}\left(i\omega _0t\right)𝐡\right)=𝔍\left(\mathrm{exp}\left(i\omega _0t\right)\overline{𝐡}\right)=0`$, see Eq. (13), the solvability of Eq. (20) requires

$$_t^{t+2\pi /\omega _0}𝑑t^{}𝐡^{}𝐠_m(X,Y,T,t^{})\mathrm{exp}\left(i\omega _0t^{}\right)=0,$$
(22)
(or its equivalent complex-conjugate) which ensures that $`𝐠_m`$ does not contain secular terms (proportional to $`\mathrm{exp}\left(i\omega _0t\right)𝐡`$ or to $`\mathrm{exp}\left(i\omega _0t\right)\overline{𝐡}`$). Once condition (22) is verified, the asymptotic solution to Eq. (20) reads

$`𝐮_m(X,Y,T,t)`$ $`=`$ $`u_m(X,Y,T)\mathrm{exp}\left(i\omega _0t\right)𝐡+\overline{u}_m(X,Y,T)\mathrm{exp}\left(i\omega _0t\right)\overline{𝐡}`$ (24)
$`+𝐮_m^{}(X,Y,T,t),`$
where $`u_m(X,Y,T)`$ is a function of the slow scales, and the last term is the particular solution. Note that the solution (24) should involve, in principle, terms proportional to all the eigenvectors of $`𝔍()`$ \[which are those of $`𝒥_0`$, see Eq. (9)\]. However all of them are damped according to $`\mathrm{exp}\left[\left|Re\mathrm{\Lambda }_i(\mu _0,0)\right|t\right]`$, since $`Re\mathrm{\Lambda }_{i3}\left(\mu =\mu _0,k=0\right)<0`$ by hypothesis, except those associated with $`(𝐡,\overline{𝐡})`$.

### A Order $`\epsilon `$

This is the first nontrivial order and reads

$$𝐠_1=0.$$
(25)
The solvability condition (22) at this order is fulfilled and, according to Eq. (24),

$$𝐮_1=u_1(X,Y,T)\mathrm{exp}\left(i\omega _0t\right)𝐡+\overline{u}_1(X,Y,T)\mathrm{exp}\left(i\omega _0t\right)\overline{𝐡}.$$
(26)
### B Order $`\epsilon ^2`$

At this order

$$𝐠_2=\left(_\mu 𝐅\right)_0\mu _2+𝐊(\mu _0,0;𝐮_1,𝐮_1).$$
(27)
The first term of the r.h.s. is null by virtue of Eq. (6). Making use of Eq. (26) and taking into account that $`𝐊`$ is symmetric and bilinear in its two last arguments, Eq. (27) can be written as

$`𝐠_2`$ $`=`$ $`2𝐊(\mu _0,0;𝐡,\overline{𝐡})\left|u_1\right|^2+𝐊(\mu _0,0;𝐡,𝐡)u_1^2\mathrm{exp}\left(i2\omega _0t\right)`$ (29)
$`+𝐊(\mu _0,0;\overline{𝐡},\overline{𝐡})\overline{u}_1^2\mathrm{exp}\left(i2\omega _0t\right).`$
The solvability condition (22) is automatically fulfilled again and, according to Eq. (24),

$`𝐮_2`$ $`=`$ $`u_2(X,Y,T)\mathrm{exp}\left(i\omega _0t\right)𝐡+\overline{u}_2(X,Y,T)\mathrm{exp}\left(i\omega _0t\right)\overline{𝐡}`$ (31)
$`+𝐯_0\left|u_1\right|^2+𝐯_2u_1^2\mathrm{exp}\left(2i\omega _0t\right)+\overline{𝐯}_2\overline{u}_1^2\mathrm{exp}\left(2i\omega _0t\right),`$
where,

$`𝐯_0`$ $`=`$ $`2𝒥_0^1𝐊(\mu _0,0;𝐡,\overline{𝐡}),`$ (32)
$`𝐯_2`$ $`=`$ $`\left(𝒥_0i2\omega _0\right)^1𝐊(\mu _0,0;𝐡,𝐡),`$ (33)
are constant vectors, and $``$ is the $`N\times N`$ identity matrix. Note that both $`𝒥_0`$ and $`\left(𝒥_0i2\omega _0\right)`$ are invertible since neither $`0`$ nor $`2i\omega _0`$ are eigenvalues of $`𝒥_0`$ by hypothesis: otherwise other eigenvalues different from $`\{i\omega _0,i\omega _0\}`$ would have null real part at the bifurcation.

### C Order $`\epsilon ^3`$

Finally, at this order we find

$`𝐠_3`$ $`=`$ $`_T𝐮_1+\mu _2\left(_\mu 𝒥\right)_0𝐮_1+𝒟_0\left(_X^2+_Y^2\right)𝐮_1`$ (35)
$`+\left(_\alpha 𝐅\right)_0\alpha _3+2𝐊(\mu _0,0;𝐮_1,𝐮_2)+𝐋(\mu _0,0;𝐮_1,𝐮_1,𝐮_1).`$
Application of the solvability condition (22) yields, after substituting Eqs. (18), (26) and (31) into Eq. (35), and making use of the symmetry and linearity properties of vectors $`𝐊`$ and $`𝐋`$,

$$_Tu_1=c_1\mu _2u_1+c_2\left(_X^2+_Y^2\right)u_1+c_3\left|u_1\right|^2u_1+c_4\alpha _{3,1},$$
(36)
where

$`c_1`$ $`=`$ $`𝐡^{}\left(_\mu 𝒥\right)_0𝐡,`$ (37)
$`c_2`$ $`=`$ $`𝐡^{}𝒟_0𝐡,`$ (38)
$`c_3`$ $`=`$ $`2𝐡^{}𝐊(\mu _0,0;𝐡,𝐯_0)+2𝐡^{}𝐊(\mu _0,0;\overline{𝐡},𝐯_2)+3𝐡^{}𝐋(\mu _0,0;𝐡,𝐡,\overline{𝐡}),`$ (39)
$`c_4`$ $`=`$ $`𝐡^{}\left(_\alpha 𝐅\right)_0,`$ (40)
are constant coefficients.

Finally note that, making use of Eqs. (19), (15) and (26), the asymptotic state of the system can be written as

$`𝐮(𝐫,t)`$ $`=`$ $`u(𝐫,t)\mathrm{exp}\left(i\omega _0t\right)𝐡+\overline{u}(𝐫,t)\mathrm{exp}\left(i\omega _0t\right)\overline{𝐡}+𝒪\left(\epsilon ^2\right),`$ (41)
$`u(𝐫,t)`$ $`=`$ $`\epsilon u_1(X,Y,T,t),`$ (42)
hence $`u`$ denotes the leading order amplitude of oscillations. Its evolution equation is obtained from Eq. (36) by returning to original scales and parameters via Eqs. (14) and (15), and reads

$$_tu(𝐫,t)=\left(\mu \mu _0\right)c_1u+c_2^2u+c_3\left|u\right|^2u+f(𝐫,t),$$
(43)
where

$$f(𝐫,t)=c_4\epsilon ^3\alpha _{3,1}(X,Y,T)$$
(44)
is a slowly varying function proportional to the complex amplitude of the fundamental component of forcing $`\epsilon ^3\alpha _{3,1}`$, see Eqs. (17,18).

Eq. (43), or (36), is a complex Ginzburg–Landau equation containing an inhomogeneous forcing term, $`f(𝐫,t)`$, which generalizes Eq. (1) for $`n=m=1`$. Eq. (43) is valid whenever $`Rec_30`$ (supercritical bifurcation) since otherwise it could lead to unbounded solutions. If $`Rec_3>0`$ the bifurcation is subcritical and the analysis must incorporate higher orders in the $`\epsilon `$–expansion.

## VI Application to two-frequency forcing

Just for the sake of illustration let us finally consider the two-frequency forcing case. In particular we assume that the forcing parameter $`\alpha `$ has the form

$$\alpha (𝐫,t)=A\left(\omega _1t\right)+A\left(\omega _2t\right),$$
(45)
where $`A`$ is a $`2\pi `$ periodic function $`\left[A\left(\theta \right)=A\left(\theta +2\pi \right)\right]`$. This means that we are dealing with the superposition of two forcings of equal amplitude but of different frequencies ($`\omega _1`$ and $`\omega _2`$). An example of this type of forcing could be the illumination of a photosensitive version of the Belousov-Zhabotinsky reaction with two equal sources, whose light intensities are periodically modulated in time at two different frequencies. Another example could be the injection into a laser cavity of two coherent fields of equal amplitudes and different frequencies.

Due to the commented periodicity one can write

$$A\left(\omega _it\right)=_{p=1}^{\mathrm{}}\left[A_p\mathrm{exp}\left(ip\omega _it\right)+\overline{A}_p\mathrm{exp}\left(ip\omega _it\right)\right].$$
(46)
(Note that the term $`p=0`$ is absent, as in Eqs. (17,18), since it would correspond to a constant bias. This bias term, if any, is implicitly considered in our analysis at the same level as any other constant parameter of the system and it would appear, in principle, in the determination of the Hopf bifurcation). Upon rewriting Eq. (46) as

$$A\left(\omega _it\right)=_{p=1}^{\mathrm{}}\left[A_p\mathrm{exp}\left(ip\delta _it\right)\mathrm{exp}\left(ip\omega _0t\right)+\overline{A}_p\mathrm{exp}\left(ip\delta _it\right)\mathrm{exp}\left(ip\omega _0t\right)\right],$$
(47)
where $`\delta _i=\omega _i\omega _0`$, we can express Eq. (45) in the form (17,18):

$$\alpha (𝐫,t)=_{p=1}^{\mathrm{}}\left[a_p\left(t\right)\mathrm{exp}\left(ip\omega _0t\right)+\overline{a}_p\left(t\right)\mathrm{exp}\left(ip\omega _0t\right)\right],$$
(48)
where

$$a_p\left(t\right)=A_p\left[\mathrm{exp}\left(ip\delta _1t\right)+\mathrm{exp}\left(ip\delta _2t\right)\right]=2A_p\mathrm{exp}\left(ip\nu t\right)\mathrm{cos}\left(p\omega t\right),$$
(49)
where

$$\nu \frac{\delta _1+\delta _2}{2}=\frac{\omega _1+\omega _2}{2}\omega _0,\omega \frac{\delta _1\delta _2}{2}=\frac{\omega _1\omega _2}{2}.$$
(50)
(We note incidentally that, in the case of single-frequency forcing, $`\omega _1=\omega _2\omega _e`$ and hence $`\nu =\omega _e\omega _0`$ ($`\omega _e=\omega _0+\nu `$), and $`\omega =0`$.) Hence all our previous analysis is valid for this type of driving, whenever forcing is weak ($`A_p\epsilon ^3`$) and the time scale along which the $`a_p`$’s vary is slow. Specifically this requires $`\nu ,\omega \epsilon ^2`$. If these relations are satisfied, the system will be described by Eq. (43) with $`f`$ given by Eq. (44) and $`\epsilon ^3\alpha _{3,1}`$ given by $`a_{p=1}\left(t\right)`$:

$$_tu(𝐫,t)=\left(\mu \mu _0\right)c_1u+c_2^2u+c_3\left|u\right|^2u+2c_4A_1\mathrm{exp}\left(i\nu t\right)\mathrm{cos}\left(\omega t\right).$$
(51)
In order to simplify this equation let us define

$$U(𝐫,t)=u(𝐫,t)\mathrm{exp}\left(i\nu t\right),$$
(52)
which, substituted into Eq. (51) yields

$$_tU(𝐫,t)=\left[\left(\mu \mu _0\right)c_1i\nu \right]U+c_2^2U+c_3\left|U\right|^2U+2c_4A_1\mathrm{cos}\left(\omega t\right).$$
(53)
The actual state of the system, Eq. (41), is given in terms of $`U`$ as

$$𝐮(𝐫,t)=2Re\left\{U(𝐫,t)\mathrm{exp}\left[i\left(\omega _0+\nu \right)t\right]𝐡\right\}.$$
(54)
Finally, note that, especially in nonlinear optics, one is used to express the state of the system in terms of the negative-frequency part of the oscillations, i.e., as

$$𝐮(𝐫,t)=2Re\left\{\overline{U}(𝐫,t)\mathrm{exp}\left[i\left(\omega _0+\nu \right)t\right]\overline{𝐡}\right\},$$
(55)
and the complex amplitude $`\overline{U}(𝐫,t)`$ verifies the complex conjugate of Eq. (53):

$$_t\overline{U}(𝐫,t)=\left[\left(\mu \mu _0\right)\overline{c}_1+i\nu \right]\overline{U}+\overline{c}_2^2\overline{U}+\overline{c}_3\left|\overline{U}\right|^2\overline{U}+2\overline{c}_4\overline{A}_1\mathrm{cos}\left(\omega t\right).$$
(56)
This absolutely trivial transformation is done here just to point out that the sign that affects the mistuning $`\nu `$ in the forced CGL equation depends on whether one uses the amplitude of the positive- or the negative-frequency part of the oscillations in order to describe the dynamics of the system.


text:
# Deutsch-Jozsa algorithm for continuous variables

## Abstract

We present an idealized quantum continuous variable analog of the Deutsch-Jozsa algorithm which can be implemented on a perfect continuous variable quantum computer. Using the Fourier transformation and XOR gate appropriate for continuous spectra we show that under ideal operation to infinite precision that there is an infinite reduction in number of query calls in this scheme.

In principle, quantum computers can have remarkable computational powers which classical computers cannot DJCV\_sl ; sl1 . In the last few years it has been shown that it is possible for quantum computers to perform certain computational tasks faster than any classical computer pb ; DJCV\_rf ; DJCV\_dd ; DJCV\_dj ; bv ; drs ; div . Quantum computation exploits quantum interference and entanglement to outperform its classical counterparts. The first algorithm promising benefits from quantum parallelism was discovered by Deutsch and Jozsa DJCV\_dj . Soon after this Shor discovered DJCV\_ps  his now famous algorithm for factoring large numbers ej . Subsequently, a fast quantum search algorithm was discovered by Grover DJCV\_lg ; DJCV\_lgr ; in addition, the time-dependent generalization of Grover’s algorithm and its stability under unitary perturbation has been studied DJCV\_akp .

It may be mentioned that all these algorithms are usually designed for qubits. However, in nature there are other classes of quantum systems whose observables form, for example, continuous spectra. Usually a continuous variable can be anything, e.g., position, momentum, energy (unbounded), the amplitudes of the electromagnetic field, etc. It is important to know how these algorithms can be generalized to continuous quantum variables. By learning how, we might make progress towards discovering new algorithms which are perhaps more naturally formulated using continuous variables.

Recently continuous quantum information has played an important role in teleportation sb  and even error-correction codes slb ; sl2  with a possible implementation using linear devices sam . Moreover, quantum computation over continuous variables has also been studied. It was found that the universal continuous variable quantum computation can be effected using simple non-linear operations with coupling provided solely by linear operations DJCV\_lb . Just as standard quantum computation can be thought of as the coherent manipulation of two-level systems qubits, continuous quantum computation can be thought of as the manipulation of qunats.

The first algorithm to have been studied for potential implementation using continuous quantum variables was the Grover virtual database search pbl . Here, we go on to generalize the Deutsch-Jozsa algorithm to continuous variables. This scheme naively gives an infinite speed-up over classical function evaluation.

To start with, let us recall the standard Deutsch-Jozsa algorithm for qubits. In this case we are given a number $`i\{0,\mathrm{},2^n1\}B^n`$ and a “black box” or “oracle query” that computes a binary function $`f(i):B^nB`$. Further, the function $`f(i)`$ which only takes values $`0`$ or $`1`$ is promised  to be either constant or balanced (with an equal number of each type of outcome over all input strings). The aim is to determine this property for $`f(i)`$, i.e., whether it is constant or balanced. On a classical computer in the worst case the oracle query requires $`O(2^n)`$ function evaluations. However, if one calculates the function using reversible quantum operations then only a single function evaluation is required to achieve the goal DJCV\_dj ; cemm .

In the continuous variable setting we pose the problem in the following way. Suppose there is a particle located somewhere along the $`x`$-axis. Since $`x`$ is a continuous variable it can take value from $`\mathrm{}`$ to $`+\mathrm{}`$ (in practice it may be from $`L`$ to $`L`$, where $`L`$ is some length scale involved). Suppose there are two persons Alice and Bob playing a game mike . Alice tells Bob a value of $`x`$ and Bob calculates some function $`f(x)`$ which takes values $`0`$ or $`1`$. Further, Bob has promised Alice that he will use a function which is either constant or balanced. A constant function is $`0`$ or $`1`$ for all values of $`x(\mathrm{},+\mathrm{})`$. For a balanced function, $`f(x)=0`$ or $`1`$ for exactly half of the cases. One can define the balanced function more precisely in the following manner. Imagine that the interval for the continuous variable $`x`$ has been divided into $`n`$ sub-intervals. Let $`\mu `$ be the Lebesgue measure on $``$. A function $`f(x)`$ is balanced provided the Lebesgue measure of the support for where the function is zero is identical to the Lebesgue measure of the support for where the function is one, i.e., $`\mu (\{x|f(x)=0\})=\mu (\{x|f(x)=1\})`$. Now, Alice wants to know whether Bob has chosen a constant or balanced function. In the classical scenario since there are an infinite number of possiblities for $`x`$ Alice needs to ask Bob (who has the oracle) an infinite number of times! However, we can show if we use a perfect continuous variable quantum computer and unitary operators that can be implemented on them, then a single function evaluation is required to know this global property of the function.

Let us consider a continuous variable system whose Hilbert space is infinite dimensional and spanned by a basis state $`|x`$ satisfying the orthogonality condition $`x|x^{}=\delta (xx^{})`$. In a continuous variable scheme a basic operation is the Fourier transformation between position and momentum variables in phase space (analog to the Walsh-Hadamard transformation for qubits). We can define the Fourier transformation as an active operation on a qunat state $`|x`$ as

$$|x=\frac{1}{\sqrt{\pi }}𝑑ye^{2ixy}|y,$$
(1)
where both $`x`$ and $`y`$ are in the position basis. This has been used in developing error correction codes slb ; sam  and Grover’s algorithm for continuous variables pbl . This Fourier transformation can be easily applied in physical situations. For example, when $`|x`$ represents quadrature eigenstate of a mode of the electromagnetic field, $`|x`$ is simply an eigenstate of the conjugate quadrature produced by a $`\pi /2`$ phase delay.

Another useful gate on a continuous variable quantum computer is XOR gate (analogous to the controlled NOT gate for qubits but without the cyclic condition) defined as sam

$$|x|y|x|x+y.$$
(2)
Further, we assume that given a classical circuit for computing $`f(x)`$ there is a quantum circuit which can compute a unitary transformation $`U_f`$ on a continuous variable quantum computer. If a quantum circuit exists that transforms

$$|x|y|x|y+f(x),$$
(3)
then by linearity it can also act on any superposition of qunat states. For example, if we evaluate the function on a state (1) along with another qunat state $`|z`$, we have

$$𝒰_f(|x|z)=\frac{1}{\sqrt{\pi }}𝑑ye^{2ixy}|y|z+f(y).$$
(4)
This shows that using quantum parallelism for idealized qunat computers one can evaluate all possible values of a function simultaneously with one application of $`𝒰_f`$.

Now, we present the Deutsch-Jozsa algorithm for a continuous variable quantum computer. The set of instructions for deciding the constant or balanced nature of function $`f(x)`$ are give below

(i). Alice stores her query in a qunat register prepared in an ideal position eigenstate $`|x_0`$ and attaches another qunat in a position eigenstate $`|\pi /2`$. So the two quants are in the state $`|x_0|\pi /2`$

(ii). She creates superpositions of qunat states by applying the Fourier transformation to the query qunat and the target qunat. The resulting state is given by

$$|x_0|\pi /2=\frac{1}{\pi }𝑑x𝑑ye^{2ix_0x+i\pi y}|x|y.$$
(5)
(iii). Bob evaluates the function using the unitary operator $`𝒰_f`$. The state transforms as

$$\frac{1}{\sqrt{\pi }}𝑑xe^{2ix_0x+i\pi f(x)}|x|\pi /2.$$
(6)
Here, the key role is played by the ancilla qunat state $`|\pi /2`$. To see how the function evaluation takes place consider the intermediate steps given by

$`𝒰_f(|x|\pi /2)={\displaystyle \frac{1}{\sqrt{\pi }}}{\displaystyle 𝑑ye^{i\pi y}𝒰_f(|x|y)}=(1)^{f(x)}|x|\pi /2.`$ (7)
If the function $`f(x)=0`$ there is no sign change and if $`f(x)=1`$ there is a sign change. After the third step performed by Alice, she has a quant state in which the result of Bob’s function evaluation is encoded in the amplitude of the qunat superposition state given in (6). To know the nature of the function she now performs an inverse Fourier transformation on her qunat state.

(iv). The qunat states after Fourier tranform is given by

$$|q=\frac{1}{\pi }𝑑x𝑑x^{}e^{2ix(x_0x^{})}(1)^{f(x)}|x^{}|\pi /2.$$
(8)
(v). Alice measures her qunat by projecting onto the original position eigenstate $`|x_0`$. In an ideal continuous variable scheme the correct projection operator is defined as cohen

$$P_{\mathrm{\Delta }x_0}=_{x_0\mathrm{\Delta }x_0/2}^{x_0+\mathrm{\Delta }x_0/2}𝑑y|yy|.$$
(9)
As has been explained in pbl ; cohen  if the observable has a continuous spectrum then the measurement cannot be performed precisely but must involve some spread $`\mathrm{\Delta }x_0`$. Therefore, the action of projection onto the qunat state after step (iv) is given by

$$P_{\mathrm{\Delta }x_0}|q=\frac{1}{\pi }𝑑x_{x_0\mathrm{\Delta }x_0/2}^{x_0+\mathrm{\Delta }x_0/2}𝑑ye^{2ix(x_0y)}(1)^{f(x)}|y|\pi /2.$$
(10)
Now consider two possibilities. If the function is constant then the above equation reduces to $`\pm |x_0|\pi /2`$. \[In simplifying we need to use the Dirac delta function $`(1/\pi )𝑑xe^{2ix(x_0y)}=\delta (x_0y)`$.\] This means that if Alice measures $`|x_0`$ she is sure that $`f(x)`$ is definitely constant. In the other case, i.e., when the function is balanced she will not get the measurement outcome to be $`|x_0`$. In fact, in the balanced case the outcome is orthogonal to the constant case as the result gives zero. Therefore, a single function evaluation (follwed by a measurement onto $`|x_0`$) in a qunat quantum computer can decide whether the promised function is constant or balanced. Unlike the qubit case, in the idealized continuous variable case the reduction in the number of query calls is from infinity to one.

In conclusion, we have generalised the primitive quantum algorithm (Deutsch-Jozsa algorithm) from the discrete case to the idealized continuous case. It may be worth mentioning that as in error correction codes for continuous-variablees slb , if one replaces the Hadamard transform and XOR gate by their continuous-variable analogs in original Deutsch-Jozsa algorithm for qubit case, then the idealized algorithm works perfectly. This theoretically demonstrates the power of quantum computers to exploit the superposition principle giving an infinite speed up compared to classical scenario. This idealized analysis has not considered the affects of finite precision in measurement or state construction and so whether it may be implemented experimentally remains an open question for further study. Part of the difficulty in extending this work in this direction is that defining an oracle for continuous variables appears to be a difficult task, one that we have carefully avoided here. An alternate way forward might be to consider some sort of “hybrid” approach involving both qunats and qubits. This is precisely what Seth Lloyd considers in the following chapter.

AKP thanks P. van Loock and R. Simon for useful feedback. AKP also thanks G. Giedke for discussions during Benasque Science Center-2000 in Spain on defining balanced function for continuous variables.


text:
# Nondecoupling Effects of Heavy Higgs Particles in Two Higgs Doublet Model 11footnote 1talk at the international conference RG2002, to be published soon in Acta Physica Slovaca

## 1 Introduction

One of the least understood features of the Standard Model (SM) of electroweak interactions is the mechanism of electroweak symmetry breaking descending from the structure of the Higgs sector. Since the measurable quantities usually depend only weakly on its particular realization, many alternative models were proposed.

Perhaps the most popular (nonsupersymmetric) extension of the SM Higgs sector is the so called Two Higgs doublet model (THDM, 2HDM) , a theory with two Higgs doublets instead of one in the usual case. Not only it is capable to reproduce all the predictions of the standard theory but it also provides a nice framework for some possible new phenomena, among other things the CP-violation in the Higgs sector.

Having two doublets in the model the number of Higgs degrees of freedom gets doubled, i.e. we are left with 5 (8 total - 3 Goldstone modes) physical Higgs states. As in the case of SM the mass of the lightest state $`h^0`$ is expected to be at the electroweak scale but the other states $`H^0`$, $`A^0`$, $`H^\pm `$ can in principle be quite massive. Although they cannot be produced in low-energy experiments it is legitimate to ask whether they could contribute to low-energy amplitudes by means of virtual effects i.e. if they decouple in the heavy-mass limit or not.

There is a famous theorem by Appelquist and Carazzone  concerning general decoupling properties of heavy degrees of freedom in field theories. Unfortunately it does not work well in the case of the couplings among the light and heavy sectors growing too fast with the heavy-sector masses. Note that the Higgs couplings are typically proportional to the masses of interacting particles and this spoils the validity of this theorem in many situations involving heavy virtual Higgses coupled to light sector in the game. This is not in general the case of SUSY theories in which the form of the Higgs potential is rather strictly dictated by supersymmetry; from this point of view the heavy Higgs particles in these theories decouple in a usual Appelquist-Carazzone manner, see

The problem of possibly large non-decoupling effects of heavy Higgs particles in THDM was discussed in several papers in 1990s, for instance  and .

As was shown in  we can expect relatively large non-decoupling effects of heavy Higgs bosons in this model for instance in cross-sections of processes involving longitudinal gauge bosons, in particular in $`e^+e^{}W_L^+W_L^{}`$. The magnitude of deviation of this quantity compared to the well-known SM value turns out to be of the order of several percent, which may (at least in principle) be measurable at future facilities.

In the calculation  several simplifications have been made:
1. The ’Equivalence theorem’  used therein works well only in the high-energy limit and therefore should not be used to estimate the non-decoupling features of the model (defined in the low-energy regime).
2. Only the ratio of total cross-sections is given which effectively washes out all the interesting (and probably larger) effects in differential quantities.

From this point of view we find it meaningful to recalculate the ratio of the differential cross-sections of $`e^+e^{}W_L^+W_L^{}`$ between THDM and SM without use of the Equivalence theorem.

## 2 General analysis

Let us define the central quantity of our interest – the ratio of the differential cross-sections of $`e^+e^{}W_L^+W_L^{}`$ in THDM and SM respectively:

$$\delta \frac{\mathrm{d}\sigma ^{\mathrm{𝗍𝗁𝖽𝗆}}(e^+e^{}W^+W^{})}{\mathrm{d}\sigma ^{\mathrm{𝗌𝗆}}(e^+e^{}W^+W^{})}1$$
(1)
Expanding now the THDM amplitude around the well-known SM value it is easy to obtain (at one-loop level)

$$\delta =2\mathrm{R}\mathrm{e}\frac{\mathrm{\Delta }_{\mathrm{𝗍𝗋𝖾𝖾}}+\mathrm{\Delta }_{\mathrm{𝟣}\mathrm{𝗅𝗈𝗈𝗉}}}{_{\mathrm{𝗍𝗋𝖾𝖾}}^{\mathrm{𝗌𝗆}}}+\frac{k_2}{k_1}dk_\gamma \frac{\left|^{\mathrm{𝗍𝗁𝖽𝗆}}\right|^2\left|^{\mathrm{𝗌𝗆}}\right|^2}{\left|_{\mathrm{𝗍𝗋𝖾𝖾}}^{\mathrm{𝗌𝗆}}\right|^2}+\mathrm{}$$
(2)
Here the symbol “$`\mathrm{\Delta }`$” denotes differences of given quantities between THDM and SM, for example $`\mathrm{\Delta }_{\mathrm{𝟣}\mathrm{𝗅𝗈𝗈𝗉}}`$ is the difference of all one-loop contributions to amplitudes between the models; $`^{\mathrm{𝗆𝗈𝖽𝖾𝗅}}`$ denote the corresponding bremsstrahlung amplitudes needed to regulate the IR divergences of $`\mathrm{\Delta }_{\mathrm{𝟣}\mathrm{𝗅𝗈𝗈𝗉}}`$ and $`k_i`$’s are some geometrical factors.

First we can get rid of the bremsstrahlung part of this expression: the IR divergent parts of $`\mathrm{\Delta }_{\mathrm{𝟣}\mathrm{𝗅𝗈𝗈𝗉}}`$ are combined with differences of bremsstrahlung terms to give a perfectly finite quantity which is suppressed by a factor of $`m_e/m_𝗐`$ (Yukawa couplings of Higgses to electrons in the initial state) in comparison with the rest of (2). Next, similar argumentation shows that the same proportionality factor occurs also in $`\mathrm{\Delta }_{\mathrm{𝗍𝗋𝖾𝖾}}`$. Neglecting such terms we are left with

$$\delta 2\mathrm{R}\mathrm{e}\frac{\mathrm{\Delta }_{\mathrm{𝟣}\mathrm{𝗅𝗈𝗈𝗉}}^{\mathrm{𝗂𝗋}\mathrm{𝖿𝗂𝗇}.}}{_{\mathrm{𝗍𝗋𝖾𝖾}}^{\mathrm{𝗌𝗆}}}$$
(3)
All we need are therefore contributions of IR-finite one-loop graphs which are not common to both models.

## 3 Two Higgs doublet model (THDM)

Let us now specify the basic features of THDM in more detail. As we already know the presence of the second doublet gives rise to 5 Higgs states in the spectrum: neutral scalars $`h^0`$ and $`H^0`$, charged scalars $`H^\pm `$ and a neutral pseudoscalar $`A^0`$. The most general form of the Higgs potential (in terms of $`SU(2)`$ doublets $`\mathrm{\Phi }_1`$, $`\mathrm{\Phi }_2`$, )

$$\begin{array}{c}V(\mathrm{\Phi }_1,\mathrm{\Phi }_2)=m_{11}^2\mathrm{\Phi }_1^{}\mathrm{\Phi }_1+m_{22}^2\mathrm{\Phi }_2^{}\mathrm{\Phi }_2m_{12}^2\mathrm{\Phi }_1^{}\mathrm{\Phi }_2m_{12}^2\mathrm{\Phi }_2^{}\mathrm{\Phi }_1+\hfill \\ \hfill +\frac{\lambda _1}{2}(\mathrm{\Phi }_1^{}\mathrm{\Phi }_1)^2+\frac{\lambda _2}{2}(\mathrm{\Phi }_2^{}\mathrm{\Phi }_2)^2+\lambda _3(\mathrm{\Phi }_1^{}\mathrm{\Phi }_1)(\mathrm{\Phi }_2^{}\mathrm{\Phi }_2)+\\ \hfill +\lambda _4(\mathrm{\Phi }_1^{}\mathrm{\Phi }_2)(\mathrm{\Phi }_2^{}\mathrm{\Phi }_1)+\frac{\lambda _5}{2}(\mathrm{\Phi }_1^{}\mathrm{\Phi }_2)^2+\frac{\lambda _5^{}}{2}(\mathrm{\Phi }_2^{}\mathrm{\Phi }_1)^2+\\ \hfill +[\lambda _6(\mathrm{\Phi }_1^{}\mathrm{\Phi }_1)+\lambda _7(\mathrm{\Phi }_2^{}\mathrm{\Phi }_2)](\mathrm{\Phi }_1^{}\mathrm{\Phi }_2)+[\lambda _6^{}(\mathrm{\Phi }_1^{}\mathrm{\Phi }_1)+\lambda _7^{}(\mathrm{\Phi }_2^{}\mathrm{\Phi }_2)](\mathrm{\Phi }_2^{}\mathrm{\Phi }_1)\end{array}$$
(4)
is less restrictive than in SUSY theories; there is enough freedom for the coupling constants $`\lambda _i`$ and the ’mass-parameters’ $`m_{ij}`$ to give rise to some new phenomena like the above mentioned CP-violation in the Higgs sector etc. Moreover, the bounds on the Higgs-mass patern are not so stringent as for example in the SUSY theories, it is not a problem to have the masses of $`A^0`$ or $`H^\pm `$ at a scale of several TeV ,.

In general there are two basic types of THDM concerning the mass generation of up- and down-types of fermions. In type-I models both the up and down fermion masses are generated by one of the doublets only in analogy with the minimal Higgs model while in type-II theories one of the doublets generates the up-type and the second one the down-type masses in a similar way as in the MSSM. However, our analysis turns out to be model-independent because there are no relevant Yukawa couplings in the one-loop leading term computation.

## 4 Calculation of $`\mathrm{\Delta }_{\mathrm{𝟣}\mathrm{𝗅𝗈𝗈𝗉}}^{\mathrm{𝗂𝗋}\mathrm{𝖿𝗂𝗇}.}`$

Although the full set of Feynman diagrams contributing to $`\mathrm{\Delta }_{\mathrm{𝟣}\mathrm{𝗅𝗈𝗈𝗉}}^{\mathrm{𝗂𝗋}\mathrm{𝖿𝗂𝗇}.}`$ is quite large, many of them can be safely neglected. It is mainly due to the presence of a suppressing Yukawa factor in all the diagrams involving Higgs couplings to electrons as indicated above.

Since we are using the on-shell renormalization scheme there is no need for renormalization of external legs. However, the vertex and propagator counterterms become nontrivial being not only the “$`1/\epsilon `$” parts of the dimensionally-regularised expressions but fixed by on-shell renormalization conditions . These additional structures need their own detailed discussion.

With these observations in mind we can classify all the one-loop topologies contributing to $`\mathrm{\Delta }_{\mathrm{𝟣}\mathrm{𝗅𝗈𝗈𝗉}}^{\mathrm{𝗂𝗋}\mathrm{𝖿𝗂𝗇}.}`$: it turns out that the only really important types of graphs are the following:

(5)
Here the dark shaded blob corresponds to loops involving at least one Higgs boson while the lighter stands for loops without any Higgs inside, all of them including finite parts of the relevant counterterms.

### 4.1 Oblique-type corrections

Although there is no direct suppression by the Yukawa factors in first two graphs in (5) we can neglect them because of the one-loop mixed propagators whose contributions in general look like $`\mathrm{\Pi }(p^2,m_i^2)(p_1+p_2)_\alpha `$. Contracting this expression with the leptonic current and using the Dirac equation one reproduces again suppressing factor $`m_e/m_𝗐`$. The case of the third topology is not so clear but it can be shown to exhibit the decoupling properties in the heavy (physical) Higgs mass limit . (This can be easily seen in the particular case of the on-shell renormalization scheme .)

### 4.2 Vertex corrections – one-loop TGV differences

Thus we are left with only the fourth topology in (5). To proceed, we need the differences of one-loop renormalized triple gauge vertex structures $`\gamma W^+W^{}`$ and $`ZW^+W^{}`$ ; let us denote them by $`\mathrm{\Delta }\mathrm{\Gamma }_{\sigma \mu \nu }^\gamma `$ and $`\mathrm{\Delta }\mathrm{\Gamma }_{\sigma \mu \nu }^𝖹`$ respectively. Relevant graphs can be divided into several clusters (charged bosons are denoted by a generic symbol $`A^\pm `$ (for example $`G^\pm `$ stands for the charged Goldtone bosons in $`R_\xi `$ gauge), while the neutral ones by $`B`$,$`C`$):

(6)
Their UV-divergences should be cancelled by counterterms descending from Ward-identity $`\delta Z_𝗐=\delta Z_g`$ connecting the $`W`$-boson wavefunction renormalization constant and the gauge coupling renormalization constant in the scheme fixed as in .

Having everything at hand we can write the leading contribution to $`\mathrm{\Delta }_{\mathrm{𝟣}\mathrm{𝗅𝗈𝗈𝗉}}^{\mathrm{𝗂𝗋}\mathrm{𝖿𝗂𝗇}.}`$ in the form

$$\mathrm{\Delta }_{\mathrm{𝟣}\mathrm{𝗅𝗈𝗈𝗉}}^{\mathrm{𝗂𝗋}\mathrm{𝖿𝗂𝗇}.}\underset{𝖵=\gamma ,𝖹}{}\overline{v}(p_1)\gamma _\lambda u(p_2)g_{ee𝗏}\frac{ig^{\lambda \sigma }}{sm_𝖵^2}g_{\mathrm{𝗏𝗐𝗐}}\mathrm{\Delta }\mathrm{\Gamma }_{\sigma \mu \nu }^𝖵\epsilon ^\mu (q_1)\epsilon ^\nu (q_2)$$
(7)
## 5 Results and conclusions

Since we are dealing with many complicated diagrams (there are 46 graphs in an $`R_{\xi =1}`$ gauge in (6) and 9 others to calculate the finite parts of on-shell counterterms), we are forced to use a computer. We have utilized Mathematica 4.0 with Feyncalc and Looptools. The figures Fig.1 and Fig.2 correspond to ratios of differential cross-sections of particular initial state helicity configurations $`e_L^+e_R^{}W_L^+W_L^{}`$ (in which $`\delta `$ is expected to be largest). Moreover, in this particular case the leading contribution to $`\delta `$ turns out to be $`\mathrm{cos}\theta ^{}`$-independent (CMS scattering angle). Although we use a slightly different Higgs mass pattern to exhibit mainly the basic features of $`\delta `$, its magnitude is in rough agreement with the expectation of . (The discontinuities in derivatives of the second curve originate from the fact, that the loop integrals in (3) acquire non-zero imaginary parts above some values of $`\sqrt{s}`$ which correspond to tresholds of productions of the loop particles in on-shell final states.)

As we see, at least in some regions of the parametric space we can expect relatively large non-decoupling (note the slope of the plot at Fig.1) effects of heavy Higgs bosons in the considered quantity; they can reach the order of several percent.
Acknowledgement
I would like to thank to Prof. Jiří Hořejší for all his comments and helpful suggestions during the preparation of this contribution. This work was supported by “Centre for Particle Physics”, project No. LN00A006 of the Ministry of Education of the Czech Republic.


text:
# A Near-Infrared Survey for Galactic Wolf-Rayet Stars

## 1. Introduction

Optical surveys within our Galaxy are severely hampered by dust obscuration; therefore complete samples must be obtained with longer wavelength observations. Here we describe a survey in the Galactic plane at K-band wavelengths, where the extinction is much lower than traditional V-band surveys. Our scientific driver is the discovery of the young stellar population in our Galaxy through the detection of evolved massive stars. These stars have strong emission lines, which makes them relatively easy to detect using narrow-band filters.

As a massive O star evolves, its spectrum becomes dominated by emission lines, arising either in a dense stellar wind, or in circumstellar material produced by mass loss. Stellar emission lines are most pronounced in Wolf-Rayet (WR) stars which have lifetimes $`<`$ 10 Myr, and thus are excellent tracers of recent star formation, and so also Galactic structure. A complete sample of Galactic WR stars will also enable us to understand the distribution of WN and WC sub-types in a high metallicity environment. Additionally, WR stars are critical components in our quest to understand how star formation proceeds. For example, most of the previously known WRs are relatively isolated, but recent searches in the IR have found a plethora of these objects in stellar clusters near the Galactic center (e.g. Figer et al. 1999).

## 2. Survey Description

We employ four filters situated on prominent lines of He I (2.06 um), C IV (2.08 um), Br $`\gamma `$ (2.166 um), He II (2.19 um), and three continuum filters covering wavelengths bluer and redder than each of the line filters (2.03, 2.14, and 2.248 um) (Blum & Damineli 1999, Homeier et al. 2000). In this way we minimize the effects of differential reddening on the continuum of each individual star. Each frame is analyzed with the set of routines collectively called the ’DoPhot’ package, a psf-fitting photometry package optimized for large images done in a pipeline process. DoPhot identifies, classifies, and performs photometry on objects in an image, making successive passes over the image to search each time for fainter objects. Candidates are selected for spectroscopic follow-up based on relative line and continuum photometry and image subtraction with visual examination of the output (Homeier et al. 2002, in prep).

In 1996 we began taking data at the 1.5m at CTIO, using the Cerro Tololo Infrared Imager (CIRIM). This data comprises approximately 1.5 sq degrees of sky near the Galactic Center, half of which we have follow-up spectroscopic observations for. With this data we have discovered two new WR stars and two objects with surrounding nebular emission.

## 3. Newly Discovered Emission-line Objects

Both of the newly discovered WRs appear to be late-type WC stars (WC7-8). One is $`8`$ arcmin from the Quintuplet cluster (WR1), or 20 pc projected distance. We speculate that it has been ejected from the Quintuplet or is a former member of an as-yet undetected massive stellar cluster in the Galactic Center region. The second WR is part of the Quintuplet (WR2, Quintuplet source), and is located just outside the region examined by Figer et al. (1999). Spectra of these two objects are presented in Figure 1, and coordinates are listed in Table 1. Finding charts will be published in an upcoming paper (Homeier et al. 2002, in prep).

Of the two objects with nebular emission, one is surrounded by a spherical nebula with strong H I lines, including Br $`\gamma `$ 2.166 $`\mu `$m, and He I 2.06 $`\mu `$m emission. This object could be a PN, LBV, or a young O star ionizing the surrounding gas. It is located on the edge of a region with extremely large extinction, and there appears to be more nebular emission a few arcseconds to the south-east. We speculate that this is a star-forming region on the edge of a molecular cloud, and further observations are needed to deduce the distance of this region and whether it contains high-mass stars.

## References

Blum, R., & Damineli, A. 1999, in K. A. van der Hucht, G. Koenigsberger, & P. R. J. Eenens (eds.), Wolf-Rayet Phenomena in Massive Stars and Starburst Galaxies, Proc. IAU 193 (San Francisco: ASP), 472

Figer, D., McLean, I., & Morris, M. 1999, ApJ, 514, 202

Homeier, N., Blum, R., Conti, P., Damineli, A. 2000 AAS, 197.0412

Homeier, N., Blum, R., Conti, P., Damineli, A., Pasquali, A. 2002, in prep for A&A