text:
# On the Global Dynamics of the Anisotropic Manev Problem

## 1 Introduction

The anisotropic Manev problem describes the motion of two point masses in an anisotropic configuration plane under the influence of a Newtonian force-law with a relativistic correction term. The isotropic case is the classical Manev problem; its origins lie in the work of Newton, who introduced it in Principia aiming to understand the apsidal motion of the moon (see ). Manev found in the 1930s that a proper choice of the constants that show up in the correction term allow the theoretical explanation of the perihelion advance of Mercury and of the other inner planets.

The first author suggested the study of the anisotropic Manev problem in 1995, hoping to find connections between classical, quantum, and relativistic mechanics. It was indeed proved in  that the rich collision-orbit manifold of the system exhibits classical, quantum, and relativistic properties. This encouraged further studies, as for example  and . In , using a suitable generalization of the Poincaré-Melnikov method (see  for the classical approach or  for a parallel, at least in part, complementary approach), we proved that chaos occurs on the zero-energy manifold, thus showing the complexity of the dynamics. Using perturbations techniques and the Poincaré continuation method, the second author investigated in  the classes of periodic solutions that arise from symmetries in the case of small values of the anisotropy parameter.

In this paper we gain a better understanding of the complicated global dynamics encountered in this problem. We first prove that negative-energy solutions are bounded and find the heteroclinic orbits that connect the equilibria of the collision manifold, which we obtain through McGehee-type transformations (see ). Physically they correspond to ejection-collision orbits. Then we employ perturbation techniques to detect possible global chaotic behaviour. As remarked in , the perturbation analysis of  cannot be used to study ejection-collision solutions. However, we surpass this difficulty with the help of McGehee-type coordinates, which allow us to view the anisotropic Manev problem as a perturbation of the classical Manev case.

Using an approach inspired by , which works in some degenerate cases—as for example those of unstable nonhyperbolic points or critical points located at infinity (see ), we develop a suitable extension of the Poincaré-Melnikov method, which we use to prove the existence of transversal homoclinic orbits to a periodic one. It is interesting to note that our result extends the one obtained in  to a non-Hamiltonian system that has negatively and positively asymptotic sets to a nonhyperbolic periodic orbit. In the present context the asymptotic sets are the stable and the unstable manifolds.

Then we return to the original coordinates and apply a variational principle for detecting periodic orbits. Using the rotation index, we divide the set of periodic paths into homotopy classes, which are Sobolev spaces. Then we use the lower semicontinuity version of Hilbert’s direct method (due to Tonelli, see ) to find a minimizer of the action in each class. According to the least action principle, the minimizer is a solution of the anisotropic Manev problem. We prove that the minimizer exists, belongs to the homotopy class, and is a solution in the classical sense. This generalizes a result obtained by the second author, , where it was shown that such orbits exist for small values of $`\mu >1`$. In the end we put into the evidence some new properties of symmetric periodic orbits.

The idea of using variational principles to obtain periodic orbits for $`n`$-body-type particle systems first appeared in  and has been recently used in connection with symmetry conditions to obtain new periodic orbits in the classical $`n`$-body problem (see ). But unlike the Newtonian case, the Manev force is “strong” (as defined in ), so the variational method is easier to apply in our situation than in the Newtonian one. This is because in the Manev case we do not have to deal with the difficulty of avoiding collision orbits, which have infinite action and therefore cannot be minimizers.

Our paper is organized as follows. In Section 2 we write the equations of motion and transform them to an equivalent system using a “blow-up” technique devised by McGehee, which allow us to introduce the concept of a collision manifold. In Section 3 we present two global results: the boundedness of the solutions for negative energy and the existence of certain symmetric ejection-collision orbits. In Section 4 we describe the anisotropic Manev problem as a perturbation of the Manev case. In Section 5 we develop a suitable generalization of the Poincaré-Melnikov method and in Section 6 we apply it to find infinitely many transverse homoclinic orbits that show that the dynamics of the problem is extremely complex, possibly chaotic. Finally, in Section 7 we use a variational principle to prove the existence of infinitely many classes of symmetric periodic orbits.

## 2 Equations of Motion

The (planar) anisotropic Manev problem is described by the Hamiltonian

$$H=\frac{1}{2}𝐩^2\frac{1}{\sqrt{x^2+\mu y^2}}\frac{b}{x^2+\mu y^2},$$
(1)
where $`\mu >1`$ is a constant, $`𝐪=(x,y)`$ is the position of one body with respect to the other considered fixed at the origin of the coordinate system, and $`𝐩=(p_x,p_y)`$ is the momentum of the moving particle. The constant $`\mu `$ measures the strength of the anisotropy and we can very well take $`\mu <1`$; but to remain consistent with the choice made in previous papers, we will consider $`\mu >1`$. For $`\mu =1`$ we recover the classical Manev problem. The equations of motion are

$$\{\begin{array}{c}\dot{𝐪}=𝐩\hfill \\ \dot{𝐩}=\frac{H}{𝐪}.\hfill \end{array}$$
(2)
The Hamiltonian provides the first integral

$$H(𝐩(t),𝐪(t))=h,$$
(3)
where $`h`$ is a real constant. Unlike in the classical Manev case, the angular momentum $`K(t)=𝐩(t)\times 𝐪(t)`$ does not yield a first integral. This is because the anisotropy of the plane destroys the rotational invariance.

Since our first goal is to study collision and near collision solutions, it is helpful to transform system (2) using a method developed by McGehee . The idea is to “blow-up” the collision singularity, replace it with a so-called collision manifold and extend the phase space to it. The collision manifold is fictitious in the sense that it has no physical meaning. However, studying the flow on it provides useful information about near-collision orbits. Consider the coordinate transformations

$$\{\begin{array}{c}r=|𝐪|\hfill \\ \theta =\mathrm{arctan}(y/x)\hfill \\ v=\dot{r}r=(xp_x+yp_y)\hfill \\ u=r^2\dot{\theta }=(xp_yyp_x),\hfill \end{array}$$
(4)
and the rescaling of time

$$d\tau =r^2dt.$$
(5)
Composing these transformations, which are analytic diffeomorphisms in their respective domains, system (2) becomes

$$\{\begin{array}{c}r^{}=rv\hfill \\ v^{}=2r^2h+r\mathrm{\Delta }^{1/2}\hfill \\ \theta ^{}=u\hfill \\ u^{}=(1/2)(\mu 1)(r\mathrm{\Delta }^{3/2}+2b\mathrm{\Delta }^2)\mathrm{sin}2\theta \hfill \end{array}$$
(6)
and the energy relation (3) takes the form

$$u^2+v^22r\mathrm{\Delta }^{1/2}2b\mathrm{\Delta }^1=2r^2h,$$
(7)
where $`\mathrm{\Delta }=\mu \mathrm{cos}^2\theta +\mathrm{sin}^2\theta `$ and the new variables $`(r,v,\theta ,u)(0,\mathrm{})\times \times S^1\times `$ depend on the fictitious time $`\tau `$. The prime denotes differentiation with respect to $`\tau `$.

The set

$$C=\{(r,v,\theta ,u)|r=0\text{ and the energy relation (}\text{7}\text{) holds}\}$$
(8)
is the collision manifold, which replaces the set of singularities $`\{(𝐪,𝐩)|𝐪=\mathrm{𝟎}\}`$. This 2-dimensional manifold, embedded in $`^3\times S^1`$, is homeomorphic to a torus and it is given by the equations

$$r=0\text{and}u^2+v^2=2b\mathrm{\Delta }^1.$$
(9)
The flow on the collision manifold was studied in detail in . Here we will briefly recall its main features. Let’s consider the restriction of system (6) to $`C`$. The solutions of the restriction lie on the level curves $`v=\text{constant}`$ of the torus $`C`$. There are eight equilibrium points. In the variables $`(r,v,\theta ,u)`$ the first four equilibria are $`A_0^\pm =(0,\pm \sqrt{2b/\mu },0,0)`$ and $`A_\pi ^\pm =(0,\pm \sqrt{2b/\mu },\pi ,0)`$. The corresponding eigenvalues are real and take the values $`\pm \sqrt{2b/\mu },0,\pm \sqrt{2b(1\mu )/\mu }`$. The other four equilibria are $`A_{\pm \pi /2}^\pm =(0,\pm \sqrt{2b},\pm \pi /2,0)`$ and the corresponding eigenvalues are $`\pm \sqrt{2b},0,\pm \sqrt{2b(1\mu )}`$, where the last two eigenvalues are purely imaginary since $`\mu >1`$. Moreover there are eight heteroclinic orbits which lie in the level sets $`v=\pm \sqrt{2b/\mu }`$. All the other solutions are periodic (see Fig. 1).

## 3 Heteroclinic Orbits and Bounded Solutions

The flow near the collision manifold was studied in , in which most of the results are essentially local. In this section we will prove two global results that extend the understanding of the problem under discussion. The first one concerns the boundedness of solutions on negative energy levels.

###### Theorem 1

For any negative value of the energy constant, $`h<0`$, there exists a positive real number $`M`$ such that any given solution $`(r(\tau ),v(\tau ),\theta (\tau ),u(\tau ))`$ of system (6) satisfies the relation $`r<M`$.

Proof: Let us assume that there is no $`M`$ with the above property. Then at least one unbounded solution exists. Since by the energy relation (7), $`u^2+v^2=2r^2h+2r\mathrm{\Delta }^{1/2}+2b\mathrm{\Delta }^1`$, and since $`h<0`$, there is some $`\overline{r}=r(\overline{\tau })`$ such that $`u^2+v^2`$ is negative—a contradiction. This completes the proof.

The next result deals with the existence of heteroclinic orbits connecting the equilibria but lying outside the collision manifold. But before stating and proving it, let us recall some facts that summarize the behavior of the flow near the collision manifold. Denote by $`P_\eta `$ the periodic orbit on $`C`$ having $`v=\eta `$. The following property was proved in .

###### Proposition 1

On the collision manifold $`C`$ the equilibria $`A_0^\pm `$ and $`A_\pi ^\pm `$ are saddles whereas the equilibria $`A_{\pm \pi /2}^\pm `$ are centers. Outside the collision manifold the equilibria $`A_0^\pm `$, $`A_{\pm \pi /2}^+`$, and $`A_\pi ^+`$ have a $`1`$-dimensional unstable analytic manifold, whereas the equilibria $`A_0^{}`$, $`A_{\pm \pi /2}`$, and $`A_\pi ^{}`$ have a $`1`$-dimensional stable analytic manifold. Each periodic orbit $`P_\eta `$ on $`C`$ with $`v=\eta >0`$($`v=\eta <0`$) has a $`2`$-dimensional local unstable analytic manifold, while the periodic orbit $`v=0`$ has both a $`2`$-dimensional local unstable and a $`2`$-dimensional local stable manifold (see Fig. 2)

The above properties are local, the following one, however, is global. We will now show that the equilibria with positive $`v`$ coordinate have a 1-dimensional global unstable manifold while the equilibria with a negative $`v`$ have a 1-dimensional stable manifold. Moreover, the equilibria are connected by heteroclinic orbits starting from an equilibrium with positive $`v`$ and ending in the symmetric one with respect to the $`(\theta ,u)`$ plane.

###### Theorem 2

There are four heteroclinic orbits outside the collision manifold $`C`$: $`\gamma _{\pi /2}`$ ,$`\gamma _0`$ ,$`\gamma _{\pi /2}`$ ,$`\gamma _\pi `$ connecting respectively $`A_{\pi /2}^+`$ with $`A_{\pi /2}^{}`$, $`A_0^+`$ with $`A_0^{}`$, $`A_{\pi /2}^+`$ with $`A_{\pi /2}^{}`$, and $`A_\pi ^+`$ with $`A_\pi ^{}`$ (see Fig. 2)

Proof: First we show that $`u=0`$ and $`\theta =0,\pi ,\pm \pi /2`$ describe four invariant sets. Consider $`\theta (0)=\theta _0=0`$ and $`u(0)=u_0=0`$, as initial conditions. Then $`\theta 0`$, $`u0`$ satisfies system (6), hence (by the uniqueness property for solutions) $`u=0`$, $`\theta =0`$ define an invariant set. The same reasoning can be applied if $`\theta =\pm \pi /2`$ or $`\pi `$.

Now let’s study the energy relation (7) when $`u=0`$ and $`\theta =0,\pi `$. After simple computations we get

$$v^2+\left(\sqrt{2|h|}r\frac{1}{\sqrt{2\mu |h|}}\right)^2=\frac{2b}{\mu }+\frac{1}{2\mu |h|}.$$
(10)
The above equation describes an ellipse whose intersections with $`r=0`$ give $`v=\pm \sqrt{2b/\mu }`$, which are exactly the equilibrium points $`A_0^\pm `$ and $`A_\pi ^\pm `$. Moreover the maximum value of $`|v|`$ is

$$v_{max}=\sqrt{\frac{1}{\mu }\left(2b+\frac{1}{2|h|}\right)}$$
(11)
and the maximum value of $`r`$, attained when $`v=0`$, is

$$r_{max}=\frac{1}{(2\sqrt{\mu }|h|)}+\frac{\sqrt{\frac{1}{\mu }(14hb)}}{2|h|}.$$
(12)
Consequently for $`u=0`$, $`\theta =0`$ ($`\theta =\pi `$) there exist heteroclinic orbits $`\gamma _0`$, ($`\gamma _\pi `$) ejecting from $`A_0^+`$ ($`A_0^{}`$) and tending to $`A_0^{}`$ ($`A_\pi ^{}`$) (see Fig. 2).

Similarly when $`u=0`$ and $`\theta =\pm \pi /2`$ the energy relation can be reduced to the form

$$v^2+\left(\sqrt{2|h|}r\frac{1}{\sqrt{2|h|}}\right)^2=2b+\frac{1}{2|h|},$$
(13)
which describes an ellipse. The intersections with $`r=0`$ are $`v=\pm \sqrt{2b}`$ and represent the equilibria $`A_{\pm \pi /2}^\pm `$. In this case

$$v_{max}=\sqrt{2b+1/2|h|}$$
(14)
and

$$r_{max}=\frac{1}{2|h|}+\frac{\sqrt{1+4b|h|}}{2|h|}.$$
(15)
Thus we found heteroclinic orbits $`\gamma _{\pm \pi /2}`$ ejecting from $`A_{\pm \pi /2}^+`$ and tending to $`A_{\pm \pi /2}^{}`$ (see Fig. 2). This completes the proof.

## 4 A Perturbative Approach

We will now write the anisotropic Manev problem as a perturbation of the classical Manev case. Consider weak anisotropies, i.e., choose the parameter $`\mu `$ close to $`1`$. Introducing the notation $`\mu 1=ϵ>0`$ with $`ϵ1`$, we can expand the equation of motion in powers of $`ϵ`$ to obtain

$$\{\begin{array}{c}r^{}=rv\hfill \\ v^{}=2r^2h+rϵ(r/2\mathrm{cos}^2\theta )\hfill \\ \theta ^{}=u\hfill \\ u^{}=ϵ(r/2+2b)\mathrm{sin}2\theta .\hfill \end{array}$$
(16)
The energy relation becomes

$$u^2+v^22r2b+ϵ(r+2b)\mathrm{cos}^2\theta =2r^2h.$$
(17)
For $`ϵ=0`$, system (16) and equation (17) yield the Manev problem. The collision manifold is the set of solutions given by

$$r=0,u^2+v^2=2b.$$
(18)
Notice that, from the geometric point of view, the collision manifold is a cylinder in the three-dimensional space of coordinates $`(u,\theta ,v)`$ and, since $`\theta [0,2\pi ]`$, it follows that this cylinder can be identified with a torus. The flow on the collision manifold is formed almost exclusively by non-hyperbolic periodic orbits, except for the upper and lower circles of the torus given by $`r=0,u=0,v=\pm \sqrt{2b}`$, which consist of equilibrium points. There is only a single orbit ejecting from each fixed point of the upper circle $`v=\sqrt{2b}`$ and a single orbit tending to the lower circle $`v=\sqrt{2b}`$ (see ). Moreover it can be easily proved (see ) that for every periodic orbit $`p_v`$ on the collision manifold with $`0<v<\sqrt{2b}`$ there exist a manifold of orbits, lying on a cylinder, which eject from $`p_v`$. Similarly it can be shown that for every orbit $`p_v`$, with $`\sqrt{2b}<v<0`$, there exists a manifold of orbits, lying on a cylinder, which tend to $`p_v`$.

If $`v=0`$ both types of manifolds exist, so $`p_0`$ has a homoclinic manifold. Indeed, the equations that describe the manifold can be found explicitly: they have $`u=\pm \sqrt{2b}`$. With the energy relation we get

$$v=\pm \sqrt{2r^2h+2r},$$
(19)
and using the equation of motion we obtain

$$r^{}=\pm r\sqrt{2r^2h+2r}.$$
(20)
By integrating equation (20) it is easy to find that

$$R(\tau \tau _0)=\frac{2}{2|h|+(\tau \tau _0)^2},R^{}=\frac{4(\tau \tau _0)}{(2|h|+(\tau \tau _0)^2)^2}$$
(21)
and

$$V(\tau \tau _0)=\frac{R^{}}{R}=\frac{2(\tau \tau _0)}{2|h|+(\tau \tau _0)^2}.$$
(22)
Furthermore

$$U(\tau \tau _0)=\pm \sqrt{2b}=\omega \text{and}\vartheta (\tau \tau _0,\theta _0)=\mathrm{\Theta }(\tau \tau _0)\theta _0,$$
(23)
where $`\mathrm{\Theta }(\tau \tau _0)=\omega (\tau \tau _0)`$. As $`\tau _0`$ and $`\theta _0`$ vary, equations (21-23) describe the entire $`2`$-dimensional homoclinic manifold. An orbit lying on the homoclinic manifold is represented in Fig. 3. Such an orbit is obtained by choosing $`\theta _0=0`$; it ejects from the equator of the collision manifold, spiraling around it and moving upwards, then changes directions, goes downwards and upwards again, spiraling towards the periodic orbit $`p_0`$.

The homoclinic manifold plays an important role in following section and is necessary for developing the generalization of the Melnikov technique.

## 5 A Generalized Melnikov Method

Let $`\chi =(R(\tau ),V(\tau ),\mathrm{\Theta }(\tau ),U(\tau ))`$ be the homoclinic orbit selected when we choose $`\tau _0=0`$ and $`\theta _0=0`$. Consider solutions of the form

$$\{\begin{array}{c}r(\tau ,\tau _0)=R(\tau \tau _0)+\stackrel{~}{r}(\tau ,\tau _0)\hfill \\ v(\tau ,\tau _0)=V(\tau \tau _0)+\stackrel{~}{v}(\tau ,\tau _0)\hfill \\ \theta (\tau ,\tau _0,\theta _0)=\mathrm{\Theta }(\tau \tau _0)\theta _0+\stackrel{~}{\theta }(\tau ,\tau _0)\hfill \\ u(\tau ,\tau _0)=U(\tau \tau _0)+\stackrel{~}{u}(\tau ,\tau _0).\hfill \end{array}$$
(24)
Let $`\stackrel{~}{𝐳}=(\stackrel{~}{r},\stackrel{~}{v},\stackrel{~}{\theta },\stackrel{~}{u})`$, then the variational equation is

$$\stackrel{~}{𝐳}^{}=A(\tau )\stackrel{~}{𝐳}+\stackrel{~}{𝐛}(\stackrel{~}{𝐳},\chi ,\tau ,\tau _0,\theta _0,ϵ),$$
(25)
where

$$A(\tau )=\left(\begin{array}{cccc}V& R& 0& 0\\ 1+4Rh& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 0& 0\end{array}\right)$$
(26)
and

$$\stackrel{~}{𝐛}(\stackrel{~}{𝐳},\chi ,\tau ,\tau _0,\theta _0,ϵ)=\left(\begin{array}{c}b_1\\ b_2\\ b_3\\ b_4\end{array}\right)=\left(\begin{array}{c}\stackrel{~}{r}(\tau \tau _0)\stackrel{~}{v}(\tau \tau _0)\\ ϵ\left(\frac{(R+\stackrel{~}{r})}{2}\mathrm{cos}^2(\mathrm{\Theta }\theta _0+\stackrel{~}{\theta })\right)\\ 0\\ ϵ\left(\frac{(R+\stackrel{~}{r})}{2}+2b\right)\mathrm{sin}2(\mathrm{\Theta }\theta _0+\stackrel{~}{\theta })\end{array}\right).$$
(27)
The general solution of the variational equation (25) is

$$\stackrel{~}{𝐳}=\mathrm{\Phi }(t)_{t_0}^t\mathrm{\Phi }^1(s)\stackrel{~}{𝐛}𝑑s,$$
(28)
(see ), where $`\mathrm{\Phi }`$ is the fundamental matrix. If we let $`c=\mathrm{\Phi }^1\stackrel{~}{𝐛}`$, the previous equation becomes

$$\stackrel{~}{z}_i(t)=\mathrm{\Phi }_{ij}_{t_0}^tc_j(s)𝑑s,$$
(29)
where $`c_j=detD_j(t)/(det\mathrm{\Phi })(t)`$ and $`D_j`$ is the matrix obtained replacing the $`j`$-th column of $`\mathrm{\Phi }`$ with $`\stackrel{~}{𝐛}`$. Furthermore the following formula for the trace holds:

$$\text{det}\mathrm{\Phi }(\tau )=Ce^{_{\tau _0}^\tau TrA(s)𝑑s}.$$
(30)
One solution of the homogeneous part of the variational equation is given by

$$\chi ^{}(\tau \tau _0,\theta _0)=(R^{}(\tau \tau _0),V^{}(\tau \tau _0),\mathrm{\Theta }^{}(\tau \tau _0),U^{}(\tau \tau _0)),$$
(31)
where

$$\begin{array}{c}R^{}=\frac{4(\tau \tau _0)}{(2|h|+(\tau \tau _0)^2)^2}\hfill \\ V^{}=\frac{2}{2|h|+(\tau \tau _0)^2}+\frac{4(\tau \tau _0)^2}{(2|h|+(\tau \tau _0)^2)^2}\hfill \\ \mathrm{\Theta }^{}=\pm \sqrt{2b}.\hfill \\ U^{}=0.\hfill \end{array}$$
(32)
It is easy to check that other two independent solutions are $`(0,0,1,0)`$ and $`(0,0,0,1)`$. Knowing three independent solutions of a linear system, it is possible to find a fourth independent solution $`\psi `$. This is achieved through the following lemma, which will be used to estimate how fast $`\psi `$ diverges.

###### Lemma 1

Let $`\stackrel{~}{𝐳}^{}=A\stackrel{~}{𝐳}`$ be the homogeneous part of (25). Given the three independent solutions above, a fourth is defined by

$$\{\begin{array}{c}\stackrel{~}{z}_1^{}=(1+4Rh)/V^{}\stackrel{~}{z}_4\hfill \\ \stackrel{~}{z}_2^{}=\frac{\sqrt{2b}}{V^{}}(1+4Rh)\stackrel{~}{z}_4\hfill \\ \stackrel{~}{z}_3^{}=0\hfill \\ \stackrel{~}{z}_4^{}=(V\frac{R^{}(1+4Rh)}{V^{}})\stackrel{~}{z}_4\hfill \end{array}\{\begin{array}{c}\psi _1=R^{}\stackrel{~}{z}_1+\stackrel{~}{z}_4\hfill \\ \psi _2=V^{}\stackrel{~}{z}_1\hfill \\ \psi _3=\pm \sqrt{2b}\stackrel{~}{z}_1+\stackrel{~}{z}_2\hfill \\ \psi _4=\stackrel{~}{z}_3.\hfill \end{array}$$
(33)
Proof: To find the fourth independent solution, we can use the “reduction to a smaller system,” (see ), whose direct application completes the proof.

In particular it is useful to remark that we can always choose $`\psi _4=0`$, since $`(\psi _1,\psi _2,\psi _3,0)`$ is a solution that is independent from the other ones. To obtain necessary and sufficient conditions such that the negatively and positively asymptotic sets intersect transversely, we first obtain conditions for the existence of solutions bounded on $``$ for the non-homogeneous linear variational equation around $`\chi `$.

For this, let $`()=\{\stackrel{~}{𝐛}:\times \times S^1\times \text{bounded, continuous}\}`$ with $`\stackrel{~}{𝐛}=\text{sup}_\tau \stackrel{~}{𝐛}(\tau )`$ for $`\stackrel{~}{𝐛}()`$. Then we have the following version of the Fredholm alternative for solutions bounded on $``$ (see  for a similar approach).

###### Lemma 2

Let $`\stackrel{~}{𝐳}\times \times S^1\times `$ and assume that $`\stackrel{~}{𝐳}0`$ in the expression of the function $`\stackrel{~}{𝐛}`$. Then the variational equation

$$\stackrel{~}{𝐳}^{}=A(\tau )\stackrel{~}{𝐳}+\stackrel{~}{𝐛}(\stackrel{~}{𝐳},\chi ,\tau ,\tau _0,\theta _0,ϵ)$$
(34)
has a bounded solution if and only if

$$_{\mathrm{}}^+\mathrm{}e^{_{\tau _0}^\tau TrA(s)𝑑s}R^{}(\tau \tau _0)b_2(\chi ,\tau ,\tau _0,\theta _0,ϵ)𝑑\tau =0.$$
(35)
The solution is unique and continuous and has the form $`\stackrel{~}{𝐳}=𝐋(\stackrel{~}{𝐛})+𝐰`$, where $`𝐋`$ is a bounded linear operator, $`𝐰=(0,0,\stackrel{~}{\theta }(\tau _0),\stackrel{~}{u}(\mathrm{}))`$, when $`\stackrel{~}{r}(\tau _0)R^{}(\tau _0)+\stackrel{~}{v}(\tau _0)V^{}(\tau _0))=0`$, and $`b_4`$ satisfies the relation below,

$$_{\mathrm{}}^+\mathrm{}b_4(\chi ,\tau ,\tau _0,\theta _0,ϵ)𝑑\tau =0.$$
(36)
Proof: Using Lemma 1 it is easy to determine the behavior of $`\psi `$ as $`\tau \pm \mathrm{}`$, precisely,

$$\tau \pm \mathrm{}\{\begin{array}{c}\stackrel{~}{z}_1\tau ^4\hfill \\ \stackrel{~}{z}_2\tau ^4\hfill \\ \stackrel{~}{z}_3\text{const.}\hfill \\ \stackrel{~}{z}_4\tau ,\hfill \end{array}\tau \pm \mathrm{}\{\begin{array}{c}\psi _1\tau \hfill \\ \psi _2\tau ^2\hfill \\ \psi _3\tau ^4\hfill \\ \psi _4\text{const.}\hfill \end{array}$$
(37)
Using (29) and (30), the general solution of the complete (non-homogeneous) equation (25) can be written in integral form as

$$\begin{array}{cc}\hfill \stackrel{~}{r}& =R^{}\left(A_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }(\psi _1b_2\psi _2b_1)𝑑s\right)\hfill \\ & +\psi _1\left(B+_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }(R^{}b_2V^{}b_1)𝑑s\right)\hfill \\ \hfill \stackrel{~}{v}& =V^{}\left(A_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }(\psi _1b_2\psi _2b_1)𝑑s\right)\hfill \\ & +\psi _2\left(B+_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }(R^{}b_2V^{}b_1)𝑑s\right)\hfill \\ \hfill \stackrel{~}{\theta }& =\pm \sqrt{2b}\left(A_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }(\psi _1b_2\psi _2b_1)𝑑s\right)\hfill \\ & +\psi _3\left(B+_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }(R^{}b_2V^{}b_1)𝑑s\right)\hfill \\ & +C_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }\left[(V^{}\psi _3\pm \sqrt{2b}\psi _2)b_1+(R^{}\psi _3\pm \sqrt{2b}\psi _1)b_2\right]\hfill \\ \hfill \stackrel{~}{u}& =D+_{\tau _0}^\tau b_4𝑑s,\hfill \end{array}$$
(38)
where, for notational convenience, we failed to mention the dependence on $`\stackrel{~}{𝐳}`$, $`\chi `$, $`\tau _0`$, etc.

Consider now the linearization of the problem (38) around the solution $`\stackrel{~}{𝐳}(\tau )0`$; in particular this amounts to deleting the high-order terms in the expression of $`\stackrel{~}{𝐛}`$ (i.e. $`b_1=0`$, etc.). Taking also into account the different behavior of the different solutions given in Lemma 1, it is easy to see that to have bounded solutions we need to require that

$$\psi _i\left(A_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }\psi _1b_2(\chi ,s,\tau _0,\theta _0,ϵ)𝑑s\right)\text{for}i=1,\mathrm{},4,$$
(39)
remains bounded as $`\tau \pm \mathrm{}`$. More precisely $`\stackrel{~}{𝐳}`$ is bounded on $`[\tau _0,\mathrm{})`$ if and only if

$$A=_{\tau _0}^{\mathrm{}}e^{_{s_0}^sTrA(\eta )𝑑\eta }\psi _1b_2𝑑s$$
(40)
and bounded on $`(\mathrm{},\tau _0]`$ if and only if

$$A=_{\mathrm{}}^{\tau _0}e^{_{s_0}^sTrA(\eta )𝑑\eta }\psi _1b_2𝑑s.$$
(41)
We also require

$$\stackrel{~}{u}(\pm \mathrm{})=\underset{\tau \pm \mathrm{}}{lim}\stackrel{~}{u}(\tau )=D+\underset{\tau \pm \mathrm{}}{lim}_{\tau _0}^\tau b_4(\chi ,s,\tau _0,\theta _0,ϵ)𝑑s,$$
(42)
where, obviously, $`\stackrel{~}{u}(\mathrm{})=\stackrel{~}{u}(\mathrm{})`$. The latter condition is not needed for the boundedness of the solution, but its role will be clear later when analyzing some properties of the negatively and positively asymptotic sets. It is easy to see that the above conditions are simultaneously satisfied both at $`\tau =\mathrm{}`$ and at $`\tau =+\mathrm{}`$ if for some $`\tau _0`$ the following Melnikov-type conditions:

$$\begin{array}{cc}& _{\mathrm{}}^+\mathrm{}e^{_{\tau _0}^\tau TrA(\eta )𝑑\eta }\psi _1b_2(\chi ,s,\tau _0.\theta _0,ϵ)ds=0\hfill \\ & _{\mathrm{}}^+\mathrm{}b_4(\chi ,s,\tau _0,\theta _0,ϵ)𝑑s=0\hfill \end{array}$$
(43)
are fulfilled. Thus we can rewrite the general solution (38) using (43) and, by neglecting to mention the dependence on $`\chi `$, $`s`$, etc., we obtain

$$\begin{array}{cc}\hfill \stackrel{~}{r}& =R^{}_{\mathrm{}}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }\psi _1b_2𝑑s+\psi _1\left(B+_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }R^{}b_2𝑑s\right)\hfill \\ \hfill \stackrel{~}{v}& =V^{}_{\mathrm{}}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }\psi _1b_2𝑑s+\psi _2\left(B+_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }R^{}b_2𝑑s\right)\hfill \\ \hfill \stackrel{~}{\theta }& =\sqrt{2b}_{\mathrm{}}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }\psi _1b_2𝑑s+\psi _3\left(B+_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }R^{}b_2𝑑s\right)\hfill \\ & +C_{\tau _0}^\tau e^{_{s_0}^sTrA(\eta )𝑑\eta }(R^{}\psi _3\pm \sqrt{2b}\psi _1)b_2\hfill \\ \hfill \stackrel{~}{u}& =\stackrel{~}{u}(\mathrm{})+_{\mathrm{}}^\tau b_4𝑑s.\hfill \end{array}$$
(44)
To obtain $`\stackrel{~}{r}(\tau _0)R^{}(\tau _0)+\stackrel{~}{v}(\tau _0)V^{}(\tau _0))=0`$ we must have

$$B=\frac{(R^2(\tau _0)+V^2(\tau _0))}{\psi _1(\tau _0)R^{}(\tau _0)+\psi _2(\tau _0)V^{}(\tau _0)}_{\mathrm{}}^{\tau _0}\psi _1(s)b_2(s)𝑑s.$$
(45)
Moreover we also get

$$C=\stackrel{~}{\theta }(\tau _0)\pm \sqrt{2b}_{\mathrm{}}^{\tau _0}e^{_{s_0}^sTrA(\eta )𝑑\eta }\psi _1b_2𝑑s$$
(46)
and

$$D=\stackrel{~}{u}(\mathrm{})+_{\mathrm{}}^{\tau _0}b_4𝑑s.$$
(47)
This uniquely defines $`B`$, $`C\stackrel{~}{\theta }(\tau _0)`$, and $`D\stackrel{~}{u}(\mathrm{})`$ as continuous linear functionals on $`()`$. From (44) we observe that the corresponding solution is of the form $`\stackrel{~}{𝐳}=𝐋(\stackrel{~}{𝐛})+𝐰`$, where $`𝐋`$ is a bounded linear operator. It follows that this operator is continuous and hence the solution $`\stackrel{~}{𝐳}=𝐋(\stackrel{~}{𝐛})+𝐰`$ is continuous on $`()`$. This completes the proof.

To obtain necessary and sufficient conditions that the negatively and positively asymptotic sets intersect, let us first consider all the solutions of (25) which are bounded as $`\tau \mathrm{}`$ and such that their angles remain close to the ones on the periodic orbit. The solution $`\stackrel{~}{𝐳}`$ is given by (38) satisfying (41) and (42) with negative sign. In particular the solutions of the variational equation that are bounded as $`\tau \mathrm{}`$ (i.e. which remain in a sufficiently small neighborhood of the periodic orbit as $`\tau \mathrm{}`$) and with perturbed angles that do not drift but remain near the angles on the periodic orbit, must be on the negatively asymptotic set. In the same way, we obtain the positively invariant set from the solutions that remain bounded as $`\tau \mathrm{}`$ and whose angles stay close to the one of the periodic orbit, which was in fact the reason why we required that condition (42) be satisfied.

Moreover it is important to remark that the solution we found are not only bounded but also such that $`\stackrel{~}{r}0`$, $`\stackrel{~}{v}0`$ as $`\tau \mathrm{}`$ and this is important since, on the collision manifold we have many periodic orbit and this condition is needed to show that the orbits are actually asymptotic to the equator.

With the preparations above, we can now prove the following result.

###### Theorem 3

System (16) has transversal homoclinic solutions if and only if there exist $`\tau _0^{}`$ and a $`\theta _0^{}`$ such that

$$\stackrel{~}{M}_1(\tau _0^{},\theta _0^{})=\stackrel{~}{M}_2(\tau _0^{},\theta _0^{})=0\text{and}\frac{\stackrel{~}{M}_1}{\tau _0}\frac{\stackrel{~}{M}_2}{\theta _0}\frac{\stackrel{~}{M}_1}{\theta _0}\frac{\stackrel{~}{M}_2}{\tau _0}|{}_{\begin{array}{c}\tau _0=\tau _0^{}\\ \theta _0=\theta _0^{}\end{array}}{}^{}0,$$
(48)
where

$$\begin{array}{cc}\hfill \stackrel{~}{M}_1(\tau _0,\theta _0)& =_{\mathrm{}}^+\mathrm{}e^{_{\tau _0}^\tau TrA(s)𝑑s}R^{}b_2(\stackrel{~}{𝐳}^{},\tau ,\tau _0,\theta _0,ϵ)𝑑\tau ,\hfill \\ \hfill \stackrel{~}{M}_2(\tau _0,\theta _0)& =_{\mathrm{}}^+\mathrm{}b_4(\stackrel{~}{𝐳}^{},\tau ,\tau _0,\theta _0,ϵ)𝑑\tau ,\hfill \end{array}$$
(49)
and $`\stackrel{~}{𝐳}^{}`$ is a solution of $`\stackrel{~}{𝐳}=𝐋(\stackrel{~}{𝐛}(\stackrel{~}{𝐳},\tau ,\tau _0,\theta _0,ϵ))+𝐰`$. Moreover if the perturbation is periodic we get infinitely many intersections.

Proof: The stable and unstable manifolds intersect if and only if the solution (38) satisfies the Melnikov-like conditions (35) and (42) of Lemma 2. This was already proved in the case when $`\stackrel{~}{𝐛}`$ did not implicitly depend depend on $`\stackrel{~}{𝐳}`$. But because of this implicit dependence we need to apply the implicit function theorem, which states that given $`\stackrel{~}{𝐳}=𝐋(\stackrel{~}{𝐛}(\stackrel{~}{𝐳},\tau ,\tau _0,\theta _0,ϵ))+𝐰`$ with $`\stackrel{~}{𝐳}𝐰=𝐋(0,\tau ,\tau _0,\theta _0,0)=0`$, there exist a $`\delta `$ and a unique solution $`\stackrel{~}{𝐳}^{}(ϵ,\tau _\mathrm{𝟎},\theta _\mathrm{𝟎})`$ (that has continuous derivatives up to order 2 in $`\tau _0,\theta _0,ϵ`$) such that $`ϵ<0`$, $`|\stackrel{~}{𝐳}|<\delta `$ if the linearized operator $`\stackrel{~}{𝐳}=𝐋(\stackrel{~}{𝐛}(0,\tau ,\tau _0,\theta _0,ϵ))+𝐰`$ is invertible. But Lemma 2 proved that such an operator is invertible. Moreover the homoclinic solutions are transversal if and only if the integrals (49) have simple zeroes, as functions of $`\tau _0`$ and $`\theta _0`$ (see ). This concludes the proof.

Unfortunately the Melnikov integrals of Theorem 3 are difficult to compute explicitly. To overcome this difficulty we need to rewrite these integrals to the first order approximation in $`ϵ`$. Hence if we let $`\stackrel{~}{𝐳}^{}=ϵ\chi `$ and $`\stackrel{~}{𝐛}=ϵ𝐝`$ with $`𝐝=(d_1,d_2,d_3,d_4)`$, the next result follows immediately.

###### Corollary 1

System (16) has transversal homoclinic solutions if and only if there exist $`\tau _0^{}`$ and a $`\theta _0^{}`$ such that

$$M_1(\tau _0^{},\theta _0^{})=M_2(\tau _0^{},\theta _0^{})=0\text{and}\frac{M_1}{\tau _0}\frac{M_2}{\theta _0}\frac{M_1}{\theta _0}\frac{M_2}{\tau _0}|{}_{\begin{array}{c}\tau _0=\tau _0^{}\\ \theta _0=\theta _0^{}\end{array}}{}^{}0,$$
(50)
where

$$\begin{array}{cc}\hfill M_1(\tau _0,\theta _0)& =_{\mathrm{}}^+\mathrm{}e^{_{\tau _0}^\tau TrA(s)𝑑s}R^{}(\tau \tau _0)b_2(\chi (\tau \tau _0),\mathrm{\Theta }(\tau \tau _0)\theta _0)𝑑\tau ,\hfill \\ \hfill M_2(\tau _0,\theta _0)& =_{\mathrm{}}^+\mathrm{}b_4(\chi (\tau \tau _0),\mathrm{\Theta }(\tau \tau _0)\theta _0)𝑑\tau .\hfill \end{array}$$
(51)
Moreover if the perturbation is periodic we get infinitely many intersections.

Corollary 1 generalizes the Melnikov integrals obtained in  to nonhyperbolic whiskered tori (periodic orbits) in non-Hamiltonian systems. We remark that the second integral in (51) converges only conditionally. This is not a new feature of this non-Hamiltonian system since the same nuisance was present in . However some authors, more recently, found a way to write the Melnikov conditions for hyperbolic whiskered tori in Hamiltonian systems using only convergent integrals see . It would be interesting to generalize those results to nonhyperbolic tori in non-Hamiltonian systems and to apply the newly developed technique to the problem under discussion in this paper. But this is not a project we aim to develop here.

## 6 The Melnikov Integrals

Now we would like to apply Corollary 1 to our problem. The Melnikov conditions take the form

$$\begin{array}{c}\hfill M_1(\tau _0,\theta _0)=_{\mathrm{}}^+\mathrm{}[e^{\frac{1}{2}_{\tau _0}^\tau V(s)𝑑s}R(\tau \tau _0)R^{}(\tau \tau _0)\\ \hfill \times \mathrm{cos}^2(\omega (\tau \tau _0)\theta _0)]d\tau =0\end{array}$$
(52)
and

$$M_2(\tau _0,\theta _0)=\frac{1}{2}_{\mathrm{}}^+\mathrm{}(R(\tau \tau _0)+2b)\mathrm{sin}(2(\omega (\tau \tau _0)\theta _0))𝑑\tau =0.$$
(53)
Let $`\stackrel{~}{\theta }_0=\theta _0\omega \tau _0`$. With this assumption we can rewrite the first Melnikov condition as

$$M_1=\mathrm{cos}^2\stackrel{~}{\theta }_0I_1^a+\mathrm{sin}^2\stackrel{~}{\theta }_0I_1^b\mathrm{sin}2\theta _0I_1^c,$$
(54)
where

$$\{\begin{array}{c}I_1^a=_{\mathrm{}}^+\mathrm{}e^{\frac{1}{2}_{\tau _0}^\tau V(s)𝑑s}RR^{}\mathrm{cos}^2\omega \tau d\tau \hfill \\ I_1^b=_{\mathrm{}}^+\mathrm{}e^{\frac{1}{2}_{\tau _0}^\tau V(s)𝑑s}RR^{}\mathrm{sin}^2\omega \tau d\tau \hfill \\ I_1^c=_{\mathrm{}}^+\mathrm{}e^{\frac{1}{2}_{\tau _0}^\tau V(s)𝑑s}RR^{}\mathrm{sin}\omega \tau \mathrm{cos}\omega \tau d\tau .\hfill \end{array}$$
(55)
The second Melnikov condition can be expressed as

$$M_2=\mathrm{cos}2\stackrel{~}{\theta }_0I_2^a+\mathrm{sin}2\stackrel{~}{\theta }_0I_b^2,$$
(56)
where

$$\{\begin{array}{c}I_2^a=\frac{1}{2}_{\mathrm{}}^+\mathrm{}(R+2b)\mathrm{sin}2\omega \tau d\tau \hfill \\ I_2^b=\frac{1}{2}_{\mathrm{}}^+\mathrm{}(R+2b)\mathrm{cos}2\omega \tau d\tau .\hfill \end{array}$$
(57)
All the integrals above can be computed using the method of residues. Straightforward computations give

$$I_1^a=I_1^b=\frac{1}{|h|}_{\mathrm{}}^+\mathrm{}\frac{(\tau \tau _0)\mathrm{cos}2\omega \tau }{(2|h|+(\tau \tau _0)^2)^2}𝑑\tau =\frac{\pi \mathrm{sin}(2\omega \tau _0)e^{2\omega \sqrt{2|h|}}}{|h|\sqrt{2|h|}}$$
(58)
and

$$I_1^c=\frac{1}{|h|}_{\mathrm{}}^+\mathrm{}\frac{(\tau \tau _0)\mathrm{sin}\omega \tau \mathrm{cos}\omega \tau }{(2|h|+(\tau \tau _0)^2)^2}𝑑\tau =\frac{\pi \mathrm{cos}(2\omega \tau _0)e^{2\omega \sqrt{2|h|}}}{|h|\sqrt{2|h|}}.$$
(59)
Particular care is needed when integrating $`I_2^a`$ and $`I_2^b`$ since they converge only conditionally. To obtain computational convergence, we choose the limits in $`I_2^a`$ such that

$$\begin{array}{cc}\hfill I_2^a& =\underset{N\mathrm{}}{lim}_{N\pi /2\omega }^{N\pi /2\omega }\left(b+\frac{1}{2|h|+(\tau \tau _0)^2}\right)\mathrm{sin}2\omega \tau d\tau \hfill \\ & =\frac{\pi \mathrm{sin}(2\omega \tau _0)e^{2\omega \sqrt{2|h|}}}{\sqrt{2|h|}}.\hfill \end{array}$$
(60)
The integral was also computed using the method of residues. Similarly, for $`I_2^b`$, we have

$$\begin{array}{cc}\hfill I_2^b& =\underset{N\mathrm{}}{lim}_{N\pi /2\omega }^{N\pi /2\omega }\left(b+\frac{1}{2|h|+(\tau \tau _0)^2}\right)\mathrm{cos}2\omega \tau d\tau \hfill \\ & =\frac{\pi \mathrm{cos}(2\omega \tau _0)e^{2\omega \sqrt{2|h|}}}{\sqrt{2|h|}}\hfill \end{array}$$
(61)
and thus

$$M_1=M_2=\mathrm{sin}(2(\omega \tau _0+\stackrel{~}{\theta }_0))\frac{\pi e^{2\omega \sqrt{2|h|}}}{\sqrt{2|h|}}.$$
(62)
We therefore have only one independent condition; this is clearly a consequence of the energy relation.

We can find simple zeroes when $`\mathrm{sin}(2(\omega \tau _0+\stackrel{~}{\theta }_0))=0`$, i.e., for $`(\omega \tau _0+\stackrel{~}{\theta }_0)=\theta _0=\pm k\pi /2`$ for $`k=0,1,2,\mathrm{}`$.

Hence, by Corollary 1, we have proved the existence of an infinite sequence of intersections on the Poincaré section of the negatively and positively asymptotic sets of the periodic orbit and the existence of homoclinic orbits leaving the equator of the collision manifold and going back to it. This situation is clearly reminiscent of the chaotic dynamics described by the Poincaré-Birkhoff-Smale theorem in terms of symbolic dynamics and the Smale horseshoe. Unfortunately this theorem cannot be directly applied, nor can the theorems proved in , since the Poincaré-Birkhoff-Smale theorem considers hyperbolic fixed points while the arguments in  apply to area-preserving diffeomorphisms. However the arguments contained in those theorems strongly suggest the occurrence of a chaotic dynamics.

Moreover it is easy to verify, and interesting to remark, that the orbits we found above are not $`\overline{S}_0`$-symmetric, where the $`\overline{S}_0`$ symmetry is defined by $`\overline{S}_0(r,v,\theta ,u,\tau )=(r,v,u,\tau )`$ (see ) and an orbit $`\gamma (\tau )`$ is said to be $`\overline{S}_0`$-symmetric if $`\overline{S}_0(\gamma (\tau ))=\gamma (\tau )`$. Indeed an orbit is $`\overline{S}_0`$-symmetric if and only if it has a point on the zero velocity curve, i.e., if there is a $`\overline{\tau }`$ such that $`v(\overline{\tau })=u(\overline{\tau })=0`$ (see ). But this cannot happen in our problem because the unperturbed solution verifies $`u\pm \sqrt{2b}`$. Thus for $`ϵ`$ small enough the perturbed orbit can never have $`u=0`$.

We can now summarize the above discussion as follows:

###### Theorem 4

Let us consider the anisotropic Manev problem given by the equation of motion (6) with the energy relation (7). Then there is an infinite sequence of intersections in the Poincaré section of the negatively and positively asymptotic sets of the periodic orbits at the equator of the collision manifold (possibly giving rise to a chaotic dynamics). Furthermore there exist the homoclinic non $`\overline{S}_0`$-symmetric orbits to the periodic orbit described above.

## 7 Periodic Solutions

We now return to the original Cartesian coordinates, which are more convenient for the purpose of finding certain periodic solutions. Let us first notice that the equations (2) admit the following symmetries:

$$\begin{array}{c}S_0(x,y,p_x,p_y,t)=(x,y,p_x,p_y,t),\hfill \\ S_1(x,y,p_x,p_y,t)=(x,y,p_x,p_y,t),\hfill \\ S_2(x,y,p_x,p_y,t)=(x,y,p_x,p_y,t),\hfill \\ S_3(x,y,p_x,p_y,t)=(x,y,p_x,p_y,t),\hfill \\ S_4(x,y,p_x,p_y,t)=(x,y,p_x,p_y,t),\hfill \\ S_5(x,y,p_x,p_y,t)=(x,y,p_x,p_y,t),\hfill \\ S_6(x,y,p_x,p_y,t)=(x,y,p_x,p_y,t),\hfill \end{array}$$

which are the elements of an Abelian group of order eight, isomorphic to $`𝐙_2\times 𝐙_2\times 𝐙_2`$, that is generated by $`S_0,S_1,S_2`$ (see ). (The symmetry $`S_0`$ is the one denoted by $`\overline{S}_0`$ in the McGehee coordinates of the previous section.) To obtain certain families of periodic solutions, we will use the symmetries $`S_0`$,$`S_1`$ and $`S_2`$ in connection with the variational principle according to which extremum values of the action integral yield periodic solutions of the equations (2). To reach this goal we first need to introduce some notations.

Let $`C^{\mathrm{}}([0,T],^2)`$ be the space of $`T`$-periodic $`C^{\mathrm{}}`$ cycles $`f:[0,T]^2`$. Define the inner products

$$\begin{array}{c}f,g_{L^2}=_0^Tf(t)g(t)𝑑t,\hfill \\ f,g_{H^1}=f,g_{L^2}+\dot{f},\dot{g}_{L^2},\hfill \end{array}$$
(63)
and let $`_{L^2}`$, $`_{H^1}`$ be the corresponding norms. Then the completion of $`C^{\mathrm{}}([0,T],^2)`$ with respect to the norm $`_{L^2}`$ is denoted by $`L^2`$ and it is the space of square integrable functions. The completion with respect to $`_{H^1}`$ is denoted by $`H^1`$ and is the Sobolev space of all absolutely continuous $`T`$-periodic paths that have $`L^2`$ derivatives defined almost everywhere (see ).

Let $`\mathrm{\Sigma }_i([0,T],^2)`$ denote the subset of $`H^1`$ formed by the $`S_i`$-symmetric paths, with $`i\{0,1,2,3,4,5,6\}`$. It is easy to see that each $`\mathrm{\Sigma }_i`$ is a subspace of $`H^1`$; in fact they are Sobolev spaces and have many interesting properties. In the following we will restrict our attention to the spaces $`\mathrm{\Sigma }_i`$ with $`i=0,1,2,6`$. Let us now prove the following result.

###### Lemma 3

Let $`H^1`$ be defined as above, then the subspaces $`\mathrm{\Sigma }_i`$ of $`S_i`$-symmetric paths with $`i=0,1,2,6`$ are closed, weakly closed, and complete with respect to the norm $`_{H^1}`$, and are therefore Sobolev spaces. Moreover

$$H^1=\mathrm{\Sigma }_1\mathrm{\Sigma }_2=\mathrm{\Sigma }_0\mathrm{\Sigma }_6.$$
(64)
Proof: We first show an interesting fact: we can write $`f=(f^1,f^2)`$ as the sum of an $`S_1`$ and an $`S_2`$-symmetric path. Indeed it is well known that we can write $`f_1`$ and $`f_2`$ as the sum of an even and an odd absolutely continuous function, i.e. as $`f_1=f_1^e+f_1^o`$ and $`f_2=f_2^e+f_2^o`$. Using this idea we can write the path $`f(t)`$ as the sum of an $`S_1`$-symmetric function, $`f_{S_1}=(f_1^e,f_2^o)`$, and an $`S_2`$-symmetric one, $`f_{S_2}=(f_1^o,f_2^e)`$. Now fix an element $`f\mathrm{\Sigma }_1`$. Then $`f,g_{H^1}=0`$ for every $`g\mathrm{\Sigma }_2`$. This is because

$$f,g_{H^1}=_0^T(f_1g_1+f_2g_2)𝑑t+_0^L(f_1^{}g_1^{}+f_2^{}g_2^{})𝑑t,$$

where the first integrand is an odd function and the second is an odd function almost everywhere. Thus the above scalar product is zero for every $`g\mathrm{\Sigma }_2`$.

Let us denote the space orthogonal to $`\mathrm{\Sigma }_1`$ by $`\mathrm{\Sigma }_1^{}=\{g\mathrm{\Sigma }_1:f,g_{H^1}=0\text{for every}g\mathrm{\Sigma }_1\}`$. It is easy to see that $`\mathrm{\Sigma }_1^{}`$ is closed and that $`S_2\mathrm{\Sigma }_1^{}`$. Now we need to show that $`S_2\mathrm{\Sigma }_1^{}`$. Assume there is $`h\mathrm{\Sigma }_1^{}`$ such that $`h0`$ and $`h\mathrm{\Sigma }_2`$. Then write $`h=h_{S_1}+h_{S_2}`$ and consider $`h_{S_1},h_{S_1}+h_{S_2}_{H^1}`$, which means that $`h_{S_1},h_{S_1}_{H^1}=h_{S_1}=\delta >0`$. But this contradicts the hypothesis that $`h\mathrm{\Sigma }_1^{}`$. Therefore $`\mathrm{\Sigma }_2=\mathrm{\Sigma }_1^{}`$. So $`\mathrm{\Sigma }_2`$ and consequently $`\mathrm{\Sigma }_1`$ are closed and such that $`H^1=\mathrm{\Sigma }_1\mathrm{\Sigma }_2`$. Moreover, since $`H^1`$ is a metric space, $`\mathrm{\Sigma }_1`$ and $`\mathrm{\Sigma }_2`$ are complete. Also $`\mathrm{\Sigma }_1`$ and $`\mathrm{\Sigma }_2`$ are weakly closed since they are norm-closed subspaces. The statements for $`\mathrm{\Sigma }_0`$ and $`\mathrm{\Sigma }_6`$ can be proved in a similar way. This completes the proof.

Let us now introduce some new definitions. We will say that a path in $`\mathrm{\Sigma }_i`$ is of class $`L_n`$, $`n=0,\pm 1,\pm 2,\pm 3,\mathrm{}`$, if its winding number about the origin of the coordinate system is $`n`$ (i.e. if it makes $`n`$ loops around the origin). The sign of $`n`$ is positive for a counterclockwise rotation and negative otherwise. Consider the sets $`\overline{\mathrm{\Sigma }}_i([0,T],^2\backslash \{0\})`$. Notice that they are open submanifolds of the spaces $`\mathrm{\Sigma }_i([0,T],^2)`$ and that the family $`(L_n)_{n𝐙}`$ provides a partition of those spaces into homotopy classes, also called components. Two periodic orbits of the isotropic Manev problem ($`\mu =1`$), one of class $`L_8`$ and the other of class $`L_9`$, are depicted in Fig. 4.

The Lagrangian $`L(𝐪,\dot{𝐪})=T(\dot{𝐪})+U(𝐪)`$ of the anisotropic Manev problem given by system (2) has the expression

$$L(x,y,\dot{x},\dot{y})=\frac{1}{2}(\dot{x}^2+\dot{y}^2)+\frac{1}{\sqrt{x^2+\mu y^2}}+\frac{b}{x^2+\mu y^2},$$
(65)
and the action integral along a path $`f`$ from time $`0`$ to time $`T`$, whose Euclidean coordinate representation is $`𝐪=𝐪(t)=(x(t),y(t))`$, takes the form

$$A_T(f)=_0^TL(𝐪(t),\dot{𝐪}(t))𝑑t.$$

According to Hamilton’s principle, the extremals of the functional $`A_T`$ are solutions of the equations (2). Hence we want to obtain periodic solutions of (2) by finding extremals of the functional $`A`$. For this we will use a direct method of the calculus of variation, namely the lower-semicontinuity method (see ). In preparation of a satisfactory theory of existence, the notion of admissible function has to be relaxed since the extremals we obtain belong to a Sobolev space. Therefore the above method provides only “weak” solutions of our problem. To show that the paths are regular enough to be classical solutions, we need the following result, proved in .

###### Lemma 4

The critical points of $`A_T|_{\overline{\mathrm{\Sigma }}_i([0,T],^2\backslash \{0\})}`$ are $`T`$-periodic solutions of equations (2).

In particular it is well known that if $`f`$ is a minimizer of the action $`A_T`$ in the space $`H^1([0,t],^2)`$ and if $`f`$ has no collisions, then $`f`$ is a $`T`$-periodic solution to (2). Collision have to be excluded because equations (2) break down at collisions and because the action is not differentiable at paths with collisions. In this paper we are interested to restrict ourself to the spaces $`\mathrm{\Sigma }_i`$ of $`S_i`$-symmetric paths for $`i=0,1,2`$. The paths that are $`S_6`$-symmetric have to be excluded in the study of periodic orbits since $`S_6`$-symmetric paths must intersect the origin and therefore encounter collisions.

Now it is not obvious that a collisionless minimizer in $`\mathrm{\Sigma }_i`$ is a periodic solution of system (2). However, according to the principle of “symmetric criticality” (see for example ) this is actually true. Indeed, it can be proved that if $`f`$ is a collision free path with $`dA_t(f)(h)=0`$ for every $`f\mathrm{\Sigma }_i`$, then $`dA_T(f)(h)=0`$ for all $`fH^1([0,T],^2)`$ and thus $`f`$ is a critical point in the bigger loop space $`H^1`$ (see ).

The only obstacle left for applying the direct method is the “noncompactness” of the configuration space. Indeed we want to exclude the possibility that the minimizer is obtained when the bodies are at infinite distance from each other or are collision paths. The first problem is solved restricting ourselves to non-simple cycles, i.e., to cycles that are not homotopic to a point and thus are not in the homotopy class $`L_0`$. The second problem is solved by the following result.

###### Lemma 5

Any family $`\mathrm{\Gamma }`$ of non-simple homotopic cycles in $`\overline{\mathrm{\Sigma }}_i([0,T],^2\backslash \{0\})`$ for $`i=0,1,2`$ on which $`J(f)=_0^T\frac{1}{2}|\dot{𝐪}(t)|^2𝑑t`$ and $`E(f)=_0^TU(𝐪(t))𝑑t`$ are bounded, is bounded away from the origin.

The proof of this result follows from  if we remark that the anisotropic Manev potential is “strong” according to Gordon’s definition and that the Lagrangian is positive.

To apply the direct method we still need to recall some properties of lower semicontinuous (l.s.c) functions. Let $`:X`$ be a real valued function on a topological space $`X`$. Then $``$ is l.s.c. if and only if $`^1(\mathrm{},a]`$ is closed for every $`a`$, in which case $``$ is bounded below and attains its infimum on every compact subset of $`X`$. Moreover when $`X`$ is Hausdorff then compact sets are necessarily closed and thus we have the following result.

###### Proposition 2

Suppose $`:X`$ is a real valued function on an Hausdorff space $`X`$ and

$$^1(\mathrm{},b]\text{ is compact for every real b}.$$

Then $``$ is l.s.c., bounded below, and attains its infimum value on $`X`$.

We can now prove the main result of this section.

###### Theorem 5

For any $`T>0`$ and any $`n=\pm 1,\pm 2,\pm 3,`$, there is at least one $`S_i`$-symmetric ($`i=0,1,2`$) periodic orbit of the anisotropic Manev problem that has period $`T`$ and winding number $`n`$ (i.e., belongs to the homotopy class $`L_n`$).

Proof: Let $`X`$ be a component of $`\overline{\mathrm{\Sigma }}_i([0,T],^2\backslash \{0\})`$ for $`i=0,1,2`$, that consist of non-simple cycles. Endow $`X`$ with the weak topology it inherits from $`\mathrm{\Sigma }_i([0,T],^2)`$. Then $`X`$ is a subset of an Hilbert space and it is weakly compact if and only if it is weakly closed.

We wish to apply Proposition 2 with $`=A_T`$ and thus we have to show that $`XA_T^1(\mathrm{},b]`$ is a bounded and weak-closed subset of $`\mathrm{\Sigma }_i([0,T],^2)`$.

Since $`J=A_TE`$ and $`U>0`$, we have $`E>0`$ and therefore

$$\begin{array}{cc}& Jb\text{on}A_{T}^{}{}_{}{}^{1}(\mathrm{},b]=A_{T}^{}{}_{}{}^{1}[0,b],\hfill \\ & E=A_TJb\text{on}A_{T}^{}{}_{}{}^{1}[0,b].\hfill \end{array}$$
(66)
Since $`Jb`$ the elements of $`X`$ are bounded in arc length, and from Lemma 5 it follows that the elements of $`X`$ are bounded away from the origin. Moreover the elements of $`X`$ are non-simple and thus bounded in the $`C^0`$ norm and hence in the $`L^2`$ norm. This last fact combined with $`Jb`$ shows that $`X`$ is bounded in the $`_{H^1}`$ norm. Thus also $`XA_{T}^{}{}_{}{}^{1}(\mathrm{},b]`$ is bounded in the $`H^1`$ norm.

Now suppose that $`\{f_n\}=\{(f_n^1,f_n^2)\}`$ is a sequence in $`XA_{T}^{}{}_{}{}^{1}[0,b]`$ that converges weakly to a cycle $`f\mathrm{\Sigma }_i([0,T],^2)`$ for $`i=0,1,2`$. From general principles, $`f_n_{H^1}`$ is bounded and $`f_n_{L^2}f_n_{L^2}`$ because weak $`\mathrm{\Sigma }_i`$-convergence implies $`C^0`$-convergence. Since $`J(f_n)=1/2f_n_{H^1}^21/2f_n_{L^2}^2`$ it means that $`J(f_n)`$ is bounded and since $`Eb`$ on $`A_{T}^{}{}_{}{}^{1}[0,b]`$ it follows that $`\{E(f_n)\}`$ is bounded. Moreover, Lemma 5 guarantees that the functions $`f_n`$ are bounded away from the origin so that $`f`$ is homotopic to the $`f_n`$ in $`^2\backslash \{0\}`$. Therefore $`fX`$.

To complete the proof we have to show that $`fA_{T}^{}{}_{}{}^{1}[0,b]`$. We know that $`E(f_n)E(f)`$ since weak convergence in $`\mathrm{\Sigma }_i`$ implies $`C^0`$-convergence. For each $`n`$ let

$$g_n(t)=\frac{1}{\sqrt{(f_n^1(t))^2+\mu (f_n^2(t))^2}}+\frac{1}{(f_n^1(t))^2+\mu (f_n^2(t))^2}$$

and denote

$$g(t)=\frac{1}{\sqrt{(f^1(t))^2+\mu (f^2(t))^2}}+\frac{1}{(f^1(t))^2+\mu (f^2(t))^2}.$$

Each $`g_n`$ is of class $`L^1`$ since $`A_T(f_n)<\mathrm{}`$. This implies that the set of all $`t`$ for which $`f_n(t)=0`$ has zero measure, otherwise the integral of $`g_n(t)`$ would be unbounded. So $`g_n(t)g(t)`$ almost everywhere. Also $`_0^Tg_n(t)𝑑t<A_T(f_n)b`$. By Fatou’s lemma it follows that $`g`$ is $`L^1`$ and that

$$_0^Tg(t)𝑑t=_0^Tlim\; infg_n(t)dtlim\; inf_0^Tg_n(t)𝑑t.$$

Now we can use the fact that the norm is weakly sequentially lower semicontinuous (see ), thus

$$\dot{f}_{L^2}^2=f_{H^1}^2f_{L^2}^2lim\; inff_n_{H^1}^2f_{L^2}^2=lim\; inf\dot{f}_n_{L^2}^2,$$

where the last equality holds since $`\{f_n\}`$ converges strongly to $`f`$ in $`L^2`$. Consequently

$$\begin{array}{cc}\hfill A_T(f)& =\frac{1}{2}\dot{f}_{L^2}^2+_0^Tg(t)𝑑t\hfill \\ & lim\; inf\frac{1}{2}\dot{f}_n_{L^2}^2+lim\; inf_0^Tg_n(t)𝑑tlim\; infA_T(f_n)b.\hfill \end{array}$$
(67)
Relation (67) now implies that $`fA_{T}^{}{}_{}{}^{1}[0,b]`$. This completes the proof.

Recall now that two intersections of every $`S_1`$-symmetric ($`S_2`$-symmetric) orbit with the $`x`$ axis ($`y`$ axis) must be orthogonal. To distinguish them from accidental orthogonal intersections, which do not follow because of the symmetry, we will call them essential orthogonal intersections. From the proof of Theorem 5 and obvious index theory considerations, the following result follows (see also Fig. 4).

###### Corollary 2

If the essential orthogonal intersections with the $`x`$-axis ($`y`$-axis) of an $`S_1`$-symmetric ($`S_2`$-symmetric) periodic orbit lie on the same side of the axis with respect to the origin of the coordinate system, then the orbit has an even winding number. If the essential orthogonal intersections are on opposite sides with respect to the origin, then the periodic orbit has an odd winding number.

Since the symmetries $`S_0,S_1`$ and $`S_2`$ generate the entire symmetry group, it is clear that Theorem 5 captures all periodic orbits with symmetries. This result, however, does not tell if other symmetric periodic orbits exist beyond the the ones with $`S_0`$, $`S_1`$ and $`S_2`$ symmetries. Let us therefore end our paper by proving that $`S_3`$-periodic orbits do indeed exist. In fact they form a rich set if compared to the one of $`S_3`$-symmetric orbits of the anisotropic Kepler problem (given by (1) with $`b=0`$), which contains only circular orbits. We will show that in our case each homotopy class $`L_n`$, $`n=4k+1`$, $`k`$ integer, contains at least one $`S_3`$-symmetric periodic orbit. Other homotopy classes may have $`S_3`$-symmetric periodic orbits, but our approach proves their existence only for winding numbers of the form $`n=4k+1`$, $`k`$ integer.

We consider the set of all paths with one end on the $`x`$ axis and the other on the $`y`$ axis of the coordinate system. As in the case of periodic cycles discussed in the first part of this section, for a given $`T^{}=T/4>0`$ this set can be endowed with a Hilbert space structure, the completion of which is a Sobolev space. We further divide this space in homotopy classes $`_n,n=0,\pm 1,\pm 2,\mathrm{}`$ according to the winding number $`n`$.

Using the boundary conditions, it is easy to see that in each class $`_n`$ the minimizer of the action is a an arc orthogonal to the $`x`$ and $`y`$ axes. Its existence and the fact that it is a solution in the classical sense can be proved in a similar way as we did for periodic cycles. Once obtaining such a solution with ends on the $`x`$ and $`y`$ axes, we can use the $`S_3`$-symmetry and the orthogonality with the axes to complete this solution arc to a periodic orbit of period $`T>0`$. The symmetry implies that the winding number is of the form $`n=4k+1`$, $`k`$ integer. This is because if, for example, a solution arc with the ends on the $`x`$ and $`y`$ axes has a loop around the origin, then the corresponding periodic orbit has four loops around the origin. We have thus obtained the following result.

###### Theorem 6

For any $`T>0`$ and any $`n=4k+1`$, $`k`$ integer, there is at least one $`S_3`$-symmetric periodic orbit of the anisotropic Manev problem that has period $`T`$ and winding number $`n`$ (i.e., belongs to the homotopy class $`L_n`$).

It is interesting to note in conclusion that if viewing the anisotropy parameter as a perturbation and the anisotropic Manev problem as a perturbation of the isotropic case (see Section 4), the $`S_i`$-symmetric ($`i=0,1,2`$) periodic orbits of the isotropic problem are deformed but not destroyed by introducing the anisotropy, no matter how large its size. This shows that the $`S_i`$ symmetries ($`i=0,1,2`$) play an important role in understanding the system and are an indicator of its robustness relative to perturbations.

Acknowledgements. Florin Diacu was supported in part by the Pacific Institute for the Mathematical Sciences and by the NSERC Grant OGP0122045. Manuele Santoprete was sponsored by a University of Victoria Fellowship and a Howard E. Petch Research Scholarship.


format: text

text:
# Berry phase in a non-isolated system

## Abstract

We investigate the effect of the environment on a Berry phase measurement involving a spin-half. We model the spin+environment using a biased spin-boson Hamiltonian with a time-dependent magnetic field. We find that, contrary to naive expectations, the Berry phase acquired by the spin can be observed, but only on timescales which are neither too short nor very long. However this Berry phase is not the same as for the isolated spin-half. It does not have a simple geometric interpretation in terms of the adiabatic evolution of either bare spin-states or the dressed spin-resonances that remain once we have traced out the environment. This result is crucial for proposed Berry phase measurements in superconducting nanocircuits as dissipation there is known to be significant.

It was recently suggested Falci00  that it should be possible to observe the Berry phase (BP)Berry84  in a superconducting nanostructure, and possibly use it to control the evolution of the quantum state Jones00Ekert00 ; Wilhelm01 . This intriguing sugestion however did not consider the coupling to the environment, which is never negligible in such structures Nakamura99 . To truly understand the feasibility of the proposed experiment, we must know the effect of the environment on the BP. Originally the BP was defined for systems whose states were separated by finite energy gaps. Here we ask whether a BP can be observed in a system whose spectrum is continuous because it is not completely isolated from its environment. All real systems are coupled, at least weakly, to their environment and as a result never have a truly discrete energy level spectrum. The usual requirement for adiabaticity is that the parameters of the Hamiltonian are varied slowly compared to the gap in the spectrum. Here there is no gap so naively one would say that adiabaticity is impossible and hence the BP could never be observed. However experiments have observed the BP, both directly and indirectly Anandan97 , so this argument must be too naive. We therefore take a simple model in which a quantum system, which when isolated exhibits a BP, is coupled to many other quantum degrees of freedom. We then ask two questions. Firstly, under what conditions can the BP be observed? Secondly, is the observed BP the same as that of the isolated system? While others have investigated systems with a BP coupled to other degrees of freedom Gaitan98 ; Ao99 ; Avron98Avron99 , we believe we are the first to explicitly address these two questions.

We distinguish between the system and the environment in the following way. We have complete experimental control over the system, but almost no control over the environment. The most that we can do to the environment is to ensure the “universe” (system $`+`$ environment) is in thermal equilibrium, with a temperature $`T`$. We will assume we have enough control over $`T`$ to take it to zero, and thus prepare the universe in its ground state. However any procedure to measure a BP in an isolated system must involve measuring a phase difference from a superposition of two states. When the system is not isolated most such procedures involve the mixing of a large number of eigenstates (of the universe), this leads to the effects that we discuss belowendnote-about-Avron .

We choose to investigate a spin-half which is coupled to both a magnetic field and an environment (a bath of harmonic oscillators). Our model is a biased spin-boson model Leggett87  with a time-dependent field. When isolated from its environment, the spin exhibits a BP if we slowly rotate the magnetic field around a closed loop. This model, chosen primarily for its simplicity, is extremely relevant to a recent proposal for observing a BP in a superconducting nanocircuit Falci00 . While we make no attempt to accurately model the true coupling between the nanocircuit and its environment, we believe our results give an excellent indication of what to expect in the real system. Our work will also be very relevant to realisations of the BP quantum computers proposed in Jones00Ekert00 .

In this Letter we concentrate on an Ohmic environment Leggett87 , with the universe initially at zero-temperature wg-tbp . We find that the spin-environment coupling causes the spin-eigenstates to become spin-resonances which have the following properties. (i) The energy distance between them is Lamb shifted by $`\delta E`$. (ii) The higher energy resonance exponentially decays to the lower one on a time-scale, $`T_1`$, and observables containing phase information exponential decay on a timescale $`T_2`$. (iii) There are adiabatic phase-shifts, which divide into two catagories with different symmetries; the phase which vanishes when the Hamiltonian is time-independent we call $`\delta \mathrm{\Phi }_{\mathrm{BP}}`$; while those phase-shifts (and amplitudes) which do not vanish we schematically refer to as $`\mathrm{\Phi }_{\mathrm{shift}}`$. The former scales with the winding number of the BP experiment, while the latter does not (see below). All of these effects go like the second power of the spin-environment coupling, see eqs. (5)-(7).

Effect (ii) means that one cannot perform an arbitrarily long experiment to measure a phase: so we must find the BP from an experiment where the system’s Hamiltonian is taken round a closed loop in a finite time period, $`t_\mathrm{p}T_2`$. In such an experiment there is typically a non-zero amplitude for returning to the initial state and this amplitude has a phase. We interpret the latter as the sum of a dynamic phase which scales linearly with $`t_\mathrm{p}`$, an adiabatic phase ($`\mathrm{\Phi }_{\mathrm{BP}}+\mathrm{\Phi }_{\mathrm{shift}}`$) which is independent of $`t_\mathrm{p}`$, and non-adiabatic contributions which are proportional to $`\mathrm{\Delta }t_\mathrm{p}`$ to some negative powernon-adiab . Here $`\mathrm{\Delta }`$ is the energy difference between the spin-resonances (we set $`\mathrm{}=1`$). Thus the BP is present for arbitrary $`t_\mathrm{p}`$, it is simply masked by the non-adiabatic contributions unless $`t_\mathrm{p}`$ is long enough. For the BP to be observed we must choose a value for $`t_\mathrm{p}`$ which is neither too short nor very long, so that it obeys $`\mathrm{\Delta }^1t_\mathrm{p}T_2`$. However we then actually observe a combination of $`\mathrm{\Phi }_{\mathrm{BP}}`$ and $`\mathrm{\Phi }_{\mathrm{shift}}`$. To distinguish between these two effects we note that when we do not rotate the Hamiltonian $`\mathrm{\Phi }_{\mathrm{BP}}=0`$ while $`\mathrm{\Phi }_{\mathrm{shift}}`$ is unchanged.

Now we ask if the environment’s effect on the BP is observable. To do this we must first decide what BP we would naively expect to observe. There are two possible cases to consider: (i) The system evolves in a magnetic field that we directly control, then we would expect the BP to be given by the solid-angle enclosed by that field, $`\mathrm{\Phi }_{\mathrm{BP}}^{(0)}`$. The deviation from this expectation is given by $`\delta \mathrm{\Phi }_{\mathrm{BP}}`$ in eq. (6). For this deviation to be observable it must be much larger than the non-adiabatic corrections at $`t_\mathrm{p}T_2`$; this means that $`\mathrm{\Delta }T_2\delta \mathrm{\Phi }_{\mathrm{BP}}1`$. The functional form of $`T_2`$ and $`\delta \mathrm{\Phi }_{\mathrm{BP}}`$, in (5) and (6), have the same dependence on the strength of the coupling to the environment, $`C`$. Thus the condition reduces to one dominated by the dependence on $`\gamma `$, where $`\gamma `$ (defined below eq. (6)) characterises the environment. We conclude that there is a wide range of values of $`\gamma `$ for which we can observe $`\delta \mathrm{\Phi }_{\mathrm{BP}}`$. (ii) The second case is more complicated, but is relevant to the superconducting nanocircuit in Falci00 . There we have no independent measure of the bare spin Hamiltonian, the control parameters (gate voltages and magnetic fluxes) enter the spin Hamiltonian in combination with unknown constants (capacitances and inductances). Thus we know nothing about the bare spin-eigenstates, or the solid angle that they enclose when we vary the experimental parameters. However we can measure the spin resonances in the presence of the environment as a function of the experimental parameters. Then one might predict the observed BP is given by the solid-angle enclosed by these spin-resonances. This prediction is given above (8); it is of a similar form to the correct result, but contains a very different function of the distribution of oscillators in the environment. The deviation from this expectation is given by $`\delta \mathrm{\Phi }_{\mathrm{BP}}^{}`$, for it to be observable we require that $`\mathrm{\Delta }T_2\delta \mathrm{\Phi }_{\mathrm{BP}}^{}1`$. Again this reduces to a function independent of $`C`$, where $`\delta \mathrm{\Phi }_{\mathrm{BP}}^{}`$ is observable over a wide range of $`\gamma `$. Finally, we assume we measure $`\mathrm{\Phi }_{\mathrm{shift}}`$ when for a time-independent Hamiltonian before carrying out the BP experiment. Then we do not require $`\delta \mathrm{\Phi }_{\mathrm{BP}}`$ (or $`\delta \mathrm{\Phi }_{\mathrm{BP}}^{}`$) to be larger than $`\mathrm{\Phi }_{\mathrm{shift}}`$ for it to be observable.

To be concrete we assume here that the BP is measured using the spin-echo methodJones00Ekert00 . We consider an experiment where we start with the field along the $`z`$-axis and the universe (spin+oscillators) is in its ground state Shnirman02 . (a) The field is then (instantaneously) prepared at its initial value, $`𝑩_0`$, which is at angle $`\theta `$ to the $`z`$-axis, At the same time the spin is (intantaneously) placed in the state $`\frac{1}{\sqrt{2}}\left(|+|\right)`$ relative to $`𝑩_0`$. Then (b) we adiabatically rotate the magnetic field, $`𝑩(t)`$, $`n`$ times around a closed loop with constant angular velocity, $`𝝎=\widehat{𝒛}\omega `$, (see Fig. 1) for a time period, $`t_\mathrm{p}=2\pi n/\omega `$. We call $`n`$ the winding number. After which (c) the spin is flipped and (d) the field is rotated with angular velocity $`𝝎`$ for time $`t_\mathrm{p}`$. Finally (e) the spin is flipped again and (f) the spin state is measured. By “flip the spin” we mean $`||`$, where the $``$ and $``$ are relative to the direction of the $`𝑩`$-field at that time. This can be achieved by applying a instantaneous $`\pi `$-pulse oriented along the $`y`$-axis. By instantaneous we mean much faster than the fastest oscillator in the environment. We ask what the probability is that the final spin-state, after carrying out (a)-(f), is in a given direction in the plane perpendicular to $`𝑩_0`$. For an isolated spin-half, the probability of the final spin-state being $`\frac{1}{\sqrt{2}}\left(\mathrm{e}^{\mathrm{i}\xi /2}|+\mathrm{e}^{\mathrm{i}\xi /2}|\right)`$ is Jones00Ekert00

$`P(\xi )=\frac{1}{2}\left[1+\mathrm{cos}\left(\xi 4\mathrm{\Phi }_{\mathrm{BP}}^{(0)}\right)\right],`$ (1)
where all spin-states are defined relative to the axis of the field $`𝑩_0`$. Measuring this probability as a function of $`\xi `$ yields the BP for an isolated spin, $`\mathrm{\Phi }_{\mathrm{BP}}^{(0)}=\pi n(1\mathrm{cos}\theta )`$. We wish to know what we observe if we carry out the same measurement for a spin which has been weakly coupled to a bath of oscillators throughout the experiment.

The Hamiltonian we consider contains the spin-half in the above time-dependent magnetic field, $`𝑩(t)`$, which is also coupled to a bath of harmonic oscillators with frequencies $`\{\mathrm{\Omega }_j\}`$. Writing it in terms of creation, $`\widehat{b}^{}`$, and annihilation, $`\widehat{b}`$, operators for the oscillators,

$`(t)`$ $`=`$ $`{\displaystyle \frac{g}{2}}𝑩(t)\widehat{𝝈}+{\displaystyle \underset{j,\alpha }{}}\mathrm{\Omega }_j\left(\widehat{b}_{j,\alpha }^{}\widehat{b}_{j,\alpha }+\frac{1}{2}\right)`$ (3)
$`{\displaystyle \frac{g}{2}}{\displaystyle \underset{j,\alpha }{}}{\displaystyle \frac{C_\alpha }{(2m\mathrm{\Omega }_j)^{1/2}}}\left(\widehat{b}_{j,\alpha }^{}+\widehat{b}_{j,\alpha }\right)\widehat{\sigma }_\alpha ,`$
where $`j`$ is summed over all oscillators and $`\alpha `$ is summed over the $`(x,y,z)`$ components of the oscillator. The number of oscillators with frequency $`\mathrm{\Omega }`$ to $`\mathrm{\Omega }+\mathrm{d}\mathrm{\Omega }`$ is $`p(\mathrm{\Omega })\mathrm{d}\mathrm{\Omega }`$. The spectral density Leggett87  is given by $`J(\mathrm{\Omega })=_j\pi (gC)^2\left(2m\mathrm{\Omega }_j\right)^1\delta (\mathrm{\Omega }\mathrm{\Omega }_j)=\pi (gC)^2p(\mathrm{\Omega })\left(2m\mathrm{\Omega }\right)^1`$. Here we restrict ourselves to $`z`$-axis spin-environment coupling with $`C_\alpha =C\delta _{\alpha z}`$ endnote-isotropic . Then for $`𝑩𝝈=B\sigma _z`$, the exact ground state of the universeShnirman02  is simply $`|_j|0_{}^j`$ where oscillator $`j`$ is in the ground state, $`|0_{}^j`$, of the harmonic potential centred at $`(0,0,\frac{1}{2}gC)`$. We consider an Ohmic bath of oscillators with $`J(\mathrm{\Omega })=\frac{\pi }{2}\stackrel{~}{C}^2\mathrm{\Omega }\mathrm{exp}\left[\mathrm{\Omega }/\mathrm{\Omega }_\mathrm{m}\right]`$, and work in the limit of small dimensionless coupling $`\stackrel{~}{C}gC(\alpha /m)^{1/2}1`$.

The time-dependence in (3) makes the problem unpleasant, however we remove this by going to the primed-basis which rotates with the $`𝑩`$-field. In this non-inertial basis the spin experiences an effective field $`\left(𝑩_0+g^1𝝎\right)`$. For our problem the effective field is $`𝑩_+`$ for $`0<t<t_\mathrm{p}`$ (shown in Fig. 1), and $`𝑩_{}`$ for $`t_\mathrm{p}<t<2t_\mathrm{p}`$, where $`𝑩_\pm =\left(𝑩_0\pm g^1\omega \widehat{𝒛}\right)`$. Having removed the time-dependent, we calculate the evolution of the system in the frame which has its z-axis parallel to the field (either $`𝑩_+`$ or $`𝑩_{}`$) these frames we call the plus\- and minus-basis respectively (the former is shown in Fig. 1). Finally we rotate back to the lab-frame to evaluate observables.

Before we give a detailled explanation of how we calculate the spin’s evolution in the presence of the environment, we giving our results. The anisotropic nature of the coupling results in $`P(\xi )`$ containing $`𝒪[\stackrel{~}{C}^2]`$-terms which go like $`\mathrm{exp}[\pm \mathrm{i}gBt_\mathrm{p}]`$. To simplify the resulting expressions we average $`t_\mathrm{p}`$ over a range $`(gB)^1`$ to remove these terms, then

$`P(\xi )`$ $`=`$ $`\frac{1}{2}[1+\mathrm{e}^{2t_\mathrm{p}/T_2}\mathrm{cos}(\xi 4\mathrm{\Phi }_{\mathrm{BP}}\kappa _1)`$ (4)
$`+\left|\kappa _2\right|\mathrm{e}^{2t_\mathrm{p}/T_2}\mathrm{cos}\left(\xi +4\mathrm{\Phi }_{\mathrm{BP}}\mathrm{arg}[\kappa _2]\right)`$
$`+\left|\kappa _3\right|\left(2\mathrm{e}^{2t_\mathrm{p}/T_2}\mathrm{e}^{4t_\mathrm{p}/T_2}\right)\mathrm{cos}\left(\xi \mathrm{arg}\left[\kappa _3\right]\right)`$
$`\kappa _4\mathrm{cos}\xi ],`$
where $`\mathrm{\Phi }_{\mathrm{BP}}=\mathrm{\Phi }_{\mathrm{BP}}^{(0)}+\delta \mathrm{\Phi }_{\mathrm{BP}}`$. For compactness we have dropped an uninteresting real $`𝒪[\stackrel{~}{C}^2t_\mathrm{p}^0]`$ term from the first exponent while retaining such terms elsewhere. The $`\kappa `$s (which were schematically refered to as $`\mathrm{\Phi }_{\mathrm{shift}}`$ above) are $`𝒪[\stackrel{~}{C}^2]`$ and so are comparable to $`\delta \mathrm{\Phi }_{\mathrm{BP}}`$; however they are independent of the time-dependence of $`(t)`$, and hence independent of the winding number, $`n`$. We find,

$`T_2^1`$ $`=`$ $`(2T_1)^1=\frac{\pi }{8}\stackrel{~}{C}^2\mathrm{\Omega }_\mathrm{m}\gamma \mathrm{e}^\gamma \mathrm{sin}^2\theta `$ (5)
$`\delta \mathrm{\Phi }_{\mathrm{BP}}`$ $`=`$ $`\frac{\pi }{8}n\stackrel{~}{C}^2\left[f^{}(\gamma )2\gamma ^1f(\gamma )\right]\mathrm{sin}^2\theta \mathrm{cos}\theta ,`$ (6)
where $`\gamma =gB/\mathrm{\Omega }_\mathrm{m}`$. The function $`f(x)=x\mathrm{e}^x\mathrm{Ei}(x)+x\mathrm{e}^x\mathrm{Ei}(x)`$ where we define $`\mathrm{Ei}(x)`$ as the principal-value of the Exponential integral, $`_x^{\mathrm{}}dt\mathrm{e}^t/t`$, and $`f^{}(x)\mathrm{d}f(x)/\mathrm{d}x`$. Eq. (6) is simply the $`\omega `$-dependent term in the Lamb shift of the energy when in the rotating frame. This generates a term of $`𝒪[nt_\mathrm{p}^0]`$ in the phase which in the laboratory frame is a contribution to the BP.

The $`n`$-independent factors are

$`\kappa _1`$ $`=`$ $`\frac{1}{4}\stackrel{~}{C}^2\pi \mathrm{e}^\gamma \mathrm{sin}\theta \left[\mathrm{cos}\theta +\frac{1}{4}\mathrm{sin}\theta \right]`$
$`\kappa _2`$ $`=`$ $`\frac{1}{16}\stackrel{~}{C}^2\left(\gamma ^1f(\gamma )+\mathrm{i}\pi \mathrm{e}^\gamma \right)\mathrm{sin}^2\theta `$
$`\kappa _3`$ $`=`$ $`\frac{1}{4}\stackrel{~}{C}^2\left[\left(\gamma ^1f(\gamma )+\mathrm{i}\pi \mathrm{e}^\gamma \right)\mathrm{sin}\theta \mathrm{cos}\theta 2\gamma ^1\mathrm{sin}\theta \right]`$
$`\kappa _4`$ $`=`$ $`\frac{1}{2}\stackrel{~}{C}^2\left[\mathrm{e}^\gamma \mathrm{Ei}(\gamma )\mathrm{sin}\theta \mathrm{cos}\theta \gamma ^1\mathrm{sin}\theta \right].`$ (7)
Now we check that the BP is not simply given by the solid-angle enclosed by the spin-resonances. If this were the case then the BP for this experiment would be $`\mathrm{\Phi }_{\mathrm{BP}}^{(0)}\frac{\pi }{4}n\stackrel{~}{C}^2\gamma ^1f(\gamma )\mathrm{sin}^2\theta \mathrm{cos}\theta `$, the correct result deviates from this prediction by

$`\delta \mathrm{\Phi }_{\mathrm{BP}}^{}`$ $`=`$ $`\frac{\pi }{8}n\stackrel{~}{C}^2f^{}(\gamma )\mathrm{sin}^2\theta \mathrm{cos}\theta ,`$ (8)
for most $`\gamma `$ this deviation is significant.

We now discuss the method we use to obtain these results. The Hamiltonian in the primed-basis is time-independent and is given by

$`^{}(𝑩_\pm )`$ $`=`$ $`{\displaystyle \frac{g}{2}}𝑩_\pm \widehat{𝝈}+{\displaystyle \underset{j}{}}\mathrm{\Omega }_j\left(\widehat{b}_j^{}\widehat{b}_j+\frac{1}{2}\right)`$ (9)
$`{\displaystyle \frac{g}{2}}{\displaystyle \underset{j}{}}{\displaystyle \frac{C}{(2m\mathrm{\Omega }_j)^{1/2}}}\left(\widehat{b}_j^{}+\widehat{b}_j\right)\widehat{\sigma }_z.`$
If we write the spin’s initial density matrix as $`\rho _0`$ and the oscillators initial density matrix as $`\rho _0^{\mathrm{osc}}`$, then we are interested in the spin density matrix at time $`t`$, after we have traced over the oscillator states, $`\rho _t=\text{tr}_{\mathrm{osc}}\widehat{U}_t\left(\rho _0\rho _0^{\mathrm{osc}}\right)\widehat{U}_t^{}`$, where $`\widehat{U}_t`$ is the evolution operator. We find it helpful to write the spin density matrix as a vector $`𝝆`$ whose elements are $`(\rho _{11},\rho _{12},\rho _{21},\rho _{22})`$. Then the spin evolution equation (after the oscillators have been traced over) can be written as $`𝝆_t=𝐊(t)𝝆_0`$, where this defines $`𝐊(t)`$ as a four-by-four matrix which gives the time evolution of the elements of the spin’s density matrix. The initial state of the oscillators enters in the functional form of the elements of $`𝐊(t)`$. For the experiment described above eq. (1) we need to calculate $`𝝆_{2t_\mathrm{p}}=𝐊^{\mathrm{flip}}𝐊(𝑩_+,t_\mathrm{p})𝐊^{\mathrm{flip}}𝐊(𝑩_{},t_\mathrm{p})𝝆_0`$. The spin-flip is assumed to be fast enough to leave all the oscillators unchanged while flipping the spin, then in the primed-basis $`𝐊^{\mathrm{flip}}`$ simply has “1”s on the off-diagonal and “0”s elsewhere. This leaves the calculation of the propagation matrix $`𝐊(𝑩_+,t_\mathrm{p})`$, we can find $`𝐊(𝑩_{},t_\mathrm{p})`$ by reversing the sign of $`\omega `$ throughout. For weak coupling to the bath it is natural to work in the plus-basis (see Fig. 1), which has its $`z`$-axis parallel to $`𝑩_+`$, in this basis $`𝐊(𝑩_+,t_\mathrm{p})`$ becomes diagonal if $`\stackrel{~}{C}0`$. Finally the coupling between spin and oscillators in the plus-basis at time $`t`$ is $`𝐂_+(t)=𝐂^{}(t)𝐑_+`$, where $`𝐑_+`$ is the SO(3) rotation from the primed-basis to the plus-basis.

Now we use the real-time transport method Schoeller94  to write the following differential equation for $`𝐊^+(t)`$,

$`_t𝐊^+(t)=\mathrm{i}𝐄^+𝐊^+(t)+{\displaystyle _0^t}d\tau 𝚺^+(\tau )𝐊^+(t\tau ),`$ (10)
where all bold symbols are $`4\times 4`$ matrices. The matrix $`𝐄^+`$ gives the evolution of the propagation matrix when there is no coupling to the bath. Because we are in the plus-basis it is diagonal with $`E_{11}^+=E_{44}^+=0`$ and $`E_{22}^+=E_{33}^+=gB^+g\left|𝑩+g^1𝝎\right|`$. The matrix $`𝚺^+(\tau )`$ is the contribution of all irreducible diagrams with one or more interactions with the bath of oscillators. Equation (10) is exact, however to proceed we treat this equation to first order in $`\stackrel{~}{C}^2`$. Thus in the integral on the left hand side of (10) we treat $`𝚺^+(\tau )`$ to first order in $`\stackrel{~}{C}^2`$ and $`𝐊^+(t\tau )`$ to zeroth order. So we can write $`𝐊^+(t\tau )𝐊^+(t)𝐊_0^+(\tau )`$ where the corrections to the approximations are $`𝒪[\stackrel{~}{C}^2]`$ and so can be ignored. Now $`𝚺^+`$, which is evaluated below, is dominated by small-$`\tau `$ so we take the upper limit on the integral to infinity. The error we make in doing so is $`𝒪\left[(\mathrm{\Omega }_\mathrm{m}t_\mathrm{p})^1\right]`$ which we neglect. This systematic approximation results in the interaction becoming local in time. Then we get $`_t𝐊^+(t)=\left[\mathrm{i}𝐄^++𝐗^+\right]𝐊^+(t)`$ where $`𝐗^+=_0^{\mathrm{}}d\tau 𝚺^+(\tau )𝐊_0^+(\tau )`$, and diagonalise the matrix $`\left(\mathrm{i}𝐄^++𝐗^+\right)`$ to find $`𝐊^+(t)`$.

Now we briefly discuss the evaluation of $`𝚺^+`$ to lowest order in $`\stackrel{~}{C}^2`$. At this order we need only consider irreducible diagrams with a single interaction with the oscillators. When the oscillators are traced out they leave an interaction between the spin at time $`t`$ and time $`t\tau `$. The resulting first-order irreducible diagrams are shown in Fig. 2. The contribution to $`𝚺^+`$ of the diagram with an interaction via $`\sigma _\alpha `$ at time $`t`$ and another via $`\sigma _\alpha ^{}`$ at time $`t^{}`$, after we have summed over the Ohmic bath, is

$`\chi \frac{g^2\alpha }{8m}\left[𝐂_+^\mathrm{T}(t)𝐂_+(t^{})\right]_{\alpha \alpha ^{}}`$
$`\times \left\{{\displaystyle \genfrac{}{}{0pt}{}{\left[\sigma _\alpha \right]_{\mu _i\mu }}{\left[\sigma _\alpha \right]_{\stackrel{~}{\mu }\stackrel{~}{\mu }_i}}}\right\}\left\{{\displaystyle \genfrac{}{}{0pt}{}{\left[\sigma _\alpha ^{}\right]_{\mu \mu _j}}{\left[\sigma _\alpha ^{}\right]_{\stackrel{~}{\mu }_j\stackrel{~}{\mu }}}}\right\}{\displaystyle \frac{\mathrm{\Omega }_\mathrm{m}^2\mathrm{e}^{\mathrm{i}gB^+(\mu \stackrel{~}{\mu })\tau /2}}{(1+\mathrm{i}\kappa \mathrm{\Omega }_\mathrm{m}\tau )^2}},`$ (11)
where $`\kappa =\pm 1`$, $`\chi =\pm 1`$, and other variables are shown in Fig. 2. The upper (lower) term in $`\left\{\mathrm{}\right\}`$ is applicable if the relevant vertex is on R (A). $`\kappa `$ is $`+1`$ ($`1`$) when the $`\alpha ^{}`$ vertex is on R (A). $`\chi `$ is $`+1`$ ($`1`$) if the interaction is R–R or A–A (R–A or A–R).

In conclusion, the BP can be observed in a non-isolated system, if the coupling to the environment is weak enough that $`gBT_2^1`$. The adiabatic phase is $`\mathrm{\Phi }_{\mathrm{BP}}^{(0)}+\delta \mathrm{\Phi }_{\mathrm{BP}}+\mathrm{\Phi }_{\mathrm{shift}}`$, but $`\mathrm{\Phi }_{\mathrm{shift}}`$ is not considered a BP because it does not vanish when $`n=0`$. So the BP differs from that of a isolated spin by $`\delta \mathrm{\Phi }_{\mathrm{BP}}`$, given in Eq. (6). The proportionality of $`\delta \mathrm{\Phi }_{\mathrm{BP}}`$ to $`n`$ hints that it has some geometric character, however it is a function of the environment’s spectrum and thus the total BP is not a simple geometric quantity.

We are extremely grateful to A. Shnirman for useful discussions, and we thank R. Fazio, F. Wilhelm and Y. Aharonov for enlightening comments. This work was commenced while RW was working at the Weizmann Institute and was supported by the U.S.-Israel Binational Science Foundation (BSF), by the Minerva Foundation, by the Israel Science Foundation, and by the German-Israel Foundation (GIF).


text:
# Evaluating holonomic quantum computation: beyond adiabatic limitation

## Abstract

The proposal of the optical scheme for holonomic quantum computation is evaluated based on dynamical resolution to the system beyond adiabatic limitation. The time-dependent Schrödinger equation is exactly solved by virtue of the cranking representation and gauge transformation approach. Besides providing rigorous confirmation to holonomies of the geometrical prediction that holds for the ideally adiabatic situation, the dynamical resolution enables one to evaluate elaborately the amplitude of the nonadiabatic deviation, so that the errors induced to the quantum computation can be explicitly estimated.

The recently proposed holonomic approach to quantum computation \-  surely predicts a striking contribution to the application of quantum physics. Transcending the traditional dynamical means for quantum computation, the holonomic approach realizes quantum information processing by endowing the quantum code with a non-trivial global topology (a gauge field potential) and the associated holonomies then allow for the universal quantum computing. Specifically, in the scheme of holonomic quantum computation, information is encoded in a degenerate eigenspace of the governing Hamiltonian and the holonomies (abelian as well as non-abelian) - are acquired by driving the system to undergo appropriate loops in the parameter space adiabatically. Besides suggesting an intriguing connection between the gauge fields and the information processing, such a geometrical means for quantum manipulation is believed to have built-in fault-tolerant features - due to its inherent stability against local perturbations. Considerable attention has been addressed to this topic recently and the all-geometrical implementation for universal quantum gates has been proposed by optical schemes, based on laser manipulation of ions confined in a Paul trap  or neutral atoms in an optical resonator .

The existent exploration for holonomic quantum computation is based on the analysis by pure geometrical fashion. It is true that in the adiabatic limit the holonomy associated with the evolving loop is determined by the path traced by the time-evolution ray and the curvature of the ray space. This involves the abelian holonomy (the Berry phase) and the non-abelian one merely known as adiabatic connection. Nevertheless, as a whole physical problem, as the dynamics of a system generates a time-dependent physical state, a specified geometrical object (the ray) is generated as well. In such a sense, dynamics determines the holonomy through determining the ray itself and its path. Moreover, in view that the realistic evolution of a physical system could not be ideally adiabatic and the nonadiabaticity shall alter the time-evolution of the ray and thus inevitably induce deviation from the adiabatic consequences. The evaluation of such deviation and the resulting errors in quantum computation is definitely a dynamical problem that goes beyond the geometrical exploration.

In this paper we employ a tractable model of the optical scheme to exploit this subject. For the appropriately chosen loops of the Hamiltonian in the parameter space, the time-dependent Schrödinger equation is exactly solved by virtue of the cranking representation and gauge transformation approach. The derived dynamical evolution of the system recovers the holonomic transformation provided by geometrical consequences, including the simple abelian phase factor and the general non-abelian operation. Thus our results provide further confirmation of the geometrical prediction, and besides, the errors caused by nonadiabatic effects for the holonomic quantum gate operation can be estimated explicitly.

For the proposed optical scheme of holonomy quantum computation , the basic idea relies on the adiabatic passage via the dark states since the dynamical evolution restricted to such a space is completely trivial. The system encoding the qubit is realized by a four-level $`\mathrm{\Lambda }`$-type trapped ion (or a similar cavity atom). The three ground levels $`|g_i`$ ($`i=1,2,3`$) are highly degenerate and each couples to the excited state $`|e`$ in a tunable way. The states $`|g_1`$ and $`|g_2`$ stand for the computational bases $`|0`$ and $`|1`$, respectively, and $`|g_3`$ is an ancillary level required for implementation of gate operations. Such a system admits two dark states that have no contribution from the excited state. Through changing the Rabi frequencies and driving the dark states to undergo appropriate cyclic evolutions in an adiabatic fashion, the universal single-bit gate operations $`e^{i\varphi |11|}`$ and $`e^{i\varphi \sigma _y}`$ can be generated due to the global geometry of the bundle of the eigenspace of the dark states.

To evaluate the gate operation $`e^{i\varphi |11|}`$ from a dynamical viewpoint, let us explore the state evolution generated by the periodic Hamiltonian

$$H(t)=\mathrm{\Omega }\mathrm{sin}\theta (\sigma _{2e}+\sigma _{e2})+\mathrm{\Omega }\mathrm{cos}\theta (\sigma _{3e}e^{i\phi }+\sigma _{e3}e^{i\phi }),$$
(1)
where $`\theta `$ is a fixed parameter and $`\phi `$ is assumed to rotate at a constant frequency $`\gamma `$ for convenience. The equation of motion for the system is

$$i\frac{}{t}|\mathrm{\Psi }(t)=H(t)|\mathrm{\Psi }(t).$$
(2)
It is known that in the adiabatic limit, the geometrical exploration shows that the dark state of the system, $`|D(t)=\mathrm{cos}\theta |g_2\mathrm{sin}\theta e^{i\gamma t}|g_3`$, shall acquire a net Berry phase  after a period $`T=2\pi /\gamma `$: $`|D(T)=e^{i\varphi }|D(0)`$ with $`\varphi =4\pi \mathrm{sin}^2\theta `$. To give a dynamical resolution to the system beyond adiabatic limitation, we use the cranking representation, that is, the Hamiltonian (1) can be regarded as a cranked one

$$H(t)=e^{i\gamma t\sigma _{33}}H_0e^{i\gamma t\sigma _{33}},$$
(3)
where

$$H_0=\mathrm{\Omega }\mathrm{sin}\theta (\sigma _{2e}+\sigma _{e2})+\mathrm{\Omega }\mathrm{cos}\theta (\sigma _{3e}+\sigma _{e3}),$$
(4)
and the unitary transformation $`e^{i\gamma t\sigma _{33}}`$ can be viewed as an element of the SU(3) group. Consequently, the dynamical invariant of the system can be shown as

$$I(t)=e^{i\gamma t\sigma _{33}}(H_0+\gamma \sigma _{33})e^{i\gamma t\sigma _{33}}=H(t)+\gamma \sigma _{33},$$
(5)
which satisfies

$$\frac{dI(t)}{dt}=\frac{I(t)}{t}i[I(t),H(t)]=0.$$
(6)
The second term of (5) accounts for an extra gauge potential since $`H(t)`$ depends on time explicitly. Now the recurrent basis $`|\psi (t)`$ of the system \[it differs from the basic solution $`|\mathrm{\Psi }(t)`$ of the Schrödinger equation (2) only by a phase factor\] can be obtained by solving the instantaneous eigensolutions of $`I(t)`$. It turns out, one only needs to solve the characteristic equation

$$x^3(\gamma /\mathrm{\Omega })x^2x+(\gamma /\mathrm{\Omega })\mathrm{sin}^2\theta =0$$
(7)
and the eigenvalues of $`I(t)`$ are given by $`E_i(\frac{\gamma }{\mathrm{\Omega }})=\mathrm{\Omega }x_i(\frac{\gamma }{\mathrm{\Omega }})`$ ($`i=1,0,1`$). It is straightforward to show that the recurrent basis $`|\psi _0(t)`$ represented by the middle number $`E_0`$ approaches to the dark state $`|D(t)`$ in the adiabatic limit $`\gamma /\mathrm{\Omega }0`$. Now the leakage error induced by the nonadiabatic effect can be conveniently estimated by the overlap \[see Fig. 1(a)\]

$$\eta (\theta ,\frac{\gamma }{\mathrm{\Omega }})=|\psi _0(0)|D(0)|^2=|\psi _0(T)|D(T)|^2.$$
(8)
Besides the leakage, the nonadiabatic evolution shall result in deviation to the desired phase factor. It follows, instead of the net Berry phase, the cyclic evolution here induces a total phase (the so-called Lewis-Riesenfeld phase)

$$\mathrm{\Phi }=_0^T\psi _0(t)|i\frac{}{t}H(t)|\psi _0(t)𝑑t=E_0\frac{2\pi }{\gamma }.$$
(9)
The detailed depiction of the deviation for the phase factor is shown in Fig. 1(b). Noting that in the adiabatic limit, the total phase

$$\mathrm{\Phi }=\underset{\gamma /\mathrm{\Omega }0}{lim}2\pi \frac{x_0(\frac{\gamma }{\mathrm{\Omega }})}{\gamma /\mathrm{\Omega }}=4\pi \mathrm{sin}^2\theta ,$$
(10)
the geometrical consequence is thus recovered.

The validity of the above evaluation is based on a presumption that the initial state $`|D(0)=\mathrm{cos}\theta |g_2\mathrm{sin}\theta |g_3`$ can be generated from the computational basis $`|g_2`$ and so the inverse process. Explicitly, such processes can be accomplished by the driven Hamiltonian (1) through changing the parameter $`\theta `$ adiabatically. Conventionally, the nonadiabatic effect here shall lead to an additional error for the quantum computation. However, such an error can be in principle avoided through appending a matching interaction to compensate the gauge potential term induced to the system. Specifically, one can use the following Hamiltonian (setting $`\phi =0`$)

$$H_{tot}(t)=H(t)+H_{ad}(t),H_{ad}(t)=i\dot{\theta }(t)(\sigma _{23}\sigma _{32}).$$
(11)
It follows that the dynamical invariant of the system $`H_{tot}(t)`$ now has a form $`I(t)=H(t)`$, thus the above state transformation can be processed exactly. Physically, the interaction $`H_{ad}(t)`$ can be realized by a microwave coupling to the two degenerate levels $`|g_2`$ and $`|g_3`$, with its intensity accurately controlled through a derivative feedback process.

Now we investigate the gate operation $`e^{i\varphi \sigma _y}`$ achieved by the holonomic means. The corresponding evolution is generated by the Hamiltonian

$`H(t)`$ $`=`$ $`\mathrm{\Omega }\mathrm{sin}\theta \mathrm{cos}\phi (\sigma _{1e}+\sigma _{e1})`$ (12)
$`+`$ $`\mathrm{\Omega }\mathrm{sin}\theta \mathrm{sin}\phi (\sigma _{2e}+\sigma _{e2})+\mathrm{\Omega }\mathrm{cos}\theta (\sigma _{3e}+\sigma _{e3}),`$ (13)
where the parameter $`\phi =\gamma t`$. As is known, the adiabatic cyclic evolution of the Hamiltonian generates a non-abelian holonomy due to its degeneracy structure of the dark states. It can be easily worked out, from the formula of Ref. , that the holonomic transformation

$$u_C=e^{i2\pi \mathrm{cos}\theta D_y},$$
(14)
where $`D_y=i(|D_2D_1||D_1D_2|)`$, and the two dark states, $`|D_1=\mathrm{cos}\theta |g_1\mathrm{sin}\theta |g_3`$ and $`|D_2=|g_2`$, span the degenerate space of the starting (ending) Hamiltonian. Note that the Hamiltonian (13) possesses an su(4) Lie algebraic structure and dynamical resolution to the system is usually very complicated. Surprisingly, as we shall show in the following, this system can be exactly solved by the gauge transformation approach , and its dynamical evolution analytically manifested thus leads to a complete understanding of the adiabatic and nonadiabatic properties for the time-dependent Hamiltonian system.

Similar to the cranking method used above, we introduce the unitary gauge transformation

$$U_g(t)=e^{\gamma t(\sigma _{12}\sigma _{21})}$$
(15)
to the equation of motion for the system, from which a covariant Schrödinger equation is stemmed

$`|\mathrm{\Psi }_g(t)`$ $`=`$ $`U_g^1(t)|\mathrm{\Psi }(t),`$ (16)
$`i{\displaystyle \frac{}{t}}|\mathrm{\Psi }_g(t)`$ $`=`$ $`H_g|\mathrm{\Psi }_g(t),`$ (17)
with the gauged Hamiltonian

$`H_g`$ $`=`$ $`U_g^1HU_giU_g^1{\displaystyle \frac{U_g}{t}}`$ (18)
$`=`$ $`\mathrm{\Omega }\mathrm{sin}\theta (\sigma _{1e}+\sigma _{e1})`$ (19)
$`+`$ $`\mathrm{\Omega }\mathrm{cos}\theta (\sigma _{3e}+\sigma _{e3})+i\gamma (\sigma _{12}\sigma _{21}).`$ (20)
In view that the above Hamiltonian is time independent, the basic solutions $`|\mathrm{\Psi }_g^n(t)`$ to the covariant equation (17) can be easily obtained and the corresponding eigenvalues are as follows

$`E_{1,2}`$ $`=`$ $`\pm {\displaystyle \frac{\sqrt{2}}{2}}\overline{\mathrm{\Omega }}\left[1\sqrt{14({\displaystyle \frac{\gamma }{\mathrm{\Omega }}})^2\mathrm{cos}^2\overline{\theta }}\right]^{1/2},`$ (21)
$`E_{3,4}`$ $`=`$ $`\pm {\displaystyle \frac{\sqrt{2}}{2}}\overline{\mathrm{\Omega }}\left[1+\sqrt{14({\displaystyle \frac{\gamma }{\mathrm{\Omega }}})^2\mathrm{cos}^2\overline{\theta }}\right]^{1/2},`$ (22)
where

$$\overline{\mathrm{\Omega }}=\mathrm{\Omega }\sqrt{1+(\gamma /\mathrm{\Omega })^2},\mathrm{cos}\overline{\theta }=\frac{\mathrm{cos}\theta }{1+(\gamma /\mathrm{\Omega })^2}.$$
(23)
The dynamical basis of the system (13) can be directly obtained as $`|\mathrm{\Psi }_n(t)=U_g|\mathrm{\Psi }_g^n(t)`$, from which one can see that $`E_n`$ has the natural implication related to the total phase. Now the time evolution operator generated by the Hamiltonian (13) can be given

$`U_C(T)`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{4}{}}}|\mathrm{\Psi }_n(T)\mathrm{\Psi }_n(0)|`$ (24)
$`=`$ $`{\displaystyle \underset{n=1}{\overset{4}{}}}e^{iE_n\frac{2\pi }{\gamma }}|\mathrm{\Psi }_n(0)\mathrm{\Psi }_n(0)|.`$ (25)
Considering the asymptotic behavior of the evolution in the adiabatic limit, it follows that $`lim_{\gamma /\mathrm{\Omega }0}\frac{E_{1,2}}{\gamma }=\pm \mathrm{cos}\theta `$, and the phase-equipped dynamical bases $`|\mathrm{\Psi }_1(t)`$ and $`|\mathrm{\Psi }_2(t)`$ have the form

$`|\mathrm{\Psi }_1(t)`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{2}}e^{i\gamma t\mathrm{cos}\theta }[(\mathrm{cos}\theta \mathrm{cos}\gamma t+i\mathrm{sin}\gamma t)|g_1`$ (27)
$`+(\mathrm{cos}\theta \mathrm{sin}\gamma ti\mathrm{cos}\gamma t)|g_2\mathrm{sin}\theta |g_3],`$
$`|\mathrm{\Psi }_2(t)`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{2}}e^{i\gamma t\mathrm{cos}\theta }[(\mathrm{cos}\theta \mathrm{cos}\gamma ti\mathrm{sin}\gamma t)|g_1`$ (29)
$`+(\mathrm{cos}\theta \mathrm{sin}\gamma ti\mathrm{cos}\gamma t)|g_2\mathrm{sin}\theta |g_3].`$
One can verify that they are the instantaneous eigenstates of the Hamiltonian (13) with a two-degeneracy eigenvalue $`0`$, and the equipped phases are just the Berry phases accordingly. Thus the cyclic evolution restricted to the space spanned by these two states is purely geometrical and can be denoted as

$$u(T)=e^{i2\pi \mathrm{cos}\theta }|\mathrm{\Psi }_1(0)\mathrm{\Psi }_1(0)|+e^{i2\pi \mathrm{cos}\theta }|\mathrm{\Psi }_2(0)\mathrm{\Psi }_2(0)|$$
(30)
with $`|\mathrm{\Psi }_1(0)=(|D_1i|D_2)/\sqrt{2}`$ and $`|\mathrm{\Psi }_2(0)=(|D_1+i|D_2)/\sqrt{2}`$. It can be easily recognized that the operator (30) is just the non-abelian holonomy (14), thus the geometrical nature is verified again.

The above dynamical resolution to the system is important. Besides offering a vivid verification to the remarkable formula of non-abelian holonomy , which holds for the ideally adiabatic situation, it enables one to evaluate elaborately the amplitude of the nonadiabaticity deviation and the resulting errors to the holonomic gate operation $`e^{i\varphi \sigma _y}`$. In detail, the population transfer from the initial state $`|\mathrm{\Psi }(0)=|g_2`$ is pictured in Fig. 2. The leakage out of the computational space can be estimated by the projection (see also Fig. 2)

$$\eta (\theta ,\frac{\gamma }{\mathrm{\Omega }})=\underset{i=1}{\overset{2}{}}|D_i|U(T)|\mathrm{\Psi }(0)|^2.$$
(31)
Similar to the former case, to transform the computational basis $`|g_1`$ into the dark state $`|D_1=\mathrm{cos}\theta |g_1\mathrm{sin}\theta |g_3`$ and to invert the process successfully, one needs to use the Hamiltonian (13) (with $`\theta `$ tunable and $`\phi =0`$) along with a matching interaction

$$H_{tot}(t)=H(t)+H_{ad}(t),H_{ad}(t)=i\dot{\theta }(t)(\sigma _{13}\sigma _{31}).$$
(32)
Up to now, we have investigated the single-qubit holonomic operations of the optical scheme and evaluated the nonadiabaticity-causing errors for quantum computation. It deserves to point out that, the proposal  of geometrical implementation for the controlled two-qubit phase shift gate, $`e^{i\varphi |1111|}`$, which is sufficient for the universal quantum computation along with the two single-qubit gates, can be explored in a similar way. The scheme is realized by two-color laser manipulation on the trapped ions . Briefly, the transition from the ground states $`|g_2`$ and $`|g_3`$ to the excited state $`|e`$ is driven by two different bi-chromatic laser beams with their amplitudes and phases of the Rabi frequencies controllable. The frequencies of the laser fields are tuned so that the two-photon process, exciting pair ions, is resonant and the single-photon excitation is off-resonant. Hence the system can be described by an effective Hamiltonian in the notation of Ref.

$`H_{eff}`$ $``$ $`|\mathrm{\Omega }_1|^2(e^{i2\phi _1}|eeg_2g_2|+h.c.)`$ (33)
$`+`$ $`|\mathrm{\Omega }_2|^2(e^{i2\phi _2}|eeg_3g_3|+h.c.),`$ (34)
where the relative intensity of the Rabi frequencies $`\mathrm{tan}\theta =|\mathrm{\Omega }_1|^2/|\mathrm{\Omega }_2|^2`$ and the phase difference $`\phi /2=\phi _1\phi _2`$ are tunable. One can see that the bases $`|g_1g_1`$, $`|g_1g_2`$, and $`|g_2g_1`$ are decoupled from the evolution, and the component $`|g_2g_2`$ serving as the code $`|11`$ evolved in an enclosed space spanned by {$`|g_2g_2,|g_3g_3,|ee`$}. Introducing the su(3) generators explicitly

$`A_{e2}`$ $`=`$ $`e^{i2\phi _1}|eeg_2g_2|,A_{e3}=e^{i2\phi _1}|eeg_3g_3|,`$ (35)
$`A_{23}`$ $`=`$ $`|g_2g_2g_3g_3|,A_{\mu \nu }^{}=A_{\nu \mu },`$ (36)
the Hamiltonian (34) can be rewritten as

$$H_{eff}=g\mathrm{sin}\theta (A_{2e}+A_{e2})+g\mathrm{cos}\theta (A_{3e}e^{i\phi }+A_{e3}e^{i\phi }).$$
(37)
Obviously this Hamiltonian possesses an su(3) algebraic structure isomorphic to that of system (1), thus all the discussions therein also hold for the present system.

It should be noted that, the effective Hamiltonian (34), respecting a second-order process of the interaction, is quite a rough description of the model. Specifically, it ignores the same second-order process induced by virtual photons excitation in the self-transitions of the states: $`|g_{2(3)}|g_{2(3)}`$ and $`|e|e`$. It can be anticipated that such self transitions shall dress the energy levels of the ions and lift the degeneracy of the ground states, which in turn affects the desired gate operation. Detailed exploration of this point shall be presented in a future report.

This work was supported in part by the Postdoctoral Science Foundation, the special funds for Major State Basic Research Project No. G001CB3095, the National Natural Science Foundation No. 10175029 of China, and the Research Grants Council of the Hong Kong Government.

Caption of Fig. 1:

Deviation induced by nonadiabaticity for abelian holonomy. (a) The overlap $`\eta `$ for the parameters $`\theta (0,\pi /2)`$ and $`\gamma /\mathrm{\Omega }[0,1]`$; (b) The total phase $`\mathrm{\Phi }`$ for $`\theta [0,\pi ]`$ and $`\gamma /\mathrm{\Omega }[0,1]`$.

Caption of Fig. 2:

Deviation induced by nonadiabaticity for non-abelian holonomy. The initial state is prepared in $`|D_2`$. The two solid curves show the results for the population of the target state on $`|D_1`$ and $`|D_2`$, as a function of $`1\mathrm{cos}\theta `$, respectively. The dashed curve depicts the total population $`\eta `$ on the computational space. Figures (a), (b), (c) and (d) correspond to $`\gamma /\mathrm{\Omega }=0.01,0.2,0.5`$ and $`0.8`$, respectively.