text:
# Relativistic Massive Vector Condensation

## I Introduction

The Quantum chromodynamics (QCD) phase diagram as a function of temperature, chemical potential and the number of light flavors is very rich and highly structured. QCD should behave as a color superconductor for a sufficiently large quark chemical potential . This leads to new phenomenological applications associated with the description of quark stars, neutron star interiors and the physics near the core of collapsing stars . Much is known about the phase structure of QCD at nonzero temperature through a combination of perturbation theory and lattice simulations while the phase structure at nonzero chemical potential and for large numbers of flavors has been less extensively explored . Standard importance sampling methods employed in lattice simulations fail at nonzero chemical potential for $`N_c=3`$ since the fermionic determinant is complex. For $`N_c=2`$, the situation is very different since the quarks are in a pseudoreal representation of the gauge group and lattice simulations can be performed .

At non zero chemical potential Lorentz invariance is explicitly broken down to the rotational subgroup $`SO(3)`$ and higher spin fields can condense, thus potentially enriching the phase diagram structure of QCD and QCD-like theories. Indeed in  it has been suggested that some of the lightest massive vectors at non zero baryon chemical potential may condense. Recent lattice simulations seem to support these predictions . In the context of the Electroweak theory vector condensation in the presence of a strong external magnetic field has also been studied in .

In this paper we explore the condensation of relativistic massive vectors when introducing a non zero chemical potential via an effective Lagrangian approach. In order for the vectors to couple to the chemical potential they must transform under a given global symmetry group. In view of the possible physical applications we consider the vectors to belong to the adjoint representation of the $`SU(2)`$ group.

To set the stage for the more complicated case of the vector field we first consider a scalar field in the adjoint representation of $`SU(2)`$ at non zero chemical potential. We choose the chemical potential to lie along one of the generator’s directions explicitly breaking $`SU(2)`$ to $`U(1)\times Z_2`$. We then investigate the dispersion relations of the 3 fields as function of the chemical potential and show the existence of a second order phase transition when the chemical potential equals the zero density mass of the scalar field. In the broken phase we observe the emergence of one gapless mode associated to the spontaneous breaking of the $`U(1)`$ symmetry possessing linear dispersion relations. This is consistent with the Chadha-Nielsen  counting rule for Goldstone’s bosons when Lorentz invariance is not assumed.

The situation becomes more involved when turning to the spin-one massive field case. When the vector field condenses the rotational $`SO(3)`$ invariance breaks spontaneously together with part of the flavor (i.e. global) symmetry of the theory. The condensate locks together space and flavor symmetries. Following the standard mean-field approach we identify our order parameter with the vector field itself.

In order to elucidate the mechanism we construct the simplest effective Lagrangian at zero chemical potential. This consists of a standard kinetic type term and a tree level potential constructed as a series expansion in the order parameter. We truncate our potential at the fourth order in the fields while retaining all the terms allowed by Lorentz and $`SU(2)`$ symmetry. The introduction of the potential term allows us to describe the phase transition in some detail while our results will be general.

We prove that the number of gapless states emerging when the flavor$`\times `$rotational symmetry breaks spontaneously is insensitive to the choice of the coefficients in the effective potential. However the momentum dependence of the associated dispersion relations is heavily affected by the choice of the potential. The latter is related to the differences between the symmetries of the potential with respect to the ones of the kinetic type terms. The full vector Lagrangian (i.e. kinetic $``$ potential) always possesses the $`Z_2\times U(1)\times SO(3)`$ (flavor$`\times `$rotational) symmetry which at high enough chemical potential breaks spontaneously to $`Z_2\times SO(2)`$. The gapless modes are linked to the 3 broken generators. One state is an $`SO(2)`$ scalar while the other two states constitute a vector of $`SO(2)`$. The scalar Goldstone possesses linear dispersion relations independently of the choice of the vector potential and it resembles the gapless mode of the scalar theory. However the $`SO(2)`$ vector state can have either quadratic or linear dispersion relations. In particular (as summarized in the Table I) the quadratic ones occur when the potential parameters are such that it possesses an enhanced $`SO(6)`$ symmetry. When condensing the vector breaks the $`SO(6)`$ symmetry of the potential down to $`SO(5)`$ while the kinetic-type term still has only a $`U(1)\times SO(3)`$ invariance. In this case the potential has 5 flat directions (i.e. 5 null curvatures) and we would count in the spectrum 2-gapless vectors and an $`SO(2)`$ scalar. However the reduced symmetry of the kinetic term prevents the emergence of two independent gapless vectors while turning the linear dispersion relations of one vector state into quadratic ones. This result is also in agreement with the Chadha-Nielsen counting scheme. However since the Goldstone modes with quadratic dispersion relations must be counted twice relatively to the ones with linear dispersion relations we saturate the number of broken generators associated with the symmetries of just the potential term which is larger than the ones associated with the symmetries of the full Lagrangian. It is certainly important to investigate, in the future, if the choice of the parameter space leading to an $`SO(6)`$ symmetry for the potential term is stable against quantum corrections.

Our model can be related immediately to the physics of 2 colors QCD which has been partially analyzed in  and is within reach of present lattice studies. Here we suggest that vector condensation can be the leading mechanism for a superfluid transition for 2 color QCD with one flavor. Indeed in this case the global symmetry is $`SU(2)`$ which remains intact at zero chemical potential and no Goldstone bosons emerge. The massive spectrum of the theory still contains an $`SU(2)`$ spin-one massive multiplet. When turning on a non zero chemical baryon potential the vector condensation can lead to a novel superfluid phase transition. Clearly our model can also be applied to QCD with 3 color at non zero isospin chemical potential .

The paper is structured as follows. In section II we first consider a toy model for a scalar field in the regular representation of $`SU(2)`$ at non zero chemical potential and investigate all of the relevant features of the model. We then propose a simple effective Lagrangian for vectors which allows us to uncover the main features related to the condensation of relativistic massive spin-one fields. We conclude and point out some physical applications in section III. In particular we suggest that new phases should set in for QCD like theories. Finally we extend our simple model to $`D1`$ number of space dimensions (with Euclidean $`SO(D)`$ Lorentz symmetry).

## II Vector Condensation at Nonzero Chemical Potential

### A Toy Model: Scalar Theory

In order to illustrate and disentangle the main problems related to the vector condensation we first consider the case of a scalar field in the adjoint representation of an $`SU(2)`$ global symmetry. This will allow us to show the feature which only depends on the global symmetries first. The toy model Lagrangian is:

$$=\mathrm{Tr}\left[_\mu \varphi ^\mu \varphi \right]m^2\mathrm{Tr}\left[\varphi ^2\right]\lambda \left(\mathrm{Tr}\left[\varphi ^2\right]\right)^2,$$
(1)
with $`\varphi =\varphi ^aT^a`$ and $`a=1,2,3`$ and $`T^a=\tau ^a/2`$. $`\varphi ^a`$ is a real field and $`\tau ^a`$ the standard Pauli matrices.

We introduce the chemical potential associated with the $`T^3`$ generator. We shall see that our choice turns out to be particularly relevant when considering 2 colors QCD and one flavor at non zero quark chemical potential. The effects of the chemical potential are included by generalizing the covariant derivative as follow:

$$D_\mu \varphi =_\mu \varphi i\mu \delta _{\mu 0}[T^3,\varphi ].$$
(2)
$`\mu `$ is the chemical potential. Expanding the kinetic term we have that the Lagrangians reads in components $`\varphi =(\phi ^\pm =(\varphi ^1i\varphi ^2)/\sqrt{2},\psi =\varphi ^3)`$:

$``$ $`=`$ $`{\displaystyle \frac{1}{2}}(\psi )^2+|\phi |^2+\mu \phi ^{(+)}i_0\phi ^{(+)}\mu \phi ^{()}i_0\phi ^{()}\left(m^2\mu ^2\right)|\phi |^2{\displaystyle \frac{m^2}{2}}\psi ^2{\displaystyle \frac{\lambda }{4}}\left(|\phi |^2+\psi ^2\right)^2.`$ (3)
Our choice of the chemical potential explicitly breaks the $`SU(2)`$ ’flavor’-symmetry to a $`U(1)`$ acting on the complex field $`\varphi `$ and a discrete $`Z_2`$ symmetry for the neutral field $`\psi `$. For $`\mu m`$, the model has a unique vacuum at $`|\phi |=\psi =0`$ and the effect of the chemical potential is simply to split the energies of the charged states as follows:

$`E_{\phi ^\pm }(\stackrel{}{p})=\mu +\sqrt{m^2+\stackrel{}{p}^2},E_\psi (\stackrel{}{p})=\sqrt{m^2+\stackrel{}{p}^2}.`$ (4)
The mass-gaps are defined as the energy of each state evaluated at zero momentum. Expanding the dispersion relations for small momenta gives:

$`E_{\phi ^\pm }(\stackrel{}{p})=(m\mu )+{\displaystyle \frac{\stackrel{}{p}^2}{2m}}+𝒪(\stackrel{}{p}^4).`$ (5)
This leads to the presence of a gapless excitation with a nonrelativistic type of dispersion relation at $`\mu =m`$. Note that due to the loss of Lorentz–covariance the mass-gaps do not coincide with the zero density masses. The latter are defined as the curvature of the potential evaluated on the vacuum and are $`m_{\phi ^\pm }^2=m^2\mu ^2`$, $`m_\psi ^2=m^2`$. Since the $`\psi `$ field is not affected at all by the chemical potential it does not condense and its mass-gap is $`m_\psi =m`$ for $`\mu m`$.

For $`\mu >m`$ the classical potential is unstable when $`\varphi =0`$ and we choose the global minimum for

$$\phi ^1=\sqrt{\frac{\mu ^2m^2}{\lambda }},$$
(6)
breaking the $`U(1)\times Z_2`$ symmetry to $`Z_2`$. Computing the curvature of the potential on the new vacuum for all the fields yields:

$`m_\psi ^2=\mu ^2,m_{\varphi ^1}^2=2\left(\mu ^2m^2\right),m_{\varphi ^2}^2=0.`$ (7)
Note that in this phase the $`\psi `$ mass is given precisely by the chemical potential. This is so since $`\psi `$ is not charged under the chemical potential. After diagonalizing the quadratic terms we obtain the following dispersion relations:

$`E_\pm (\stackrel{}{p})=\sqrt{\stackrel{}{p}^2+\overline{\mu }^2\pm \sqrt{\overline{\mu }^4+4\stackrel{}{p}^2\mu ^2}},`$ (9)
with $`\overline{\mu }^2=3\mu ^2m^2`$. Note that with our choice of the vev $`\phi ^\pm `$ are no longer eigenstates. We see that a gapless mode persists and the massive one has the mass gap $`\sqrt{2}\overline{\mu }`$. At small momenta the dispersion relation of the gapless (Goldstone-) mode is

$`E_{}^2(\stackrel{}{p}){\displaystyle \frac{\mu ^2m^2}{\overline{\mu }^2}}|\stackrel{}{p}|^2+{\displaystyle \frac{2\mu ^4}{\overline{\mu }^6}}|\stackrel{}{p}|^4.`$ (11)
Hence in the broken phase we obtain the standard superfluid Goldstone boson with an energy depending linearly on the momentum, $`E_{}v_s(\mu )|\stackrel{}{p}|`$. However the velocity $`v_s(\mu )`$ vanishes at the transition point $`v_s(\mu =m)=0`$ turning the dispersion relation of the gapless excitation into quadratic ones (i.e. $`E_{}|\stackrel{}{p}|^2/2m`$). This behavior can be better understood by rewriting the velocity as function of $`m_{\varphi ^1}^2`$ as follows:

$$v_s^2=\frac{m_{\varphi ^1}^2}{m_{\varphi ^1}^2+4\mu ^2}.$$
(12)
We observe that the superfluid velocity depends on the curvature of the potential in the direction orthogonal to the gapless excitation. At the second order phase transition point this curvature is zero, reflecting the conformal invariance of the potential for the $`\phi `$–fields. In the absence of Lorentz breaking we have two massless modes at this point both with linear dispersion relations. At nonzero chemical potential only one mode is gapless while the information about the conformal invariance is now inherited by the quadratic dispersion law. We shall see in the next section, that this behavior is quite general, i.e. whenever the potential has more symmetries than the derivative terms, this information is encoded nontrivially into the momentum dependence of the dispersion relations. Furthermore we can now smoothly connect the dispersion relations of each state when traversing the different regions of the phase diagram. Our result is in agreement with the Chada-Nielsen  counting scheme as well as with recent studies related to the physics of color superconductivity .

### B The Vector Story

There are different ways to introduce vector fields at the level of the effective Lagrangian (for example the hidden local gauge symmetry of Ref. , or the antisymmetric tensor field of Ref. ) and they are all equivalent at tree-level. We choose to introduce the massive vector fields following the method outlined in . This method also allows a straightforward generalization of our model to an arbitrary number of space dimensions. We adopt the following model Lagrangian for describing relativistic spin one fields (in $`3+1`$ dimensions) belonging to the regular representation of $`SU(2)`$

$``$ $`=`$ $`{\displaystyle \frac{1}{4}}F_{\mu \nu }^aF^{a\mu \nu }+{\displaystyle \frac{m^2}{2}}A_\mu ^aA^{a\mu }{\displaystyle \frac{\lambda }{4}}\left(A_\mu ^aA^{a\mu }\right)^2+{\displaystyle \frac{\lambda ^{}}{4}}\left(A_\mu ^aA^{a\nu }\right)^2,`$ (13)
with $`F_{\mu \nu }^a=_\mu A_\nu ^a_\nu A_\mu ^a`$, $`a=1,2,3`$ and metric convention $`\eta ^{\mu \nu }=\mathrm{diag}(+,,,)`$. Here $`m^2`$ is the tree level mass term and $`\lambda `$ and $`\lambda ^{}`$ are positive dimensionless coefficients with $`\lambda >\lambda ^{}`$. Other possible choices of the parameters do not guarantee stability of the potential and will not be considered in the following. This Lagrangian is not renormalizible, contrary to the scalar theory presented above. However it possesses the same global symmetries as the toy model.

The effect of a nonzero chemical potential associated to a given conserved charge - related to the generator (say $`B`$) - can be readily included  by modifying the derivatives acting on the vector fields:

$$_\nu A_\rho _\nu A_\rho i[B_\nu ,A_\rho ],$$
(14)
with $`B_\nu =\mu \delta _{\nu 0}BV_\nu B`$ where $`V=(\mu ,\stackrel{}{0})`$. The vector kinetic term modifies according to:

$`\mathrm{Tr}\left[F_{\rho \nu }F^{\rho \nu }\right]`$ $``$ $`\mathrm{Tr}\left[F_{\rho \nu }F^{\rho \nu }\right]4i\mathrm{Tr}\left[F_{\rho \nu }[B^\rho ,A^\nu ]\right]2\mathrm{T}\mathrm{r}\left[[B_\rho ,A_\nu ][B^\rho ,A^\nu ][B_\rho ,A_\nu ][B^\nu ,A^\rho ]\right].`$ (15)
The terms due to the kinetic term, after integration by parts, yields

$`_{kinetic}`$ $`=`$ $`{\displaystyle \frac{1}{2}}A_\rho ^a\left\{\delta _{ab}\left[g^{\rho \nu }\mathrm{}^\rho ^\nu \right]4i\gamma _{ab}\left[g^{\rho \nu }V{\displaystyle \frac{V^\rho ^\nu +V^\nu ^\rho }{2}}\right]+2\chi _{ab}\left[VVg^{\rho \nu }V^\rho V^\nu \right]\right\}A_\nu ^b`$ (16)
with

$$\gamma _{ab}=\mathrm{Tr}\left[T^a[B,T^b]\right],\chi _{ab}=\mathrm{Tr}\left[[B,T^a][B,T^b]\right].$$
(17)
For $`B=T^3`$ we have

$$\gamma _{ab}=\frac{i}{2}ϵ^{ab3},\chi _{11}=\chi _{22}=\frac{1}{2},\chi _{33}=0.$$
(18)
The chemical potential induces a “magnetic-type” mass term for the vectors at tree-level. The similarities with respect to the scalar toy model as well as the symmetries of the vect Lagrangian are more easily understood using the following Euclidean notation:

$$\phi _M^a=(A_M^1,A_M^2),\psi _M=A_M^3,$$
(19)
with $`A_M=(iA_0,\stackrel{}{A})`$ and metric signature ($`+,+,+,+`$). In these variables the potential reads:

$`V_{Vector}`$ $`=`$ $`{\displaystyle \frac{m^2}{2}}\left[|\stackrel{}{\phi }_0|^2+\psi _M^2\right]+{\displaystyle \frac{m^2\mu ^2}{2}}|\stackrel{}{\phi }_I|^2+{\displaystyle \frac{\lambda }{4}}\left[|\stackrel{}{\phi }_M|^2+\psi _M^2\right]^2{\displaystyle \frac{\lambda ^{}}{4}}\left[\stackrel{}{\phi }_M\stackrel{}{\phi }_N+\psi _M\psi _N\right]^2`$ (20)
with $`I=X,Y,Z`$ while $`M,N=0,X,Y,Z`$ and repeated indices are summed over. At zero chemical potential $`V_{Vector}`$ is invariant under the $`SO(4)`$ Lorentz transformations while only the $`SO(3)`$ symmetry is manifest at non zero $`\mu `$. Clearly the $`SU(2)`$ symmetry is also explicitly broken, at non zero chemical potential, to $`U(1)`$ as for the scalar model.

For the reader’s convenience we summarize the quadratic term in the fields Lagrangian in the new variables.

$`_{kinetic}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\psi _M\left[\delta _{MN}\mathrm{}_E_M_N\right]\psi _N+{\displaystyle \frac{1}{2}}\phi _M^a\left[\delta _{MN}\mathrm{}_E_M_N\right]\phi _N^a\phi _M^a\epsilon ^{ab3}\left[\delta _{MN}V{\displaystyle \frac{V_M_N+V_N_M}{2}}\right]\phi _N^b`$ (21)
$``$ $`\phi _M^a\left[VV\delta _{MN}V_MV_N\right]\phi _N^a`$ (22)
According to the value of the chemical potential we distinguish two phases.

#### 1 The Symmetric Phase: $`\mu m`$

Here the $`SO(4)`$ Lorentz and $`SU(2)`$ symmetries are explicitly broken to $`SO(3)`$ and $`U(1)`$ respectively by the presence of the chemical potential but no condensation happens (i.e. $`<\stackrel{}{\phi }_M>=<\psi _M>=0`$). The masses at zero chemical potential are:

$`M_{\phi _0^a}^2=M_{\psi _M}^2=m^2,M_{\phi _I^a}^2=m^2\mu ^2.`$ (23)
The dispersion relations obtained by diagonalizing the 12 (4 space-time $`\times `$ 3 $`SU(2)`$ states) by 12 quadratic eigenvalue problem leads to 3 physical vectors (i.e. each of the following states has 3 components) with the following dispersion relations:

$`E_\phi ^{}=\pm \mu +\sqrt{\stackrel{}{p}^2+m^2},E_\psi `$ $`=`$ $`\sqrt{\stackrel{}{p}^2+m^2}.`$ (24)
This shows that when approaching $`\mu =m`$ the 3 physical components associated with $`E_\phi ^{}`$ will become massless signaling an instability. Indeed we now show that for higher values of the chemical potential a vector type condensation sets in.

#### 2 The Spin-Flavor Broken Phase: $`\mu >M`$

In this phase the global minimum of the potential is for $`<\phi _0^a>=<\psi _M>=0`$, while we can choose the vev to lie in the (Spin-Flavor) direction:

$$<\phi _X^1>=\sqrt{\frac{\mu ^2m^2}{\lambda \lambda ^{}}}.$$
(25)
We have a manifold of equivalent vacua which are obtained rotating the chosen one under a $`Z_2\times U(1)\times SO(3)`$ transformation. The choice of the vacuum partially locks together the Lorentz group and the internal symmetry while leaving unbroken only the subgroup $`Z_2\times SO(2)`$. Two generators associated to the Lorentz rotations are now spontaneously broken together with the $`U(1)`$ generator.

To compute the full dispersion relations of the theory we first need to provide the curvatures evaluated on the vacuum:

$`\psi \mathrm{Sector}`$ (26)
$`M_{\psi _0}^2`$ $`=`$ $`M_{\psi _Y}^2=M_{\psi _Z}^2=m^2+\lambda {\displaystyle \frac{\mu ^2m^2}{\lambda \lambda ^{}}},M_{\psi _X}^2=\mu ^2,`$ (27)
$`\phi ^1\mathrm{Sector}`$ (28)
$`M_{\phi _0^1}^2`$ $`=`$ $`\mu ^2,M_{\phi _X^1}^2=2(\mu ^2m^2),M_{\phi _Y^1}^2=M_{\phi _Z^1}^2=0,`$ (30)
$`\phi ^2\mathrm{Sector}`$ (31)
$`M_{\phi _0^2}^2`$ $`=`$ $`m^2+\lambda {\displaystyle \frac{\mu ^2m^2}{\lambda \lambda ^{}}},M_{\phi _X^2}^2=0,M_{\phi _Y^2}^2=M_{\phi _Z^2}^2=\lambda ^{}{\displaystyle \frac{\mu ^2m^2}{\lambda \lambda ^{}}}.`$ (32)
In general three states have null curvature (specifically $`M_{\phi _X^2}^2=M_{\phi _Y^1}^2=M_{\phi _Z^1}^2=0`$) however for the special case $`\lambda ^{}=0`$ we have 5 zero curvature states. To explain this behavior we not that for $`\lambda ^{}=0`$ the potential possesses $`SO(6)`$ global symmetry which breaks to $`SO(5)`$ when the vector field condenses. The associated 5 states would correspond to the ordinary Goldstone modes in the absence of an explicit Lorentz breaking. Using these curvatures we compute the dispersion relations. Four of the $`\phi `$ states have (2 per each sign) the following dispersion relations:

$`E_{\phi _V^{}}^2`$ $`=`$ $`\stackrel{}{p}^2+\mathrm{\Delta }^2\sqrt{4\mu ^2\stackrel{}{p}^2+\mathrm{\Delta }^4}`$ (33)
with

$$\mathrm{\Delta }^2=2\mu ^2+\frac{\lambda ^{}}{2}\frac{(\mu ^2m^2)}{\lambda \lambda ^{}},$$
(34)
and to each sign we associate a vector (with two components) with respect to $`SO(2)`$. In the limit of small momenta:

$`E_{\phi _V^{}}^2`$ $`=`$ $`{\displaystyle \frac{\lambda ^{}}{2\mathrm{\Delta }^2}}{\displaystyle \frac{\mu ^2m^2}{(\lambda \lambda ^{})}}\stackrel{}{p}^2+{\displaystyle \frac{2\mu ^4}{\mathrm{\Delta }^6}}|\stackrel{}{p}|^4+𝒪(p^6)=v_{\phi _V^{}}^2\stackrel{}{p}^2+\mathrm{},`$ (35)
$`E_{\phi _V^+}^2`$ $`=`$ $`2\mathrm{\Delta }^2+\left[1+{\displaystyle \frac{2\mu ^2}{\mathrm{\Delta }^2}}\right]\stackrel{}{p}^2{\displaystyle \frac{2\mu ^4}{\mathrm{\Delta }^6}}|\stackrel{}{p}|^4+𝒪(p^6).`$ (36)
For the last two $`\phi `$ physical modes (scalars of $`SO(2)`$) we show directly the dispersion relations as a momentum expansion.

$`E_{\phi _S^{}}^2`$ $`=`$ $`{\displaystyle \frac{\mu ^2m^2}{3\mu ^2m^2}}\stackrel{}{p}^2+𝒪(p^4)=v_{\phi _S^{}}^2\stackrel{}{p}^2+\mathrm{},E_{\phi _S^+}^2=2(3\mu ^2m^2)+\gamma \stackrel{}{p}^2+𝒪(p^4),`$ (37)
where the subscripts $`V,S`$ stand for vector and scalar respectively, and

$$\gamma =\frac{(3+8c)\mu ^6+(716c)m^2\mu ^4+2(5c4)m^4\mu ^22(c1)m^6}{\mu ^2(3\mu ^2m^2)(c\mu ^2(c1)m^2)}>0,$$
(38)
and $`c=\lambda /(\lambda \lambda ^{})`$. Finally for completeness we show the three $`\psi `$ type state dispersion relations:

$$E_{\psi _I}=\sqrt{\stackrel{}{p}^2+M_{\psi _I}^2},$$
(39)
with $`M_{\psi _I}`$ defined in eq. (27). As for the scalar toy model, the velocities of our gapless modes can be expressed directly in terms of the curvatures in the directions orthogonal to the gapless modes as follows:

$$v_{\phi _V^{}}^2=\frac{M_{\phi _Y^2}^2}{M_{\phi _Y^2}^2+4\mu ^2},v_{\phi _S^{}}^2=\frac{M_{\phi _X^1}^2}{M_{\phi _X^1}^2+4\mu ^2}.$$
(40)
We find that the number of gapless states emerging when the flavor$`\times `$rotational symmetry spontaneously breaks is insensitive to the choice of the coefficients in the effective potential. However the momentum dependence depends crucially on the symmetries of the potential. The full vector Lagrangian (i.e. kinetic $``$ potential) possesses the $`Z_2\times U(1)\times SO(3)`$ (flavor$`\times `$rotational) symmetry which breaks, at non zero chemical potential, when the vector condenses to $`Z_2\times SO(2)`$. The gapless modes are associated to the 3 broken generators. One state is a scalar of $`SO(2)`$ while the other two fields constitute a vector of $`SO(2)`$. The scalar Goldstone possesses linear dispersion relations independently of the choice of the vector potential. This is the Goldstone associated to the $`U(1)`$ breaking of the previously investigated scalar toy model. However for the dispersion relations of the $`SO(2)`$ vector state we have either quadratic or linear dispersion relations. The quadratic ones occur when the potential parameters are such that it possesses an enhanced $`SO(6)`$ symmetry. When condensing the vector breaks the potential symmetries $`SO(6)`$ to $`SO(5)`$ while the kinetic term (i.e. the derivative part of the Lagrangian) is still invariant under the continuous $`U(1)\times SO(3)`$ only. In this case the potential has 5 flat directions (i.e. 5 null curvatures (occurring for $`\lambda ^{}=0`$)) and we would count in the spectrum 2-massless vectors and a scalar of $`SO(2)`$. However the reduced symmetry of the kinetic term prevents the emergence of another gapless vector type mode while turning the linear dispersion relations of one of the vector states into quadratic ones. As for the toy model case, at the second order phase transition point, all the velocities of the gapless modes vanish regardless of the value of the coupling constants. Again the information of the symmetries (conformal and/or global) of the potential, at nonzero chemical potential, are nicely transferred from the curvature to the momentum dependence of the dispersion relations via the velocities of the gapless modes. Our result is in agreement with the Chadha-Nielsen counting scheme.

## III Conclusions and Physical Applications

We summarized our findings in the Table I and conclude our work by considering some possible physical applications of massive relativistic vector condensation at high chemical potential.

### A 2 Color QCD

Two color QCD at non zero chemical potential recently attracted a great flurry of interest since lattice simulations can be reliably performed at non zero quark chemical potential. In  vector condensation at non zero quark chemical potential has been predicted for two color QCD with an even numbers of flavors. The predictions are supported  by lattice studies. Here we confirmed the previous results and provided a more detailed study of the vector condensation phenomenon at non zero chemical potential.

An interesting direct application of our results is for QCD with two colors and one light flavor. This theory has global quantum symmetry group $`SU(2)`$. The extra classical axial $`U_A(1)`$ symmetry is broken by the Adler-Bell-Jackiw anomaly. At zero chemical potential Lorentz invariance cannot be broken and the simplest bilinear condensate (see Ref.  for conventions) is of the type:

$$ϵ^{c_1c_2}ϵ^{\alpha \beta }Q_{\alpha ,c_1}^IQ_{\beta ,c_2}^JE_{IJ},$$
(41)
with $`E=2iT^2`$ the antisymmetric matrix in the flavor space and $`Q_{\alpha c}^I`$ is a Weyl spinor with $`\alpha =1,2`$ the spinorial index, $`c=1,2`$ the color and $`I=1,2`$ the flavor indices. The generators of $`SU(2)`$ in the fundamental representation are $`T^a=\tau ^a/2`$ with $`\tau ^a`$ the standard Pauli matrices. It is easy to check that $`SU(2)`$ remains unbroken since we have $`T_{}^{a}{}_{}{}^{T}E+ET^a=0`$ for all of the $`SU(2)`$ generators. So we have no Goldstone bosons at all in the theory. However we do have physical vectors (i.e. spin 1 fields) and massive scalar states at zero density. These states should be comparable in mass (recall the $`\eta ^{}`$ versus the octet of ordinary vector bosons in 3 color QCD). The massive spin one fields (i.e. the hadrons) must transform according to the regular representation of $`SU(2)`$ and carry non trivial baryon number which is proportional to $`T^3`$.

Our analysis can be immediately used to predict that when turning on a non zero baryon chemical potential a superfluid phase transition will set in for $`\mu `$ larger then the mass of the vector fields. This is so since we always have at least one gapless mode with linear dispersion relations. We also predict the spontaneous breaking of rotational invariance and the emergence of a vector state (i.e. doublet of $`SO(2)`$) with either linear or quadratic dispersion relations. It is straightforward to extend and apply our analysis to the case of different number of flavors and to ordinary QCD at finite isospin chemical potential. .

### B Extra-dimensions

Finally we extend our model to an $`SU(2)`$ charged vector in $`D1`$ number of space dimensions (with $`SO(D)`$ Euclidean Lorentz symmetry) with an associated non zero chemical potential. In this case the chemical potential first breaks explicitly $`SO(D)`$ to the spatial $`SO(D1)`$. At sufficiently high chemical potential $`SO(D1)`$ breaks spontaneously to $`SO(D2)`$. The gapless states consist of a scalar field under the $`SO(D2)`$ transformations with linear dispersion relations and one vector of $`SO(D2)`$ with either linear or quadratic dispersion relations.

###### Acknowledgements.

It is a pleasure for us to thank M.P. Lombardo, H.B. Nielsen, and P. Olesen for helpful discussions. We are also indebt to J. Schechter for discussions and careful reading of our manuscript. The work of F.S. is supported by the Marie–Curie fellowship, while W.S. acknowledges support by DAAD.


format: text

text:
# Universal Attractors of Reversible Aggregate-Reorganization Processes

## Abstract

We analyze a general class of reversible aggregate-reorganization processes. These processes are shown to exhibit globally attracting equilibrium distributions, which are universal, i.e. identical for large classes of models. Furthermore, the analysis implies that for studies of equilibrium properties of any such process, computationally expensive reorganization dynamics such as random walks can be replaced by more efficient, yet simpler methods. As a particular application, our results explain the recent observation of the formation of similar fractal aggregates from different initial structures by diffusive reorganization \[Filoche and Sapoval, Phys. Rev. Lett. 85, 5118 (2000)\].

Many complex systems like microbial colonies, colloidal aggregates, or particle-adsorbate structures are composed of a large number of similar units Bunde\_et\_al . For the growths of such particle clusters, the model of diffusion-limited aggregation (DLA) has been introduced Witten  and intensively investigated in several generalizations Bunde\_et\_al ; Meakin88 ; Erzan . In these non-equilibrium models, particles are permanently injected into the system and stick irreversibly to an existing cluster which grows towards a fractal structure. In real systems, however, particles are often rearranged, dynamically transforming the aggregate structure. This rearrangement is appropriately modelled by processes exhibiting reversible dynamics Bunde\_et\_al ; Meakin88 ; FilochePRL .

Recently, Filoche and Sapoval have proposed a specific model of diffusion-reorganized aggregation (DRA) FilochePRL  in which particles detach from the boundary of a given cluster of fixed size and reaggregate after a random walk on the underlying square lattice. In numerical simulations, they observe the formation of fractal structures similar to those found in various irreversible models. The DRA process transforms different initial aggregates, e.g. simple one-dimensional and two-dimensional structures, into statistically similar fractal aggregates after many reorganization steps.

In this Letter, we analyze a general class of reversible processes that dynamically reorganize a fixed number of particles constituting one connected aggregate on some lattice. Particles are allowed to reorganize by any reversible process that retains a connected aggregate, including diffusion as a particular case. We prove that all these processes exhibit equilibrium distributions that are globally attracting. Our general mathematical perspective reveals that such attractors are universal, i.e. independent of the dynamical model details: In discrete-time models, this distribution only depends on the mechanism of particle disaggregation but is independent of the subsequent particle repositioning dynamics. In particular, this implies that ordinary diffusion and different mechanisms like surface diffusion Bunde\_et\_al  or Lévy-flights Levy  are equivalent in this context. Furthermore, we show that in continuous-time models equilibria are independent of both the mechanism of repositioning and that of disaggregation. On a $`d`$-dimensional cubic lattice, if all connected aggregates are accessible by the repositioning dynamics, this leads to the uniform distribution that is known as the ensemble of lattice animals (cf. Janse ; BotetAndWessel ; Lubensky\_et\_al ), currently also investigated in number theory under the name polyominoes (Golomb\_et\_al  and refs. therein).

In addition to these analytical results, we present strong numerical evidence that, for sufficiently large aggregates, there is no significant difference between equilibrium distributions obtained from discrete- or continuous-time models. Our results are valid for a general class of processes, including reorganization of aggregates consisting of different kinds of particles as well as repositioning of whole clusters instead of single particles. Here, for simplicity of presentation, we focus on processes in which single, identical particles are reorganized.

We consider reorganization processes of aggregates – certain connected collections of a fixed number $`N`$ of particles on some lattice. Starting from an aggregate $`A`$, one step of the discrete-time dynamics is defined by the following two rules (Fig. 1):

* Randomly select a particle from a set of free particles $`F(A)`$ with equal probability, where $`F(A)`$ is chosen such that the aggregate remains connected.
* Reposition this particle by a reversible process at some site $`x`$ of possible aggregation sites $`D(A,B)`$ of an intermediate object. By reversibility, this particle is again a free particle of the new aggregate $`B`$.

The above two rules are repeatedly applied to the evolving aggregate. Note that physically the step-by-step application of the two rules reflects a separation of time scales that is realized in many processes Bunde\_et\_al : The disaggregation of a particle occurs much slower than its repositioning. To model a specific reorganization process, the sets $`F(A)`$ of free particles (i), the repositioning (ii) as well as the underlying lattice can be defined as desired as long as every aggregate under consideration can be reached iterating the above two rules.

To be specific we focus on $`d`$-dimensional cubic lattices, but the analysis can be analogously applied for arbitrary lattices. We fix aggregates on a lattice of $`(N+1)^d`$ sites with periodic boundary conditions and view translated aggregates of the same shape as different states. A given process is identified as a Markov chain Norris  with state space $`S_N`$ consisting of all connected aggregates of $`N`$ particles, which have at least one free particle, i.e. non-empty $`F(A)`$. The state of this discrete-time process is then described by the distribution $`p_t(A)`$, giving the probability that after $`t`$ reorganization steps the aggregate $`A`$ is present. Starting with initial distribution $`p_0(A)`$, a transition matrix $`T`$ specifies the time evolution

$$𝐩_t=T𝐩_{t1}=T^t𝐩_0$$
(1)
of the state vector $`𝐩_t:=(p_t(A_1),\mathrm{},p_t(A_{|S_N|}))^\mathrm{T}`$. As $`S_N`$ is finite, the long-time ($`t\mathrm{}`$) behaviour of the discrete-time Markov chain is determined by a steady state distribution $`𝐩_d^{}=T𝐩_d^{}`$. If the chain is irreducible, i.e. every aggregate $`AS_N`$ can be reached by the process from every other, $`𝐩_d^{}`$ is the global attractor towards which every initial distribution converges Norris . This uniqueness property depends on the choice of free particles $`F(A)`$. In general, one reasonably restricts the state space $`S_N`$ to an appropriate subset of aggregates consistent with $`F(A)`$. For instance, if $`F(A)`$ is the set of all one-bond particles BotetAndWessel , only tree-like aggregates are considered. In such cases our analysis applies to the respective irreducible subset of $`S_N`$ (also simply denoted by $`S_N`$ in the following). With this convention, the stationary distribution $`𝐩_d^{}`$ is the unique attractor for every initial distribution $`𝐩_0`$ on $`S_N`$. We now specify the transition matrix $`T`$ which completely determines the dynamics of the process.

If a one-step transition from aggregate $`A`$ to $`B`$ is possible, denoted by $`AB`$, the particle which is disaggregated and its reaggregation site are determined uniquely aa . Thus the underlying lattice does not explicitly enter the analysis. The transition probabilities are given by the matrix elements

$$T_{BA}=\{\begin{array}{cc}\frac{1}{|F(A)|}r_{BA}\hfill & \text{if}AB\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(2)
which, if non-zero, consist of two factors: According to rule (i) a free particle is randomly selected from $`F(A)`$ with uniform probability $`1/|F(A)|`$. Then, the particle is repositioned at the intermediate aggregate with probability $`r_{BA}`$ to obtain aggregate $`B`$ according to rule (ii). The reversibility of the repositioning procedure is characterized by the symmetry $`r_{BA}=r_{AB}`$. Thus, the dynamics is described by the master equation

$`p_{t+1}(A)p_t(A)`$ $`=`$ $`{\displaystyle \underset{BS_N}{}}\left[T_{AB}p_t(B)T_{BA}p_t(A)\right]`$ (3)
$`=`$ $`{\displaystyle \underset{BS_N}{}}\left[{\displaystyle \frac{p_t(B)}{|F(B)|}}{\displaystyle \frac{p_t(A)}{|F(A)|}}\right]r_{AB}.`$
A steady state solution of this equation is obviously given by the equilibrium distribution

$$p_d^{}(A)=|F(A)|/\underset{BS_N}{}|F(B)|,$$
(4)
which fulfills the condition of detailed balance, i.e. $`\frac{p_d^{}(B)}{|F(B)|}=\frac{p_d^{}(A)}{|F(A)|}`$ for all $`A,BS_N.`$ From the above we know that $`𝐩_d^{}`$ is the unique attractor of the dynamics (3). Hence, in equilibrium the probability to observe an aggregate $`A`$ is proportional to its number $`|F(A)|`$ of free particles.

It seems natural to realize repositioning by a random walk of the disaggregated particle, because it models diffusion in real systems FilochePRL ; BotetAndWessel . This random walk determines the transition probabilities $`r_{AB}=r_{BA}`$ between aggregates $`A`$ and $`B`$. However, independent of the actual values of $`r_{AB}`$, and hence the repositioning procedure (ii), the process leads to the same equilibrium distribution $`𝐩_d^{}`$. Thus, to study the equilibrium properties of the reorganization process, one can replace the random walk in numerical simulations: For instance, after selection of a particle, it may just be repositioned at every aggregation site with equal probability $`r_{AB}=1/|D(A,B)|`$, where $`D(A,B)`$ is the set of possible aggregation sites of the intermediate aggregate (cf. Fig. 1 and ref. Janse ). This straightforward procedure is much simpler to implement than random walks and significantly reduces the computational effort. Furthermore, one expects a faster convergence to the equilibrium distribution Guruswami : While for diffusive repositioning the particle re-aggregates with high probability near or even at its disaggregation site, the alternative process leads to a fast spreading of particles so that the state space is likely to be sampled faster. Moreover, diffusing particles in $`d2`$ dimensions have an infinite average return time leading to computational complications FilochePRL ; BotetAndWessel  that can be avoided using alternative processes.

Up to now, we have considered discrete-time evolution rules, which arise when performing computer simulations. However, for real reorganization processes, continuous-time modelling, where the free particles $`F(A)`$ disaggregate independently at a rate $`\gamma `$, often seems more appropriate. This leads to a continuous-time Markov chain with completely symmetric transition rates $`\gamma r_{AB}`$ occurring in the master equation

$$\frac{d}{dt}p_t(A)=\underset{BS_N}{}\gamma \left[p_t(B)p_t(A)\right]r_{AB}.$$
(5)
Thus, the attractor of these continuous-time processes is

$$p_c^{}(A)=1/|S_N|.$$
(6)
This is a universal equilibrium distribution which is independent of all dynamical model details. It is important to note that, although $`F(A)`$ does not enter (6), it implicitly determines the set $`S_N`$ of accessible aggregates. In the case of $`d`$-dimensional cubic lattices and for every choice of $`F(A)`$ for which the process can reach all possible connected aggregates of $`N`$ particles (cf. our simulations below), this uniform distribution $`𝐩_c^{}`$ is known as the ensemble of lattice animals Lubensky\_et\_al ; Golomb\_et\_al .

The question arises, how relevant the differences between discrete- and continuous-time modelling are. To answer this question, we consider the fraction

$$f_N(A)=|F(A)|/N$$
(7)
of free particles of an aggregate $`A`$ of $`N`$ particles. This observable determines the difference in the equilibrium distributions of discrete-time and their associated continuous-time models because it directly enters (4) in the form $`p_d^{}(A)=f_N(A)/_Bf_N(B)`$. In equilibrium, the average fractions of free particles obviously satisfy the inequality

$`f_N_{𝐩_d^{}}={\displaystyle \underset{A}{}}f_N(A)p_d^{}(A)={\displaystyle \frac{_Af_N^2(A)}{_Bf_N(B)}}`$ $`f_N_{𝐩_c^{}}`$ (8)
because $`f_N_{𝐩_c^{}}=_Bf_N(B)/|S_N|`$ by definition. Nevertheless, Monte-Carlo simulations of a specific process simulations  show that the averages $`f_N`$ clearly converge towards the same value for discrete- and continuous-time models (Fig. 2).

Differences become insignificant already for moderately sized aggregates.

This is quantified by considering the widths of the equilibrium distributions $`\sigma (f_N):=\left(f_Nf_N\right)^2^{\frac{1}{2}}`$ (Fig. 2, inset). We find a power law decay $`\sigma (f_N)N^\alpha `$ with $`\alpha =0.50\pm 0.02`$ for both time scales as would be expected for large random structures. This has important consequences: Calculations similar to the above (8) show immediately that

$$g_N_{𝐩_d^{}}=g_N_{𝐩_c^{}}\left(1+𝒪(N^\alpha )\right)$$
(9)
for every observable $`g_N(A)`$. Thus, for sufficiently large aggregates, measurements of arbitrary observables, and, equivalently, the equilibrium distributions, in discrete- and continuous-time models cannot be distinguished. This reflects the fact that the number $`|F(A)|`$ of free particles of almost all aggregates $`AS_N`$ is close to its uniform ensemble average $`|F(A)|`$ for sufficiently large $`N`$. Since this fact is of purely configurational origin, we expect a similar convergence (9) for other models, also in $`d>2`$ dimensions.

As one important observable, we have studied the fractal dimensions of aggregates from both processes. Sample aggregates as well as box-counting dimensions $`D_N`$ are shown in Fig. 3. Even for moderately sized aggregates, no significant differences have been observed. With increasing aggregate sizes, the $`D_N`$ approach from below scaling  the well-estimated fractal dimension $`D1.56`$ Lubensky\_et\_al ; Golomb\_et\_al  of square-lattice animals.

The above results predict that the uniform distribution among all connected aggregates which can be reached by a specific process is the globally attracting, unique equilibrium distribution for a wide class of models. In particular, our results explain the recent observation that fractal aggregates are reached from arbitrary initial structures by a certain diffusive reorganization process FilochePRL . Nevertheless, our results strongly indicate that the equilibrium aggregates from that study are in the same universality class as lattice animals and thus have a fractal dimension $`D1.56`$ and not, as suggested, $`D=1.74\pm 0.02`$ FilochePRL . A possible explanation of this discrepancy is that simulations in that study did not reach equilibrium at the time of measurement, despite a huge number of simulation steps. In fact, in our simulations simulations  aggregates of $`N=10000`$ particles only reached equilibrium after about $`2\times 10^8`$ particle rearrangements.

It is important to note, that in certain systems such a slow convergence may imply that the equilibrium state is not reached during the time of real or numerical experiments (cf. FilochePRL ). Thus although an equilibrium distribution is guaranteed to be reached, the long-term behavior of such systems may be governed by transient states with long lifetimes and non-universal statistical properties.

In summary, we have shown that large classes of reversible aggregate-reorganization processes possess universal attractors which are independent of the model details but may in general depend on whether discrete- or continuous-time modelling is used. Furthermore, our numerical investigations indicate that statistical properties are essentially identical for both time scales. These results imply that commonly used diffusion modelling can be replaced by computationally fast methods which in addition lead to a faster convergence towards the (identical) equilibrium distribution. In general, one specific process of aggregate-reorganization can be investigated in a number of different alternative models. Which alternative is used can be freely chosen depending on which one is more convenient to study analytically or by computer simulations.

In particular, our results explain the recently discovered “statistical attractors” in diffusive reorganization processes FilochePRL  but predict a different final equilibrium state. More generally, the results can serve as a basis for investigating equilibrium statistical properties of a general class of aggregate-reorganization processes. Reversible reorganization is always characterized by the symmetry $`r_{AB}=r_{BA}`$, which, in turn, has been the main requirement for universality. In particular, aggregates consisting of different kinds $`\alpha `$ of particles disaggregating with particle-specific rates $`\gamma _\alpha `$ are also described because a transition $`AB`$ determines the disaggregating particle uniquely. Moreover, the theory may be generalized to reorganization of whole sub-clusters rather than just single particles (cf. also Meakin88 ; Janse ). Hence more complex models of aggregate-reorganization can also be captured.

We have imposed only a small number of requirements on the dynamics which appear to be natural assumptions in many physical situations, in particular reversibility of repositioning and separation of time-scales. Thus our main results might not only be of importance to aggregate-reorganization processes. We hope that the basic ideas can also find applications in studies of other complex systems exhibiting reversible reorganization.

We thank M. Biehl, D. Brockmann, T. Geisel, L. Hufnagel, T. Kottos, and M. Prähofer for useful discussions. This work has been supported by the Max Planck Society.


text:
# Exact Persistence Exponent for One-dimensional Potts Models with Parallel Dynamics

## Abstract

We obtain $`\theta _p(q)=2\theta _s(q)`$ for one-dimensional $`q`$-state ferromagnetic Potts models evolving under parallel dynamics at zero temperature from an initially disordered state, where $`\theta _p(q)`$ is the persistence exponent for parallel dynamics and $`\theta _s(q)=\frac{1}{8}+\frac{2}{\pi ^2}[cos^1\{(2q)/q\sqrt{2}\}]^2`$ \[PRL, 75, 751, (1995)\], the persistence exponent under serial dynamics. This result is a consequence of an exact, albeit non-trivial, mapping of the evolution of configurations of Potts spins under parallel dynamics to the dynamics of two decoupled reaction diffusion systems.

Consider a stochastic variable, such as the position of a Brownian particle $`x(t)`$ moving in one dimension as a function of time ($`t`$) or the state $`\sigma (t)`$ ($`+1`$ or $`1`$) of an Ising spin evolving under a prescribed probabilistic dynamics. The following question arises naturally in many physical contexts: What is the probability that some attribute of the system, such as the state of the Ising spin or, for the Brownian particle, the sign of its $`x`$-coordinate, is unchanged in an interval $`[0,t]`$? Such questions relate to the “persistence” properties of stochastic processes and have attracted both theoretical and experimental attention in recent years.

This Letter studies persistence, under parallel dynamics, in one-dimensional ferromagnetic $`q`$-state Potts models with nearest neighbor interactions, described by the Hamiltonian

$$H=\underset{i}{}\delta _{\sigma _i,\sigma _{i+1}},$$
(1)
where $`\sigma _i`$, the Potts spin at site $`i`$, takes values $`0\mathrm{}(q1)`$. In terms of the conventional (serial) zero-temperature dynamics of this system, the problem of persistence is posed in the following way: Start from initial configurations in which each spin takes any of the allowed $`q`$ values with equal probability. Then relax the system by allowing one spin at a time to lower its energy by aligning itself with one of its neighbours, chosen from right or left with equal probability. Such a dynamics is serial (equivalently, sequential) since a single spin update is performed at each step. Now fix any spin and ask: What is the probability, averaged over initial configurations, that this spin has not changed its state up to time $`t`$? This quantity, the persistence probability $`P(t)`$, was first found numerically to decay as a power law,

$$P(t)\frac{1}{t^{\theta _s(q)}},$$
(2)
with $`\theta _s(q)`$ a new, $`q`$dependent non-trivial persistence exponent; the subscript denotes serial dynamics. A subsequent analytic tour de force derived $`\theta _s(q)`$ exactly,

$$\theta _s(q)=\frac{1}{8}+\frac{2}{\pi ^2}[cos^1\{(2q)/q\sqrt{2}\}]^2.$$
(3)
The persistence exponent is argued to be a new exponent characterizing the dynamics.

Barring this remarkable result, we know of very few other exact results in the persistence problem for many interacting variables. Further, these results, besides being restricted to one dimensional problems, have so far, with the exception of Ref., been confined to serial dynamics. The issue of the dynamics is important here, since the persistence exponent does not share the universal character of the static and dynamic exponents. It is thus of obvious interest to find other exact results for persistence in interacting many-body systems, to understand the origins of such scale invariant behaviour.

This Letter obtains the following result connecting the persistence exponent under parallel dynamics $`\theta _p(q)`$ for one-dimensional $`q`$-state ferromagnetic Potts models with the corresponding exponent $`\theta _s(q)`$ for serial dynamics:

$$\theta _p(q)=2\theta _s(q),$$
(4)
where $`P(t)1/t^{\theta _p(q)}`$ in the parallel case. We demonstrate the validity of this result first numerically and then provide a proof via a sequence of mappings to the dynamics of the following reaction diffusion system:

$`A+A`$ $`\varphi `$ $`\mathrm{Probability}1/(q1)`$ (5)
$`A`$ $`\mathrm{Probability}(q2)/(q1)`$ (6)
$`B+B`$ $`\varphi `$ $`\mathrm{Probability}1/(q1)`$ (7)
$`B`$ $`\mathrm{Probability}(q2)/(q1),`$ (8)
describing particles of two species, labelled $`A`$ and $`B`$, which diffuse and interact in one dimension. Apart from the doubling of species, this reaction scheme is identical to the one used in a domain wall mapping for persistence in the Potts model with serial dynamics, with the particles representing the domain walls. The initial distribution of $`A`$ and $`B`$ particles can be assigned given any initial configuration of Potts spins through the mapping we describe below. However, as we show here, the $`A`$ and $`B`$ particles which enter the analysis of the case of parallel dynamics are not the same as the domain walls in the serial case.

In terms of this mapping, whether a given site is persistent or not at time $`t`$ is given by the probability that it is not crossed by either a particle of type A or B upto time $`t`$. Since the two reactions are independent, provided there are no significant correlations in the initial placement of A and B particles, the joint probability that a given site is persistent with respect to the motion of both A and B particles simply factors into the product of independent probabilities at long times, as

$$P(t)\frac{1}{t^{\theta _s(q)}}\frac{1}{t^{\theta _s(q)}}\frac{1}{t^{2\theta _s(q)}}\frac{1}{t^{\theta _p(q)}},$$
(9)
yielding immediately the result of Eq. 4.

We first present results from our numerical simulations of this problem. Each spin can take one of $`(q1)`$ values, which we take, for concreteness, to be represented by the integers $`\{0\mathrm{}(q1)\}`$. The zero temperature dynamics evolves a spin configuration $`\{\sigma (t)\}`$ at time $`t`$ to a configuration $`\{\sigma (t+1)\}`$ at time $`t+1`$ through the following simple rule: Each spin at time $`t+1`$ assumes the value of one of its neighboring spins at time $`t`$, chosen from right or left with equal probability. Each such step in time constitutes a single Monte Carlo step. The parallel nature of the dynamics follows from the fact that all spins are updated together in a single time step.

We have simulated parallel dynamics using the above rules on Potts systems with different $`q`$ values, on lattices of linear size $`L=10^5`$ sites and for times $`t10^4`$, applying periodic boundary conditions. We average over a fairly large number of initial conditions, typically $`10^210^3`$, starting from configurations in which each spin is independently assigned a value in the range $`\{0\mathrm{}(q1)\}`$ with equal probability. We compute the standard persistence probability $`P(t)`$, defined as the probability that the spin at a given site has not changed its state up to time $`t`$, averaged over all sites and over an ensemble of initial conditions.

Fig. 1 shows the persistence probability $`P(t)`$ for ferromagnetic Potts systems, evolving under parallel dynamics, for $`q`$ = 3, 5, 8 and 20. Given Eqs. 3 and 4, we expect a power law tail to $`P(t)`$ with exponents $`\theta _p(3)1.076,\theta _p(5)1.386,\theta _p(8)1.588`$ and $`\theta _p(20)1.821`$. As can be seen from the fits displayed, the data are consistent with our proposal.

In a recent paper, coauthored with Shukla, we studied the Ising version of the problem of persistence with parallel dynamics, obtaining the result of Eq. 4 specialized to $`q=2`$ i.e. $`\theta _p(2)=3/4=2\theta _s(2)=2(3/8)`$. We briefly sketch the mapping which enabled us to show this exactly and then describe how such a mapping is extended to solve the more general Potts case.

For the one-dimensional Ising model with serial dynamics, domain walls obey the reaction diffusion scheme $`A+A0`$. Whether a site is persistent or not at some late time $`t`$, is equivalent to asking whether it has been crossed by a domain wall at an earlier time. As shown in Ref, the diffusing particles whose motion and annihilation govern the persistence properties of the Ising model with parallel dynamics are not the domain walls of the serial case, but entities which we termed as zero-field sites i.e. sites in which the spin is aligned with either one, but not both, of its neighbours.

Given any configuration of the Ising model, all sites can be classified (uniquely) as stable (S), unstable (U) or zero-field (Z) depending on whether they will flip with probability $`0,1`$ or $`1/2`$ respectively in the next time step. While pairs of permanently bound zero-field sites can be identified with the domain walls of the serial case, the two sites constituting such a pair can move independently under parallel dynamics. A region of stable sites is always bounded by zero-field sites; a site embedded in an initially stable region is persistent upto time $`t`$ if it has not been crossed by a zero-field site in the interval $`[0,t]`$. The doubling of the persistence exponent is a consequence of the fact that zero-field sites on A and B sublattices decouple from each other under the dynamics. Identifying zero-field sites on the $`A`$ (or $`B`$) sublattice as particles labelled $`A`$ or $`B`$ respectively, the following equations describe the interaction of these particles: $`A+A0`$ and $`B+B0`$; these equations can be obtained from Eq. 8 upon setting $`q`$ to 2. As shown in Ref., zero-field sites on the A(B) sublattice at time $`t`$ hop to sites on the B(A) sublattice at time $`t+1`$. The $`A`$ and $`B`$ labels are identified with respect to an arbitrary initial time, chosen for convenience to be $`t=0`$.

It is tempting to conjecture that a similar sublattice decoupling is responsible for the doubling of the persistence exponent with parallel dynamics vis. a vis the serial case, for general $`q`$. This idea, as we will show, is essentially correct, but for the following caveat: The A and B particles are not simply the zero-field sites of the Potts model. This result follows from the observation that the classification of all spins as stable, unstable or zero-field, as outlined in the Ising case, is incomplete when applied to the problem of persistence in Potts models with $`q>2`$.

Arguments against a naive generalization of ideas from the Ising case follow from inspecting the pictures of Potts configurations evolving under parallel dynamics shown in Figs. 2(a) and 2(b). In these pictures, obtained from simulations of a 3-state Potts model, the configuration corresponding to the earliest time is shown on the top row, while each subsequent row represents the configuration at the following time step. We have identified and labelled zero-field sites in terms of the $`A`$ and $`B`$ particles described above. Note that while A(B) particles hop at each time step to sites on the B(A) sublattice in most of the updates shown, this strict independance of the two sublattices is violated in the penultimate row of Fig. 2(a). Evidently, the dynamics of the zero-field sites alone, unlike the case for the one-dimensional Ising model, does not decouple the two sublattices in the more general Potts case. A related difficulty is apparent from inspection of the configurations of Fig. 2(b). Note that an $`AB`$ pair appears to be created upon evolution. Such a reaction cannot be described by Eq. 8, which yields a strictly non-increasing number of $`A`$ and $`B`$ particles.

Our mapping to the reaction diffusion process described by Eq. 8 in terms of configurations of the $`q`$state Potts model proceeds as follows: For the purposes of illustration, we consider a 3-state Potts model in which the states are labelled by 0, 1 and 2. If the spin at site $`i`$ is in a particular Potts state $`q`$ while its neighbours are in the same Potts state, for example the central spin in the sequence $`\mathrm{}000\mathrm{}`$, that site is labelled a stable ($`S`$) site. If a given spin at site $`i`$ has one spin (on either side) in the same state, while the other spin is in a different state, for example the central spin in $`\mathrm{}002\mathrm{}`$, we will term such a spin as a zero-field ($`Z`$) site. As in the Ising case, we will label these zero-field sites as A or B, depending on the sublattice to which they can be assigned. Note that a region of stable sites is bordered by zero-field sites.

The feature which makes this problem non-trivially different from the Ising case is the following: We must distinguish two categories of unstable sites: unstable sites of type I, labelled U1, have both neighbours in the same state which, however, differs from the state of the central spin. For example, the central spin in the configuration $`\mathrm{}202\mathrm{}`$, is an unstable spin of type I. Finally, we will define unstable sites of type II, labelled U2, in which the spins neighbouring a central spin belong to Potts states which differ from each other as well as from the state of the central spin, the central spin in the configuration $`\mathrm{}012\mathrm{}`$ being an example. Note that unstable sites of type II are absent in the Ising model.

We have now provided a prescription for mapping, given the rules above, any configuration of Potts spins to configurations of S,Z, U1 and U2. In the Ising case, U2 sites being absent, we proceeded to identify the zero-field sites ($`Z`$) as particles of type A or B and went on to show that this mapping yielded precisely Eq. 8 above in the $`q=2`$ limit, as a consequence of the decoupling of the two sublattices. The central idea which allows us to extend this intuition to the Potts model is the following: We relabel unstable sites of type II, i.e. sites labelled U2, as A or B particles, depending on the sublattice to which they belong at $`t=0`$. Given this relabelling, it is now a trivial, though tedious, exercise to verify that the rules for the interaction of A and B particles are precisely those defined by Eq. 8 and that the decoupling of the dynamics of A and B particles follows exactly. The chief result of this paper, our Eq. 4, is then a natural consequence of the arguments which led to Eq. 9.

To make this remarkable simplification explicit, we have taken the configurations of Figs. 2(a) and (b) and added the new $`A`$ and $`B`$ particles obtained by appropriately labelling unstable sites of type II as described above. These configurations are shown in Figs 3(a) and (b). Note that the incorporation of these new sites immediately solves both problems with the earlier mapping (Figs. 2(a) and (b)) which used only the zero-field sites: The configuration of Fig. 2(a) which showed an A particle remaining on the same sublattice at the following time step is now interpreted in terms of the hopping of both the $`A`$ particle and an adjacent $`B`$ particle associated with an unstable site of type II. The difficulty of the apparent creation of an $`AB`$ pair (Fig. 2(b)), is removed by recognizing that this simply arises due to the interconversion of a pair of unstable sites of type II and a pair of zero-field sites.

We now demonstrate, for a simple case, how the rates of Eq. 8 are obtained as a consequence of this mapping. Consider the time evolution of the following spin configuration in a Potts model under parallel dynamics: $`aaaabcccc`$, where Potts states $`a`$ and $`c`$ form stable domains. If $`b`$ is distinct from both $`a`$ and $`c`$ (with $`ac`$), then the mapping described above leads to the following description of the configuration: $`\mathrm{}\mathrm{}.ABA\mathrm{}..`$, where the first zero-field site has been taken to belong, arbitrarily, to sublattice $`A`$. Note that the two $`A`$ particles are associated with zero-field sites, whereas the $`B`$ particle is associated with an unstable site of type II. The following are the possible configurations at the next time step $`aaaabccccaaaaacccc,aaaacbccc,aaaaabccc,aaaaccccc,aaabacccc,aaabccccc,aaababccc,aaabcbccc`$. For distinct $`a,b`$ and $`c`$, fixing $`a`$ would give $`(q1)`$ choices for $`c`$ and $`(q2)`$ possibilities for $`b`$, where $`q`$ is the number of Potts states, yielding a total of $`(q1)(q2)`$ possible variations on b and c, given a.

Now notice that precisely two of these final states, $`aaaaacccc`$ and $`aaaaccccc`$, correspond to the reaction $`A+AA`$, with the $`A`$ particle transforming into a B particle at the next time step; see Fig. 3 for an illustration. All other states can be interpreted in terms of the diffusion of the $`A`$ and $`B`$ particles to neighbour sites, keeping in mind the interchange of A and B labels at succeeding time steps. Now if $`a=c`$ while $`ba,c`$, as in the initial configuration $`aaaabaaaaa`$, fixing $`a`$ would give $`(q1)`$ choices for $`b`$. However, if $`a=c`$, the spin flip at the “b” site yields, for the final configurations corresponding to the reaction $`A+AA`$ above, the reaction $`aaaabaaaaaaaaaaaaaaa`$, interpreted as $`A+A\varphi `$ in terms of our mapping.

Starting from an initial disordered configuration and averaging over all initial configurations, one would expect configurations of the type ($`aaaabcccc`$), with $`abc`$, and ($`aaaabaaaa`$) with $`ab`$, to appear in the relative ratio $`(q2)(q1):(q1)`$. In terms of the dynamics of the $`A`$ and $`B`$ particles, this translates into the reaction $`A+AA`$ occuring with probability $`(q2)/(q1)`$ vis a vis the reaction $`A+A\varphi `$ which occurs with probability $`1/(q1)`$. This justifies the rates which appear in our mapping to the system of decoupled reaction diffusion equations illustrated in Eq. 8.

A straightforward corollary to these results is the following: The persistence of a given spin, measured through configurations separated by two time steps, should decay as $`P(t)1/t^{\theta _s(q)}`$ i.e. with the serial exponent rather than the parallel one. This is a simple consequence of the fact that only one of the sublattices contributes to persistence in this case. This exact result is interesting in view of recent work on persistence with discrete sampling, in which an underlying continuous time (as opposed to discrete in the case studied here) process is sampled at regular time intervals and its persistence properties studied.

In conclusion, we have described a mapping between configurations of a one-dimensional $`q`$state ferromagnetic Potts model evolving under parallel dynamics and a system of two decoupled reaction-diffusion equations representing the motion and interaction of the boundary sites of stable domains. Using this mapping, we derive a general result connecting persistence exponents obtained with serial and parallel dynamics for this model. This mapping relied on identifying the crucial and somewhat unexpected role of what we have called unstable sites of type II. While the mapping here addressed the specific problem of the persistence exponent, it may be of more general use in the computation of other statistical properties of Potts models evolving under parallel dynamics.

We thank P. Shukla for useful dicussions in the context of Ref .


text:
# Bailout Embeddings, Targeting of KAM Orbits, and the Control of Hamiltonian Chaos

(version of )

## Abstract

We present a novel technique, which we term bailout embedding, that can be used to target orbits having particular properties out of all orbits in a flow or map. We explicitly construct a bailout embedding for Hamiltonian systems so as to target KAM orbits. We show how the bailout dynamics is able to lock onto extremely small KAM islands in an ergodic sea.

Control of chaos in nonlinear dynamical systems has been achieved by applying small perturbations that effectively change the dynamics of the system around the region — typically a periodic orbit — that one wishes to stabilize Ott et al. (1990). This method has been successful in dissipative systems, but its extensions to control and targeting in Hamiltonian systems Lai et al. (1993); Schroer and Ott (1997); Wu et al. (1998); Oloumi and Teychenné (1999); Macau (1998); Bolotin et al. (1999); Zhang et al. (2000); Zhang and Yao (2000); Xu et al. (2001) have met various difficulties not present in the dissipative case: the absence of attracting sets, for instance, makes it hard to stabilize anything. In addition, these methods require that one know in advance what one wants to do; in particular, the orbit to be stabilized may have to be known in advance to a fair accuracy.

In this Letter we present a novel technique allowing us to control Hamiltonian chaos, in such a way as to keep the original dynamics intact, but which shifts the stability of different kinds of orbits in the dynamics. We do so by embedding our Hamiltonian system within a larger space, meaning we augment the number of degrees of freedom, keeping an intact copy of the original system on one privileged slice; all of the control is achieved through use of the perpendicular directions to this intact copy. We call this method a bailout embedding for reasons that will promptly become clear. We apply the technique to an extremely hard and hitherto close to intractable problem in chaos control of Hamiltonian systems: selecting small KAM (Kolmogorov–Arnold–Moser) islands within chaotic seas in systems that are almost ergodic. We show below how our technique is able to find and render asymptotically stable minute islands of order within a map. While formerly this could be done by sophisticated algorithms on the basis of complex logic, our method distills the complexity of this calculation into a simple, forward-iterates dynamical system. Thus our method allows us to stabilize KAM islands without knowing their location in advance.

A striking application of Eq. (4) is to a case rather close to the original setting of Eq. (1): an incompressible fluid flow $`𝒖`$ translates in dynamical-systems terms to a divergence-free flow, of which Hamiltonian systems are the most important class. A classical problem in Hamiltonian dynamics is locating KAM tori. Hamiltonian systems live between two opposite extremes, of fully integrable systems and fully ergodic ones. Fully integrable systems are characterized by dynamics unfolding on invariant tori. The KAM theorem asserts that as a parameter taking the system away from integrability is increased, these tori break and give rise to chaotic regions in a precise sequence; for any particular value of this parameter in a neighborhood of the integrable case, there are surviving tori. The problem with finding them is that, the dynamics being volume-preserving, merely evolving trajectories either forwards or backwards does not give us convergence onto tori Gutzwiller (1990), and since for large values of the nonlinearity, they cover a very small measure of the phase space, locating them becomes an extremely difficult problem. Moreover, there is the further problem that even if we start on an island, we should be able to recognize it as such. In fact, several sophisticated analytical and semianalytical techniques had to be developed to assist in the search and characterization of KAM tori in highly nonlinear Hamiltonian systems (see for instance Percival (1979)). We shall now show that the bailout embedding solves these problems by transforming the KAM trajectories into global attractors of the embedded system; finding them is now independent of the choice of initial conditions.

Our technique is based upon the dynamics of a small neutrally buoyant sphere — a passive scalar — in an incompressible fluid flow $`𝒖`$ Babiano et al. (2000). Under assumptions allowing us to retain only the Bernoulli, Stokes drag, and Taylor added mass contributions to the force exerted by the fluid on the sphere, the equation of motion for the sphere at the position $`𝒙`$ is

$$\frac{d}{dt}\left(\dot{𝒙}𝒖(𝒙)\right)=(\lambda +\mathbf{}𝒖)\left(\dot{𝒙}𝒖(𝒙)\right)$$
(1)
so that the difference between the particle velocity and the velocity of the surrounding fluid is exponentially damped with negative damping constant $`(\lambda +\mathbf{}𝒖)`$. In the case in which the flow gradients reach the magnitude of the viscous drag coefficient, there is the possibility that around hyperbolic stagnation points the Jacobian matrix $`\mathbf{}𝒖`$ may acquire a positive eigenvalue in excess of the drag coefficient. In these instances, the trajectories of these passive scalars, instead of converging exponentially onto $`\dot{𝒙}=𝒖`$, detach from such trajectories. The result is that the passive scalars explore practically all of the flow, but tend to stay away from regions of high shear.

Equation 1 is but one instance of a more general structure. Let us consider a differential equation of the form $`\dot{x}=f(x)`$. If we take its time derivative, we obtain a different differential equation, $`\ddot{x}=f^{}(x)\dot{x}`$. It is different because, being second order, it exists in a larger space and has many solutions that are not solutions of the smaller one. Still, the original equation is contained within the larger system, in the sense that every solution of $`\dot{x}=f(x)`$ is a solution of $`\ddot{x}=f^{}(x)\dot{x}`$. We may say that $`\dot{x}=f(x)`$ is embedded within $`\ddot{x}=f^{}(x)\dot{x}`$. There are infinitely many ways to embed; for example, $`\ddot{x}=f^{}(x)f(x)`$ is also an embedding, but it is clearly inequivalent to the first.

Of course, embedding a system changes notions of stability, because stability refers to perturbations, and in a larger system there are all of the old perturbations plus a batch of new ones. So, even though all of the solutions of the original system are preserved, by adding new directions away from the old solutions we may transform formerly stable solutions into unstable ones in the larger setting. See, for example, studies of the manifold bubbling transition Venkataramani et al. (1996). The trivial way to embed a system is through a cross product; for instance, $`\dot{x}=f(x),x𝕄`$ is embedded within $`𝕄\times `$ as

$`\dot{x}`$ $`=`$ $`f(x)+g(x,y),`$
$`\dot{y}`$ $`=`$ $`\alpha y,`$ (2)
where $`g(x,y)`$ is arbitrary except for requiring that $`g(x,0)0`$, which guarantees that for $`y=0`$ we have the original system. If $`\alpha <0`$ then $`y`$ always dies out, and so we always recover the original object; in this case, we can call the embedding itself stable, in the sense that any motion away from the embedded object takes us back. Both of the “derivative” embeddings in the previous paragraph were unstable; stable versions can be constructed rather simply, for instance

$$\frac{d}{dt}\left(\dot{x}f(x)\right)=k\left(\dot{x}f(x)\right),$$
(3)
of which the previous examples were the $`k=0`$ limit. This embedding ensures that the distance between the actual trajectory and the embedding diminishes exponentially with time for any initial condition. Equation (3) begins to resemble Eq. (1).

We define a bailout embedding as one of the form

$$\frac{d}{dt}\left(\dot{x}f(x)\right)=k(x)\left(\dot{x}f(x)\right),$$
(4)
where $`k(x)<0`$ on a set of orbits that are unwanted, and $`k(x)>0`$ otherwise. Thus the natural behavior of a bailout embedding is that the trajectories in the full system tend to detach or bail out from the embedded subsystem into the larger space, where they bounce around. If, after bouncing around for a while, these orbits reach a stable region of the embedding, $`k(x)>0`$, they will once again collapse onto the original dynamical system. In this way we can create a larger version of the dynamics in which specific sets of orbits are removed from the asymptotic set, while preserving the dynamics of another set of orbits — the wanted or targeted one — as attractors of the enlarged dynamical system. For the special choice of $`k(x)=(\lambda +f)`$, as in the fluid dynamics of a passive scalar described by Eq. (1), these dynamics were shown in Babiano et al. (2000) to detach from saddle points and other unstable regions in conservative dynamics.

It is not hard to extend flow bailout, Eq. (4), to maps in the obvious fashion. Given a map $`x_{n+1}=T(x_n)`$ the bailout embedding is given by

$$x_{n+2}T(x_{n+1})=K(x_n)(x_{n+1}T(x_n))$$
(5)
provided that $`|K(x)|>1`$ over the unwanted set. (In the map system, almost any expression written for the ordinary differential equation translates to something close to an exponential; in particular, stability eigenvalues have to be negative in the ordinary differential equation case to represent stability, while they have to be smaller than 1 in absolute value in the map case). The particular choice of the gradient as the bailout function $`k(x)=(\lambda +f)`$ in a flow translates in the map setting to $`K(x)=e^\lambda T`$. A classical testbed of Hamiltonian systems is the standard map, an area-preserving map introduced by Taylor and Chirikov. The standard map is given by

$`x_{n+1}`$ $`=`$ $`x_n+{\displaystyle \frac{k}{2\pi }}\mathrm{sin}(2\pi y_n)`$
$`y_{n+1}`$ $`=`$ $`y_n+x_{n+1}`$ (6)
where $`k`$ is the parameter controlling integrability.

In order to embed the standard map, we only need to replace $`T`$ and $`K(x)`$ in Eq. (5) with the appropriate expressions. $`T`$ stems directly from Eq. (6) and, in accordance with the previous definitions, $`K(x)`$ becomes:

$$K(x)=\left(\begin{array}{cc}1& k\mathrm{cos}2\pi y_n\\ 1& k\mathrm{cos}2\pi y_n+1\end{array}\right)$$
(7)
Notice that due to the area preserving property of the standard map, the two eigenvalues of the derivative matrix must multiply to one. If they are complex, this means that both have an absolute value of one, while if they are real, generically one of them will be larger than one and the other smaller. We can then separate the phase space into elliptic and hyperbolic regions corresponding to each of these two cases. If a trajectory of the original map lies entirely on the elliptic regions, the overall factor $`\mathrm{exp}(\lambda )`$ will damp any small perturbation away from it in the embedded system. But for chaotic trajectories, which inevitably must visit some hyperbolic regions, there exists a value of $`\lambda `$ such that perturbations away from a standard map-trajectory are amplified instead of dying out in the embedding. As a consequence, trajectories are effectively expelled from the chaotic regions to finally settle in the safely elliptic KAM islands. This process can be observed clearly in Figure 1. As the value of $`\lambda `$ is decreased, the number of trajectories starting from random initial conditions which eventually settle into the KAM tori increases; see Figure 2.

This process is sensitive enough to work even at hitherto intractably high nonlinearities. For $`k=7`$, the standard map is almost ergodic: it covers almost the entire torus with a single chaotic orbit. Only a tiny island of irreductible order resists this invasion. It is located around $`(0,0.6774)`$ and covers an area approximately $`2\times 10^5`$ of the torus. Thus, from random initial conditions one would expect to see this island only once every 50000 attempts. Figure 3 shows how easily our bailout method finds this island from simple forward iterates.

A Hamiltonian system does not usually just satisfy volume conservation, but also will conserve the Hamiltonian itself. Given an ordinary differential equation $`\dot{x}=f(x)`$ with a conserved quantity $`E0`$, then $`fE=0`$. However, building a bailout embedding by the procedure above does not lead to dynamics that satisfy $`\dot{E}0`$, because the bailout embedding should be $`2n2`$ dimensional. This is clearly undesirable in the case of Hamiltonian systems, so we show now how to derive a bailout embedding that will obey a conservation law. The bailout equation can be written

$$\ddot{x}=(f\lambda )(\dot{x}f)+f\dot{x}$$

We need to correct this acceleration so that it stays on the second tangent space of the $`E0`$ surface. Let us call the raw bailout acceleration $`u`$. The second derivative $`\ddot{x}`$ has to satisfy $`\ddot{x}E+\dot{x}E\dot{x}=0`$, so we can modify $`u`$ to

$$\ddot{x}=u\frac{uE}{|E|^2}E+\frac{E}{|E|}\dot{x}E\dot{x}$$

This equation, given that we start on $`\dot{x}E=0`$, will then preserve this property.

We have presented a novel method in the control and targeting of chaos in nonlinear dynamical systems, the bailout embedding. While potentially useful in any dynamical-system setting, this method is especially suited to Hamiltonian systems. Unlike other chaos-control strategies, this method does not obliterate the original dynamics of the system, but rather preserves it in a privileged slice of phase space embedded in a higher-dimensional space, and merely shifts around the stability of its orbits. A suitable choice of a bailout function allows this strategy to target a complex set of orbits. We have shown, in particular, the targeting of KAM orbits, a case well-known from classical studies to be especially hard.

We should like to thank Leo Kadanoff and Marcelo Viana for useful discussions. JHEC acknowledges the financial support of the Spanish CSIC, Plan Nacional del Espacio contract ESP98-1347. MOM acknowledges the support of the Meyer Foundation. OP acknowledges the Spanish Ministerio de Ciencia y Tecnologia, Proyecto CONOCE, contract BFM2000-1108.