text:
# Bubbling and Regularity Issues in Geometric Non-linear Analysis

## 1. Energy quantization phenomena for harmonic maps and Yang-Mills Fields

### 1.1. The archetype of energy quantization: the $`𝜺`$-regularity

Let $`B^m`$ be the flat $`m`$-dimensional ball and $`N^n`$ be a compact, without boundary Riemannian manifold. By the Nash embedding Theorem, we may assume $`N\mathrm{lR}^k`$. Let $`W^{1,2}(B^m,N^n)`$ be the maps in $`W^{1,2}(B^m,\mathrm{lR}^k`$ that take values almost everywhere in $`N^n`$. A map $`uW^{1,2}(B^m,N^n)`$ is called stationary harmonic if it is a critical point for the Dirichlet energy

$$E(u)=_{B^m}|u|^2𝑑x$$

for both perturbations in the target (of the form $`\pi _N(u+t\varphi )`$ for any $`\varphi `$ in $`C_0^{\mathrm{}}(B^m,\mathrm{lR}^k)`$ where $`\pi _N(y)`$ is the nearest neighbor of $`y`$ in $`N`$) and in the domain (of the form $`u(id+tX)`$ for any vector field $`X`$ in $`C_0^{\mathrm{}}(B^m,\mathrm{lR}^m)`$). There are relations between these two conditions, in particular a smooth weakly harmonic map $`u`$ (i.e. critical for perturbations in the target) is automatically stationary (i.e. critical for perturbations in the domain). This is not true in the general case : there exists weakly harmonic maps which are not stationary (see \[HLP\]) and that can even be nowhere continuous (see \[Ri1\]). On the one hand, as a consequence of being weakly harmonic, one has the Euler Lagrange equation (harmonic map equation)

$$\mathrm{\Delta }u+A(u)(u,u)=0,$$
(1.1)
where $`A(u)`$ is the second fundamental form of $`N`$ embedded in $`\mathrm{lR}^k`$. Stationarity, on the other hand, implies the monotonicity formula saying that for any point $`x_0`$ in $`B^m`$ the density of energy $`r^{2m}_{B_r(x_0)}|u|^2`$ is an increasing function of $`r`$. Among stationary harmonic maps are the minimizing harmonic maps (minimizing $`E`$ for their boundary datas). In \[ScU\] R. Schoen and K. Uhlenbeck proved that there exists $`\epsilon (m,N)>0`$ such that, for any minimizing harmonic map, $`u`$ in $`W^{1,2}(B^m,N^n)`$ and any ball $`B_r(x_0)B^m`$ the following holds

$$\frac{1}{r^{n2}}_{B_r(x_0)}|u|^2𝑑x\epsilon (m,N)u_{L^{\mathrm{}}(B_{r/2}(x_0))}\frac{C}{r}$$

where $`C`$ only depends on $`m`$ and $`N^n`$. In other words, there exists a number $`\epsilon >0,`$ depending only $`m`$ and $`N`$, making a strict splitting between the almost linear regime for minimizing harmonic maps in which derivatives of the maps are under control and the totally non-linear regime where the map may be singular. Since this theorem in the early eighties, and previous partial regularity result of that type for minimal surfaces and other elliptic systems by Almgreen, DeGiorgi, Giusti, Miranda, Morrey, $`\mathrm{}`$, $`ϵ`$regularity theorems has been found in various problems of geometric analysis ( geometric flows, Yang-Mills Fields, $`\mathrm{}`$, etc). In particular it was a natural question whether the above result was true for arbitrary stationary harmonic maps. This was done 10 years later by L.C. Evans (\[Ev\]) in the case where $`N^n`$ is a standard sphere. Evans’ result was extended to general target by F. Bethuel in \[Be\]. Evans benefitted from a deeper understanding of the non-linearity in (1.1) developed in Hélein’s proof that any weakly harmonic map from a two dimensional domain is $`C^{\mathrm{}}`$ ( see \[He\] and \[CLMS\]). This fact will be discussed in the next subsection.

Combining the $`\epsilon `$regularity result with a classical Federer-Ziemer covering argument yields the following upper-bound on the size of the singular set of stationary harmonic maps. For a stationary harmonic map $`u`$ let Sing$`u`$ be the complement of the largest open subset where $`u`$ is $`C^{\mathrm{}}`$. Then

$$^{m2}(\text{Sing}u)=0,$$
(1.2)
where $`^{m2}`$ denotes the $`m2`$dimensional Hausdorff measure. For stationary maps and general targets, this is the best estimate available. For minimizing harmonic maps Schoen and Uhlenbeck proved the optimal result : dim(Sing$`u)m3`$. The reason for this improvement is the striking fact that weakly converging minimizing harmonic maps in $`W^{1,2}(B^m,N)`$ are in fact strongly converging in this space. This feature is very specific to the minimizing map case. When we blow up the map at a singular point strong convergence then implies roughly that the “singular set of the limit is not smaller than the limit of the singular set”.The Federer dimension reduction argument from the theory of minimal surfaces then gives by induction the result. Here one sees the strong connection between understanding the singular set and compactness properties of sequences of solutions and bubblings as presented in the next subsection.

### 1.2. Bubblings of harmonic maps

The following result of J.Sacks and K.Uhlenbeck \[SaU\] is perhaps the earliest example of energy quantization in non-linear analysis. Their aim was to extend the work of Eells and Sampson \[ES\] to targets of not necessarily negative sectional curvatures, that is : find in a given 2-homotopy class of an arbitrary riemannian manifold $`N^n`$ a more “natural” representant that would minimize the area and the Dirichlet energy $`E`$ in the given homotopy class . This raises the question of compactness and the possible reasons for the lack of compactness for harmonic maps from a 2-sphere (or even more generally a Riemann surface) into $`N^n`$. We should mention not only the original contribution of J.Sacks and K.Uhlenbeck that focused on minimizing sequences but also the later works by J.Jost \[Jo\] for critical points in general, by M. Struwe \[St\] for it’s heat flow version and the more recent contributions by T.Parker \[Pa\], W. Ding and G. Tian \[DiT\] and F.H.Lin and C. Wang \[LW\]. The following result has influenced deeply the non-linear analysis of the eighties, from the concentration-compactness of P.L.Lions to the compactness of J-holomorphic curves by M.Gromov and the analysis of self-dual instantons on 4-manifolds by S.K. Donaldson and K.Uhlenbeck.

Theorem 1. \[SaU\], \[Jo\] Let $`u_n`$ be a sequence of weakly harmonic maps from a surface $`\mathrm{\Sigma }`$ into a closed manifold $`N^n`$ having a uniformly bounded energy. Then a subsequence $`u_n^{}`$ weakly converges in $`W^{1,2}(\mathrm{\Sigma },N^n)`$ to a harmonic map $`u`$ into $`N^n`$. Moreover, there exist finitely many points $`\{a_1\mathrm{}a_k\}`$ in $`\mathrm{\Sigma }`$ such that the convergence is strong in $`\mathrm{\Sigma }\{a_1\mathrm{}a_k\}`$ and the following holds

$$|u_n^{}|^2dvol_\mathrm{\Sigma }|u|^2dvol_\mathrm{\Sigma }+\underset{j=1}{\overset{k}{}}m_j\delta _{a_j}\text{in Radon measure}$$

where $`m_j=_{i=1}^{P_j}E(\varphi _i^j)`$ and $`\varphi _i^j`$ are nonconstant harmonic 2-spheres of $`N^n`$ (harmonic maps from $`S^2`$ into $`N^n`$).

The loss of energy during the weak convergence is not only concentrated at points but is also quantized : the amount is given by a sum of energies of harmonic 2-spheres of $`N^n`$, the bubbles, that might sometimes be even explicitly known (for instance if $`N=S^2`$, $`E(\varphi _i^j)8\pi \mathrm{ZZ}`$). The striking fact in this result is that it excludes the possibility of losing energy in the neck between $`u`$ and the bubbles or between the bubbles themselves. This energy identity is quite surprising. It is a-priori conceivable, for instance, that, in a tiny annulus surrounding a blow-up point, an axially symmetric harmonic map into $`N`$, that is a portion of geodesic of arbitrarily small length in $`N`$, breaking the quantization of the energy, would appear. The energy identity disappears if instead of exact solutions we consider in general Palais-Smale sequences for $`E`$ in general (see the work of Parker \[Pa\] that was following his bubbling picture he established in collaboration with J.Wolfson for pseudo-holomorphic curves).

Only relatively recently a first breakthrough was made by F.H. Lin in the attempt of extending Sacks Uhlenbeck result beyond the conformal dimension.

Theorem 2. \[Li\] Let $`u_n`$ be a sequence of stationary harmonic maps from $`B^m`$ into a closed manifold $`N^n`$ having a uniformly bounded energy, then there exists a subsequence $`u_n^{}`$ weakly converging in $`W^{1,2}(B^m,N^n)`$ to a map $`u`$ and there exists a $`m2`$ rectifiable subset $`K`$ of $`B^m`$ such that the convergence is strong in $`B^mK`$ and moreover

$$|u_n^{}|^2dvol_\mathrm{\Sigma }|u|^2dvol_\mathrm{\Sigma }+f(x)^{m2}K\text{in Radon measure}$$

where $`f`$ is a measurable positive function of $`K`$.

This result, establishing the regularity of the blow-up set of weakly converging stationary harmonic map, is related to the resolution by D.Preiss \[Pr\] of the Besicovitch conjecture on measures admitting densities. Given a positive Borel regular measure $`\mu `$ such that there exists an integer $`k`$ for which $`\mu `$-a.e. the density $`lim_{r0}\mu (B_r)/r^k`$ exists and is positive, it is proved in \[Pr\] that there is a rectifiable $`k`$dimensional rectifiable set $`K`$ such that $`\mu (B^mK)=0`$. In the present situation calling $`\mu _n=|u_n|^2dvol_\mathrm{\Sigma }`$ converging to the Radon measure $`\mu =|u|^2dvol_\mathrm{\Sigma }+\nu `$, from the monotonicity formula one deduces easily that the defect measure $`\nu `$ fulfills the assumptions of Preiss theorem for $`k=m2`$ and the rectifiability follows at once. It has to be noted that the original proof of Theorem 2 in \[Li\] is self-contained and avoids D. Preiss above result.

Beyond the regularity of the defect measure $`\nu `$, the question remained whether the whole picture established in the conformal dimension could be extended to higher dimensions and if the energy identity still holds. F.H. Lin and the author brought the following answer to that question.

Theorem 3. \[LR2\] Let $`f`$ be the function in the previous theorem, if $`N^n=(S^n,g_{stand})`$ or if there is a uniform bound on $`u_n^{}_{W^{2,1}(B^m)}`$, then, for $`^{m2}`$ a.e. $`x`$ of $`K`$, there exists a finite family of harmonic 2-spheres in $`N^n`$ $`(\varphi _x^j)_{j=1\mathrm{}P_x}`$ such that

$$f(x)=\underset{j=1}{\overset{P_x}{}}E(\varphi _x^j).$$

The proof of this energy identity in higher dimension is of different nature from the one provided previously in dimension 2 where the use of objects relevant to the conformal dimension only, such as the Hopf differential, was essential. The idea in higher dimension was to develop a technique of slicing, averaging method combined with estimates in Lorentz spaces $`L^{2,\mathrm{}}L^{2,1}`$. This technique seems to be quite general as it have some genericity and permitted to solve problems of apparently different nature as we shall expose in the next section. The requirement of the $`W^{2,1}`$ bound for the case of a general target seems to be technical and should be removed. In the particular case when the target is the round sphere it was proved in \[He\] and \[CLMS\] that the non-linearity $`A(u)(u,u)`$ in equation (1.1) is in the Hardy space $`_{loc}^1`$ which immediately implies the desired $`W^{2,1}`$ bound for the maps $`u`$. Whether this fact and this $`W^{2,1}`$ bound can be extended to general targets is still unknown.

The 2 previous results suggest to view the loss of compactness in higher dimension as being exactly the one happening in the conformal $`2`$-dimensional case in the plane normal to $`K`$, locally invariant in the remaining $`m2`$ dimensions tangent to $`K`$. This understanding of the loss of compactness through creation of 2-bubbles for stationary harmonic maps has consequences in regularity theory. Indeed, in the case where $`N^n`$ admits no harmonic 2-spheres (take for instance $`N^n`$ a surface of positive genus), no blow-up can arise ; weakly converging stationary harmonic maps are strongly converging and the Federer dimension reduction argument discussed at the end of the previous subsection combined with the analysis in \[Si\] may be applied to improve the bound of the size of the singular set to dim Sing$`um4`$ (see \[Li\]).

### 1.3. High dimensional gauge theory

The work of S.K. Donaldson and R. Thomas \[DoT\] has given a new boost in the motivation for developing the non-linear analysis of high dimensional gauge theory.

Bubbling of Yang-Mills Fields

Theorems 1,2,3 and their proofs are transposable to many others geometric variational problems (such has Yang-Mills fields, Yamabe metrics…) having a given conformal invariant dimension $`p`$ (p=2 for harmonic maps, 4 for YM…etc). For instance in \[Ti\], G.Tian established the result corresponding to theorem 2 for Yang-Mills fields.

Consider a vector bundle $`E`$ over a Riemannian manifold $`(M^m,g)`$ and assume $`E`$ is issued from a principal bundle whose structure group is a compact Lie group $`G`$. Yang-Mills fields are connections $`A`$ on $`E`$ whose curvature $`F_A`$ solves

$$d_A^{}F_A=0$$
(1.3)
where $`d_A^{}`$ is the adjoint to the operator $`d_A`$, acting on $`^{}MEndE`$, with respect to the metric $`g`$ on $`M`$ and the Killing form of $`G`$ on the fibers (recall that the Bianchi identity reads instead $`d_AF_A=0`$). Yang-Mills fields are the critical points of the Yang-Mills functional $`|F_A|^2𝑑vol_g`$ for perturbations of the form $`A+ta`$ where $`a\mathrm{\Gamma }(^1MadQ)`$. It is proved in \[Ti\] that, taking a sequence of smooth Yang Mills connections $`F_{A_n}`$ having a uniformly bounded Yang-Mills energy, one may extract a subsequence (still denoted $`A_n`$) such that for some $`m4`$ rectifiable closed subset $`K`$ of $`M`$ the following holds : in a neighborhood of any point of $`MK`$ there exist good choices of gauges such that $`A_n`$, expressed in these gauges, converges in $`C^k`$ topology (for any $`k`$) to a limiting form $`A`$ that defines globally a smooth Yang-Mills connection on $`i^{}E`$, where $`i`$ is the canonical embedding of $`MK`$ in $`M`$. Moreover, one has the following convergence in Radon measure :

$$|F_{A_n}|^2dvol_M|F_A|^2dvol_M+f(x)^{m4}K,$$
(1.4)
where $`f`$ is some non-negative measurable function on $`K`$. For $`m=4`$, or under the assumption that the measures $`|_{A_n}_{A_n}F_{A_n}|dvol_M`$ remain uniformly bounded it is proved by the author in \[Ri4\] that $`f(x)`$ is, $`^{m4}`$ almost everywhere on $`K`$, quantized and equal to a sum of energies of Yang-Mills connections over $`S^4`$ : the energy identity holds. The proof use again the Lorentz space duality $`L^{2,\mathrm{}}L^{2,1}`$ applied this time to the curvature. These estimates are consequences of $`L^{4,2}`$ estimates on the connections $`A`$ for $`m=4`$. It is interresting to observe that such an $`L^{4,2}`$ estimate of the connection plays a crucial role in the resolution by J.Shatah and M. Struwe of the wave map Cauchy problem in 4 dimension for small initial datas in the “natural” space $`H^2\times H^1`$ (see \[ShS\]).

To complete the description of the blow-up phenomena several open questions remain both for harmonic maps and Yang-Mills fields. First, what are the exact nature of the limiting map $`u`$ and connection $`A`$ ? Is $`u`$ still a stationary harmonic map on the whole ball $`B^m`$ ? Is $`A`$ still a smooth Yang-Mills connection of some new bundle $`E_0`$ over $`(M,g)`$ ? Can one expect in both cases more regularity than the rectifiability for $`K`$ ? What is the exact nature of the concentration set $`K`$ ?

It happens that these two questions are strongly related to eachother : For instance, in the context of harmonic maps, if the weak limiting map $`u`$ was also a stationary harmonic map, then one would deduce from the monotonicity formula that $`K`$ is a stationary varifold and therefore inherits nice regularity properties. There is a similar notion of stationary Yang-Mills field (see \[Ti\]). The limiting $`A`$ being stationary, as expected, likewise would imply the stationarity of $`K`$. These questions in both cases are still widely open and seem to be difficult. Nevertheless in the case of Yang-Mills, in a work in preparation, T. Tao and G. Tian are proving a singularity removability result saying that the weak limit $`A`$ of smooth Yang-Mills $`A_n`$ can be extended to a smooth Yang-Mills connection aside from a $`^{m4}`$ measure zero set. Their proof is based on an $`ϵ`$-regularity result for weak Yang-Mills in high dimension which was also independently obtained by Y. Meyer and the author in \[MR\].

An important case where these questions have been solved is :

The case of $`𝛀`$-anti-self-dual Instantons

This notion extends the 4-dimensional notion of instantons to higher dimension. Assuming $`m4`$, and given a closed $`m4`$-form $`\mathrm{\Omega }`$, we say that a connection $`A`$ is an $`\mathrm{\Omega }`$anti-self-dual connection when the following equation holds

$$(F_A\mathrm{\Omega })=F_A.$$
(1.5)
In view of the Bianchi identity $`d_AF_A=0`$ and the closedness of $`\mathrm{\Omega }`$ it follows that $`d_A^{}F_A=0`$ that is $`A`$ is a particular solution to Yang-Mills equations. It is then shown in \[Ti\] that, if we further assume that the co-mass of $`\mathrm{\Omega }`$ is less than 1, then $`K`$ in (1.4) is a minimizing current in $`(M,g)`$ calibrated by $`\mathrm{\Omega }`$ (i.e. $`\mathrm{\Omega }`$ restricted to $`K`$ coincides with the volume form on $`K`$ induced by $`g`$) and therefore minimizes the area in it’s homology class.

A special case of interest arises when $`(M,g,\omega )`$ is a Calabi-Yau 4-fold with a global “holomorphic volume form” $`\theta `$ (i.e. $`\theta `$ is a holomorphic $`(4,0)`$form satisfying $`\theta \overline{\theta }=dvol_M`$). Let $`\mathrm{\Omega }`$ be the form generating the $`SU(4)`$ holonomy on $`M`$ (i.e. $`\mathrm{\Omega }`$ is a unit holomorphic section of the canonical bundle $`^{(4,0)}T^{}M`$) and assume $`M`$ is the product $`T\times V^3`$ of a 2-torus $`T`$ with a Calabi Yau 3-fold $`V^3`$, for a given $`SU(2)`$bundle E, $`T`$invariant solutions to (1.5) are pairs of connections $`B`$ on $`E`$ over $`V^3`$ and sections $`\varphi `$ of $`EndE`$ that solve a vortex equations (see 6.2.2 in \[Ti\]). The moduli space of these solutions should be related to the so-called holomorphic Casson invariants (see also R.Thomas work \[Th\]). The loss of compactness of these solutions arises along $`T`$ times holomorphic curves in $`V^3`$. To relate this moduli space with the holomorphic Casson invariants it is necessary to perturb, in a generic way, the complex structure into a not necessarily integrable one. This generates several non-linear analysis questions related to bubbling and regularity which have to be solved :

* Show that the weak limits of the vortex equations in the almost complex setting in $`V^3`$ have only isolated point singularities located on the pseudo-holomorphic blow-up set (see the extended conjectures on the singular sets of solutions to (1.5) in \[Ti\]).
* Construct solutions to the vortex equations which concentrate on some given choice of pseudo-holomorphic curves.
* Show that an arbitrary 1-1 rectifiable cycle (i.e. a 2 dimensional cycle whose tangent plane is invariant under the almost complex structure action) in an almost complex manifold (that may arise as a blow-up set in the above case) is a smooth surface aside from eventually isolated branched points.

In collaboration with G.Tian, in the direction of iii), the following result was established.

Theorem  \[RT\] Let $`(M^4,J)`$ be a smooth almost complex 4-real manifold. Then any 1-1 rectifiable cycle is a smooth pseudo-holomorphic curve aside from isolated branched points.

## 2. Ginzburg-Landau line vortices

### 2.1. The strongly repulsive asymptotic

The free Ginzburg-Landau energy on the 3-dimensional ball reads

$$GL(u,A)=_{B^3}|duiAu|^2+\frac{\kappa ^2}{2}(1|u|^2)^2+|dA|^2$$
(2.1)
where $`u`$ is a complex function satisfying $`|u|1`$ and $`A`$ is a 1-form on $`B^3`$. The study of the above variational problem was initiated by A. Jaffe and C. Taubes in 2 dimensions and by J. Fröhlich and M. Struwe in \[FS\] in 3 dimensions. The connected components of the zero set of the order parameter $`u`$ are the so called vortices that must generically be lines (for 0 being a regular value of $`u`$). The parameter $`\kappa `$ plays a crucial role in the theory. Depending on the value of $`\kappa `$, one expects different types of “behaviors” of vortices for minimizing configurations relative to various constraints (boundary datas, topological constraints…) (see \[JT\] and \[Ri3\] for a survey on these questions). In particular for large $`\kappa `$, in the so called strongly repulsive limit, the vortices of minimizing configurations tend to minimize their length (the first order of the energy is $`2\pi \mathrm{log}\kappa `$ times the length of the vortices) and repel one another, which is not the case for small $`\kappa `$. The magnetic field $`dA`$ plays no role in this mechanism ; to simplify the presentation we henforth simply impose that $`A=0`$. We are then looking at the functional

$$E_\epsilon (u)=_{B^3}|u|^2+\frac{1}{2\epsilon ^2}(1|u|^2)^2=_{B_3}e_\epsilon (u)$$
(2.2)
in the strongly repulsive limit $`\epsilon 0`$ whose critical points satisfy the non longer gauge invariant Ginzburg-Landau equation

$$\mathrm{\Delta }u+\frac{u}{\epsilon ^2}(1|u|^2)=0\text{ in }𝒟^{}(\mathrm{lR}^3).$$
(2.3)
Because we are interested in observing finite length vortices, in view of the above remark, we have to restrict to critical points satisfying $`E_\epsilon (u_\epsilon )=O(\mathrm{log}\frac{1}{\epsilon })`$. Equation (2.3) gives, in particular, that $`\mathrm{\Delta }u`$ is parallel to $`u`$. Therefore, assuming that $`|u|>\frac{1}{2}`$ in a ball $`B_r(x_0)`$ and writing $`u=|u|e^{i\varphi }`$, we deduce the nice scalar elliptic equation

$$div(|u|^2\varphi )=0\text{ in }𝒟^{}(B_r(x_0)).$$
(2.4)
Assuming $`|u_\epsilon |>1/2`$ on the whole $`B^3`$, the compactness of solutions to (2.3), in $`W^{1,2}`$ say, is reduced to the compactness of the boundary data. It is then clear that, in general, the study of the compactness of solutions $`u_\epsilon `$ to (2.3), in the limit $`\epsilon 0`$, involves the study of the compactness of the sets $`V_\epsilon =\{x;|u_\epsilon (x)|<1/2\}`$ as well as a control of the degree of $`u/|u|`$ on arbitrary closed curves in the complement of $`V_\epsilon `$ approximating the vortices. (The number $`1/2`$ may be of course replaced by an arbitrary number between 0 and 1).

The approximate vortices $`T_\epsilon `$

In their study of the 2-dimensional version of the present problem ($`B^3`$ replaced by $`B^2`$), F. Bethuel, H. Brezis and F. Hélein in \[BBH\], established that for solutions to (2.3) satisfying a fixed boundary condition from $`B^2`$ into $`S^1`$ of non-zero degree, $`V_\epsilon `$ can be covered by a uniformly bounded number of balls of radius $`O(\epsilon )`$ around which $`u/|u|`$ has a uniformly bounded topological degree. Taking then the distribution $`T_\epsilon `$ given by the sum of the Dirac masses at the center of the balls with the multiplicity given by the surrounding degrees of $`u/|u|`$, because of the uniform bound of the mass of $`T_\epsilon `$, we may extract a converging subsequence of atomic measures $`T_\epsilon ^{}`$. The compactness of the maps $`u_\epsilon `$ is a consequence of the convergence of $`T_\epsilon ^{}`$ via a classical elliptic argument.

Going back to 3 dimensions, the idea remains the same : we introduce an approximation of the vortices and try to prove compactness. More precisely, we consider a minimal 1-dimensional integer current $`T_\epsilon `$ in the homology class of $`H_1(V_\epsilon ,B^3,`$ $`\mathrm{ZZ})`$ given by the pre-image by $`u_\epsilon `$ of a regular point in $`B_{\frac{1}{4}}^2`$. Then, as in 2-dimensions, the key to the compactness of the critical points of $`E_\epsilon `$ is now the compactness of the familly of 1-dimensional rectifiable current $`T_\epsilon `$ constructed above.

The energy quantization result exposed in the following subsection shows the compactness of the $`T_\epsilon `$.

### 2.2. Quantization for G-L vortices : the $`𝜼`$-compactness

The following energy quantization result introduced by the author in \[Ri2\] says that there exists an absolute constant $`\eta >0`$ such that if the Ginzburg-Landau energy of a critical point in a given ball, suitabely renormalized, lies below this number, then there is no vortex passing through the ball of half radius. Precisely we have.

Theorem 5. \[LR1\] There exists a positive number $`\eta `$, such that for $`\epsilon `$ small enough and for any critical point of $`E_\epsilon `$ in $`B_1^3`$ satisfying $`|u|1`$, the following holds. Let $`B_r(x_0)B_{\frac{1}{2}}^3`$, $`r>\epsilon `$ then

$$\frac{1}{r}_{B_r(x_0)}e_\epsilon (u_\epsilon )\eta \mathrm{log}\frac{r}{\epsilon }|u|>\frac{1}{2}\text{ in }B_{\frac{r}{2}}(x_0).$$

The above result is optimal in the following sense : there are critical points and arbitrary $`r>\epsilon `$ such that $`\frac{1}{r}_{B_r(x_0)}e_\epsilon (u_\epsilon )2\pi \mathrm{log}\frac{r}{\epsilon }`$ and $`u(x_0)=0`$. Such a result is reminiscent of the $`\epsilon `$regularity result, although one difference here is that one has to handle blowing up energy $`E_\epsilon `$ which is not a-priori bounded as $`\epsilon `$ tends to 0.

The $`\eta `$compactness was introduced first for minimizers in 3-D in \[Ri2\] and then obtained for any critical points, still in 3-D, in \[LR1\]. Regarding dimension 2, the idea of the $`\eta `$compactness is implicit in the works before \[Ri2\]. Making this idea explicit in \[Ri5\] helped to substantialy simplify the existing proofs in the 2-D case.

It is a striking fact, whose explanation is beyond the scope of this lecture, that the original proof of theorem in \[LR1\] is based on the same techniques used by the authors to prove the energy identity (theorem 3) for harmonic maps. The key of the proof is to obtain a bound independent of $`\epsilon `$ for the density of energy $`\frac{1}{\rho }_{B_\rho (x_0)}e_\epsilon (u)`$ (for some $`\rho `$ between $`\epsilon `$ and $`r`$). The main obstacle in establishing such a bound is to control, independently of $`\epsilon `$, the part of the energy in $`B_\rho (x_0)V_\epsilon `$ coming from the vortices $`T_\epsilon `$ in $`B_\rho (x_0)`$. This is obtained by controlling, on a generic good slice $`B_{\rho _1}(x_0)V_\epsilon `$ ($`\rho /2<\rho _1<\rho `$), the $`L^{2,\mathrm{}}`$ and the $`L^{2,1}`$ norms of $`\psi `$, where $`\psi `$ is the vortex potential $`\psi :=\mathrm{\Delta }^1T_\epsilon `$. Recently in \[BBO\] F.Bethuel, H.Brezis and G.Orlandi made a nice observation by deriving, from the monotonicity formula associated to the elliptic variational problem $`E_\epsilon `$, an $`L^{\mathrm{}}`$ bound in terms of the local density of energy $`\frac{1}{\rho }_{B_\rho (x_0)}e_\epsilon (u)`$ for the vortex potential $`\psi `$ in $`B_\rho (x_0)V_\epsilon `$. This bound enabled them to replace in our original proof the use of the $`L^{2,\mathrm{}}L^{2,1}`$ duality for $`\psi \psi `$ on a slice by a more standard, in calculus of variation, $`L^{\mathrm{}}L^1`$ duality for $`\psi \mathrm{\Delta }\psi `$ on $`B_\rho (x_0)`$ itself. As explained in \[BBO\], this modification is a real improvement in the high dimensional case since it allows to end the proof of the $`\eta `$compactness without having to go through an argument based on a good slice extraction that would have been certainly more painful for $`n4`$. This modification of the duality used to control the energy of the vortex potential also allowed C. Wang recently to extend the $`\eta `$compactness result from \[LR3\] for the Ginzburg-Landau heat flow equation from 3 to 4 dimensions.

The $`\eta `$-compactness is a compactness result

Considering, as mentioned above, critical points $`u_\epsilon `$ to $`E_\epsilon `$ satisfying $`|u|1`$ and whose energy are of the order of $`\mathrm{log}\frac{1}{\epsilon }`$ ($`E_\epsilon (u_\epsilon )=O(\mathrm{log}\frac{1}{\epsilon }))`$, a Besicovitch covering argument (see \[LR1\]) combined with the $`\eta `$compactness above gives

$$V_\epsilon _{j=1}^{N_\epsilon }B_\epsilon (x_j)\text{ with }N_\epsilon =O\left(\frac{1}{\epsilon }\right).$$

This estimate, combined with the $`L^{\mathrm{}}`$ bound of $`u_\epsilon `$ deduced from the equation and the assumption $`|u|1`$, yields the following fact: $`M(T_\epsilon )=O(1)`$, (i.e. the mass of the approximated vortices is uniformly bounded). Combining this estimate with the fact that, by definition, $`T_\epsilon B^3=0`$, as a direct application of Federer-Fleming compactness theorem, we deduce the compactness of the approximated vortices $`T_\epsilon `$ in the space of integer rectifiable 1-dimensional currents.

Thus, the $`\eta `$compactness is the compactness of the approximated vortices $`T_\epsilon `$, which are rectifiable currents rather than a compactness of maps $`u_\epsilon `$. The latter compactness, in $`W^{1,p}`$ for $`p<3/2`$, nevertheless automatically follows from the previous one by the means of classical elliptic arguments for a suitable class of boundary conditions for $`u_\epsilon `$ (see \[LR1\] page 219). This class has been recently extended to arbitrary $`H^{\frac{1}{2}}(B^3,S^1)`$ conditions independent of $`\epsilon `$ in \[BBBO\]. This idea of getting convergence of maps going through convergence of currents is in the spirit of the theory of cartesian currents developed by M. Giaquinta, G. Modica and J. Soucek in \[GMS\].


format: text

text:
# Dynamical Stability of Earth-Like Planetary Orbits in Binary Systems

## 1 Introduction

Most stars have companions. Recent discoveries of extrasolar planetary systems have shown that Sun-like stars often have planetary mass companions and that these extrasolar planets reside in a wide variety of orbital configurations (Butler et al. 1999, Marcy et. al 2001; Fischer et al. 2002; Tinney et al. 2002; see also the extrasolar planet almanac111http://exoplanets.org/almanacframe.html). In addition, most solar type stars are known to reside in binary systems (Abt 1983) and thus often have stellar mass companions; the distributions of orbital parameters for these binaries are relatively well known (e.g., Duquennoy & Mayor 1991; hereafter DM91). If terrestrial planets are also present in these solar systems, the companions can affect their long term orbital stability. This paper addresses this issue of planetary stability with the overall goal of estimating the fraction of binary systems that allow an Earth-like planet to remain stable over the current age of our solar system.

Although terrestrial planets have not been detected in extrasolar systems (with main-sequence primaries) due to their small masses, they are expected to form in planetary systems alongside their Jovian counterparts (e.g., Ruden 1999; Lissauer 1993). The probable existence of such planets motivates this present work and underlies recent efforts to built the Terrestrial Planet Finder in the near future.222http://planetquest.jpl.nasa.gov/TPF/tpf\_index.html To study the stability of Earth-like planets in these systems, we must undertake a numerical investigation of the three-body problem. Although the dynamics of systems containing three gravitationally attracting bodies was considered over a century ago by Poincaré, an exact analytic solution to the general problem is not possible. The robust nature of these celestial systems is still being explored numerically. Indeed, it is now understood that multi-body gravitational systems generally exhibit chaotic behavior, thereby making it impossible to draw specific universal conclusions about their dynamics. A complete understanding of how a specific system evolves can only be garnered through a thorough numerical investigation. Furthermore, due to sensitive dependence on the initial conditions, the results of any numerical study must be presented statistically, in terms of the full distribution of outcomes resulting from effectively equivalent starting conditions. In this context, we consider equivalent starting conditions to be those with the same masses, eccentricities, and semi-major axes for the orbits, but with differing choices of initial phase angles (and other orbital elements – see Murray & Dermott 2000). As we illustrate in greater detail below, the ejection time displays a distribution of values for an ensemble of effectively equivalent initial conditions.

A great deal of work, both analytical and numerical, has already been done on stability (e.g., Szebehely 1980) and the development of chaos in celestial mechanics (e.g., Lecar, Franklin & Holman 2001). But many of the recent investigations of planetary stability focus on specific astronomical systems and are relatively narrow in scope (e.g., Benest 1996; Wiegert & Holman 1997; Laughlin & Adams 1999, 2000; Rivera & Lissauer 2000; Rivera & Haghighipour 2002; Dvorak et al. 2003). More general investigations of planetary stability in binary systems have been carried out (e.g., Rabl & Dvorak 1988; Holman & Wiegert 1999, hereafter HW99), but previous studies have not found the fraction of binaries that allow for stable Earth-like planets. In this paper, we re-examine the three-body problem by considering the stability of Earth-like planets in the presence of a companion with mass in the range 0.001 $`<M_C/(1M_{})<`$ 0.5. This mass range includes companions as small as Jupiter and as large as K stars. The portion of parameter space (the $`aϵ`$ plane) is chosen so we can combine the numerical stability results with the observed distributions of binary orbital parameters to estimate the fraction of systems that allow Earth-like planets to remain stable over long time scales.

The first result of this investigation is a determination of the ejection time for Earth-like planets over an extensive range of the companion’s initial orbital eccentricity $`ϵ`$ and semi-major axis $`a`$. In this context, the Earth-like planet is assumed to have the mass of Earth and starts in a circular orbit with a radius of 1 AU. To a reasonable approximation, we find that the ejection time $`\tau _{\mathrm{ej}}`$ depends primarily on the periastron distance $`R_{\mathrm{min}}=a(1ϵ)`$ of the companion (for a given companion mass $`M_C`$). Over the sampled regime of parameter space, the dependence of the mean ejection time $`\tau _{\mathrm{ej}}(R_{\mathrm{min}})`$ is well characterized by straight lines in a semi-log plot in the $`\tau _{\mathrm{ej}}R_{\mathrm{min}}`$ plane, even though the ejection time displays a distribution of values for a given value of periastron. The ejection time thus shows an exponential dependence on $`R_{\mathrm{min}}`$ over a particular range of periastron values. For higher mass companions, the plots also display an inner regime, where $`R_{\mathrm{min}}`$ is small and the ejection times are consistently short (a few hundred years – much shorter than any astrophysical \[or geological\] time scales of interest). Although this work focuses on the stability of Earth-like planets, scaling laws allow the results to be applied in a broader context.

The second result of this investigation addresses the possible habitability of extrasolar terrestrial planets (e.g., Kasting, Whitmire & Reynolds 1993; Rampino & Caldeira 1994). Approximately two-thirds of main-sequence stars (of solar type) are found in multiple systems (Abt 1983), and the binary frequency is even higher for pre-main-sequence stars (Ghez, Neugebauer, & Matthews 1993). These companions can disrupt the orbit of an Earth. Although this issue has been explored by several groups (e.g., Gehman, Adams, & Laughlin 1996; Laughlin, Chambers, & Fischer 2002; Menou & Tabachnik 2002), our numerical simulations shed further light on the subject. Our numerical results indicate that distant stellar companions (specifically, those with sufficiently large values of $`R_{\mathrm{min}}`$) will not disrupt the orbits of Earth-like planets. This work extends the range of parameter space studied previously and provides estimates for the fraction of binary systems that allow habitable planets.

This paper is organized as follows. The numerical methodology is described in §2, along with considerations of both Hill stability and scaling laws. The results from this suite of numerical experiments are presented in §3, along with an analytic characterization of the dependence of the mean ejection time $`\tau _{\mathrm{ej}}`$ on the periastron value $`R_{\mathrm{min}}`$. In §4, we use the numerically determined ejection times in conjunction with observed binary parameters to estimate the fraction of binary systems that can contain habitable planets. The paper concludes in §5 with a summary and discussion of these results, including an application to the stability of Earth-like planets in known extrasolar planetary systems.

## 2 Methodology

Although gravitationally interacting, three-body systems have been studied at length, their deceptively complicated nature has stymied efforts toward a complete solution. In this paper, we explore one special case of the three-body problem (the results can be scaled to other parameter choices – see below). Specifically, we consider the possible ejection of an Earth-like planet that starts in an initially circular orbit with radius 1 AU (around a primary star with one solar mass). Through long term dynamical interactions with an outer companion, the orbital elements of the Earth-like planet evolve, generally in chaotic fashion, until the planet is ejected from the system.

In order to explore this stability issue on intermediate time scales, we use two different types of numerical codes to perform a series of three-body experiments. For the first method, we use a symplectic mapping code (written for Laughlin & Adams 1999, following the lead of Wisdom & Holman 1991). Using a relatively small time step of 4 days, we can maintain high accuracy over the course of our 10 Myr integrations. In the second method, Newton’s equations of motion are integrated directly using a Bulirsh-Stoer (B-S) scheme (Press et al. 1992). Although direct integration is computationally more expensive, it is accurate and explicit. For the systems at hand, our B-S scheme incurs errors in relative accuracy of order 1 part in $`10^{11}`$ per total time step (where each time step in the three-body problem is variable, but has a typical value of about 10 days).

The Earth’s initial orbit is always set to be circular ($`ϵ_E`$ = 0) with radius $`r=a_E=1`$ AU. The mass $`M_C`$, eccentricity $`ϵ`$ and semi-major axis $`a`$ of the companion body are then specified for each run. For the vast majority of our numerical experiments, both the Earth and the companion are started in the same plane. (We briefly explore the case of non-coplanar orbits in follow-up simulations.) With given values for the masses, semi-major axes, and eccentricities, the initial conditions for the simulations must also specify the orbital phases. We consider numerical experiments with the same binary orbital parameters ($`M_C,a,ϵ`$) and a random distribution for the remaining phase angles (see Murray & Dermott 2000).

For each set of initial conditions, we integrate the system forward in time. For the sake of definiteness, and in order to cover a large range of parameter space, we use ten million years as the upper limit for the integration time. These experiments give us either a time scale for instability or a lower limit of ten million years on the possible ejection time. The planet is considered to be ejected if any of the following conditions are met: The energy of the planet becomes positive; the eccentricity of the planet exceeds unity; the periastron of the planet becomes smaller than the stellar radius (assumed to be 1 $`R_{}`$) so the planet is accreted; or the semi-major axis of the planet exceeds a maximum value (taken here to be 100 AU).

For numerical experiments with the same binary properties ($`M_C,a,ϵ`$), the ejection time $`\tau _{\mathrm{ej}}`$ varies with the choices of orbital phases. Because the systems are chaotic, this variation is not smooth, i.e., small differences in the starting phase angles can lead to large differences in the resulting ejection times. An exploration of parameter space shows that the ejection time displays a log-normal distribution for an ensemble of different phases, i.e., different realizations of the same underlying problem (see Figure 1). We have found the distribution of ejection times using both the symplectic and the B-S code and find the same distribution. Figure 1 shows the resulting distribution for the case with binary companion mass $`M_C`$ = 0.1 $`M_{}`$, eccentricity $`ϵ`$ = 0.5, and semi-major axis $`a`$ = 5 AU. Also shown is a normal (gaussian) distribution (in $`\mathrm{log}\tau _{\mathrm{ej}}`$) with the same mean and width. Figure 1 illustrates several properties of this stability problem that guide the rest of this investigation: First, the width of the distribution is substantial ($`\sigma `$ = 0.51 for the case shown here). Second, both numerical methods give essentially the same result. Third, the distribution is log-normal so that we average our ejection times in $`\mathrm{log}\tau _{\mathrm{ej}}`$ for the remainder of the paper.

In addition to the two main numerical codes, a limited number of our results were verified and extended using the MERCURY symplectic integration package (developed by Chambers 1999). First, twelve sets of simulations were repeated for a Jupiter-mass companion, with each set containing three runs at constant periastron distance $`R_{\mathrm{min}}`$ (with varying eccentricities and semi-major axes); these runs were done to make sure that the different codes give the same results. Next, taking advantage of the increased speed of the MERCURY code, we used it to extend a limited number of runs out to integration times of 100 Myr. A more detailed examination was then performed for a range of systems that remained stable over a 100 Myr time scale by extending the integration time to 500 Myr and sampling over a wider range of Jupiter’s initial eccentricity.

The numerical results obtained here can be compared with analytical results obtained previously. As a reference point, we use standard theory to specify the condition for Hill stability of the system (following Gladman 1993). We first define the dimensionless quantities

$$\gamma _j=[1ϵ_j^2]^{1/2},\eta _j=\frac{m_j}{m_1+m_2},\mu _{12}=\frac{m_1+m_2}{(1M_{})},$$
(1)
where the subscripts $`j`$=1,2 refer to the Earth-like planet and the companion, respectively. The condition for Hill stability can then be written

$$[\eta _1+\eta _2/a_2][\eta _1\gamma _1+\eta _2\gamma _2\sqrt{a_2}]^2>1+3^{4/3}\eta _1\eta _2\mu _{12}^{}{}_{}{}^{2/3},$$
(2)
where $`a_2`$ is the semi-major axis of the companion, expressed in units where the semi-major axis of Earth is unity. In the present context, we start Earth with an initially circular orbit so that $`ϵ_1=0`$ and $`\gamma _1=1`$. In the limit $`m_1m_2`$, the above expression simplifies to the form

$$\gamma _2=[1ϵ_2^2]^{1/2}>\eta _2^{3/2}\frac{1}{2}a_2\eta _1\eta _2^{5/2}\eta _1a_2^{1/2}\eta _2^1.$$
(3)
Saturating this bound, we thus obtain the eccentricity $`ϵ_2`$ required for the system to become Hill unstable (for a companion with a given mass $`\eta _2`$ and semi-major axis $`a_2`$).

This criterion predicts that three-body systems will become Hill unstable for modest values of the eccentricity. For example, with a Jupiter-like companion, the Earth becomes Hill unstable for $`ϵ_2>0.1`$, i.e., a periastron distance of 4.7 AU. As many previous authors have found for this regime of parameter space (e.g., Valsecchi et al. 1984; Milani and Nobili 1983; Gladman 1993), systems that are unstable according to the criterion  may live for an extraordinarily long time (e.g., longer than the age of the universe) before showing any signs of instability (e.g., see also HW99; Levison, Lissauer, & Duncan 1998; Duncan & Lissauer 1998; Wisdom & Holman 1991; Laughlin & Adams 1999). In this context, if the system is unstable according to the Hill criterion , then it has a chance to decay (by ejecting a planet) in the long term. The time scale for decay can be described in terms of a half-life (e.g., Adams & Laughlin 2003); in a large sample of equivalent solar systems, half of the systems will decay by ejecting a planet in a well-defined time, but any particular system could decay over a wide range of times. If the half-life (or expectation value of the decay time333The expectation value is the mean decay time averaged over the probability distribution. For an exponential decay law, the expectation value $`\tau `$ is related to the half-life $`t_{1/2}`$ via $`t_{1/2}`$ = $`\tau \mathrm{ln}2`$.) is longer than the expected lifetime of the star (typically billions of years) or the age of the universe (12–14 Gyr), then the system can be considered as stable for evaluating the prospects for the survivability of Earth-like planets. In this paper, we study systems that are Hill unstable with half-lifes in the approximate range $`10^2`$ yr $`<t_{1/2}<`$ $`10^9`$ yr.

The results of this investigation can be scaled to other parameter choices. Here, we concentrate on the case of an Earth-like planet and fix its starting radius at $`r=a=1`$ AU. The equations of motion, however, have a radial scale invariance. As illustrated by the Hill stability criteria in equations , all length scales in the problem can be scaled by the semi-major axis $`a_1`$ of the Earth’s orbit. If we rescale the problem by changing $`a_1`$ by a factor $``$, then the resulting dynamics are the same, except that the time scales differ by a factor of $`^{3/2}`$. Similarly, we perform simulations using a value of 1.0 $`M_{}`$ for the mass of the primary. If we rescale the problem for an arbitrary mass $`mM_{}/(1.0M_{})`$, the calculated time scales change by a factor of $`m^{1/2}`$ and are applicable to companion masses that are rescaled according to $`M_CmM_C`$. The mass of the ‘Earth’ that is being modeled by the simulations also changes such that $`M_EmM_E`$. Because $`M_EM_C,M_{}`$, however, the Earth acts essentially like a test particle and its mass is of little consequence.

## 3 Results

Using the methodology outlined above, we have performed a large number of simulations of the Sun-Earth-companion system using both the symplectic and B-S numerical integration schemes. The long term evolution of these numerical experiments follows the same general trend. Over the first several thousand years, the companion drives Earth into an orbit with ever higher eccentricity. Because the companion is much more massive than Earth, its orbital elements change far less than those of the planet. In addition, Earth’s orbital eccentricity does not show a slow and steady increase, but rather shows a cyclic and often chaotic pattern (for examples of this type of behavior, see Laughlin & Adams 1999; Rivera & Lissauer 2000; Lecar et al. 2001). As a result, Earth’s eccentricity is best described as a distribution of possible values. As the Earth continues to orbit its star, it samples all of the eccentricity values in this distribution, but the distribution itself evolves with time. Figure 2 illustrates this behavior for one particular case near the center of our parameter space, i.e., a binary with companion mass $`M_C`$ = 0.1 $`M_{}`$, eccentricity $`ϵ`$ = 0.5, and semi-major axis $`a`$ = 5 AU. The figure shows the distributions of eccentricity for Earth near the beginning of the integration, at two intermediate times, and near the ending of the integration (just before ejection of the Earth-like planet). In general, the Earth samples the distribution more quickly than the distribution changes; this result is consistent with the Lyapunov time being correlated with (but much shorter than) the ejection time (see Lecar, Franklin, & Murison 1992). Over longer times, however, the distribution of eccentricities grows broader and samples larger and larger values of $`ϵ`$. When the distribution becomes sufficiently wide, the Earth stands a fair chance of entering into an orbit of extremely high eccentricity, which ultimately leads to instability. The Earth then either plunges into the star or is ejected from the solar system altogether.

Using the symplectic integration scheme, we have performed an ensemble of approximately 35,000 simulations of the Sun-Earth-companion system. The semi-major axis of the companion is chosen to vary in even increments from 1 AU (where Earth is ejected almost immediately) out to 15 AU (for the 1 $`m_J`$ companion), 20 AU (for the 10 $`m_J`$ companion), and 80 AU (for the 0.1 and 0.5 $`M_{}`$ stellar mass companions). For each choice of semi-major axis, the eccentricity is chosen to sample a range of periastron values from $`R_{\mathrm{min}}`$ = 1 AU out to values so large that the ejection time always exceeds the 10 Myr range of our integrations. For each choice of companion mass, eccentricity, and semi-major axis, the remaining phase angles are chosen from a random distribution. This sampling of parameter space fills a broad band in the $`aϵ`$ plane, as shown by the gray scale plots in Figures 3 – 6. For the portion of the $`aϵ`$ plane to the upper left of the chosen region, the companion crosses the Earth orbit and would lead to rapid ejection. For the portion of the $`aϵ`$ plane to the lower right of the sampled region, the ejection times are longer than 10 Myr.

The results obtained using the symplectic code are presented by ascending order of companion mass in Figures 3 – 6. The top panel of each figure depicts a grayscale plot of ejection times in the $`aϵ`$ plane, where the darker shades correspond to longer ejection times. Contour lines of constant periastron $`R_{\mathrm{min}}=a(1ϵ)`$ are included for reference. A natural dependency of the ejection time on $`R_{\mathrm{min}}`$ is clearly illustrated by the grayscale plots, with the largest gradients occurring perpendicular to contour lines. To elucidate this trend, we plot the ejection time versus $`R_{\mathrm{min}}`$ for each companion mass in the corresponding lower panel. Although the ejection time depends most sensitively on the periastron value, the relevant initial conditions that describe the binary continue to be the eccentricity and semi-major axis. In the Figures, the star symbols represent the average value of $`\mathrm{log}\tau _{\mathrm{ej}}`$ (for each given value of $`R_{\mathrm{min}}`$) and the vertical bars represent one standard deviation about the mean. The width of the distribution (at constant periastron) arises from two sources: (A) For a given type of binary (i.e., given $`ϵ`$ and $`a`$), the ejection time has a distribution of values as shown in Figure 1. (B) For binaries with differing values of $`(a,ϵ)`$, but constant periastron $`R_{\mathrm{min}}=a(1ϵ)`$, the ejection time has additional variation (e.g., see HW99). Although the ejection times vary by an order of magnitude at a given value of $`R_{\mathrm{min}}`$, a clear functional dependence for the mean value can nonetheless be extracted from the results. Over the sampled range of periastron values $`R_{\mathrm{min}}`$, the mean ejection time follows a straight line in the semi-log plots. In the following section, we will use this property to extrapolate our results out to longer ejection time scales.

We note that the importance of the periastron distance has been found in related dynamical investigations. In a study of scattering of Trans-Neptunian objects by the planet Neptune (Holman, Grav, & Gladman 2001), the borders of the region where the bodies become chaotic are approximately described by lines of constant perihelia (that study also discusses the departures from this trend). In a related work (Duncan, Levison, & Budd 1995), objects with perihelia less than about 35 AU were found to be unstable, apparently due to cumulative effects of random forces from Neptune (exerted near perihelion). Finally, the ejection time in systems of terrestrial planets can be modeled with an equation similar to our equation  (see Chambers, Wetherill, & Boss 1996).

In an alternate set of simulations, we have performed an ensemble of approximately 8000 numerical integrations using the B-S code. The results from the sympletic code suggested that the periastron is the most important variable for determining the ejection time. For the B-S integrations, we chose even increments in periastron $`R_{\mathrm{min}}`$ of the companion. For each value of $`R_{\mathrm{min}}`$, we sample the eccentricity over the full range $`0ϵ1`$ and perform multiple realizations of the experiments using different (random) choices for the remaining phase angles. This coverage of parameter space is complementary to that used for the symplectic code. The symplectic code was used to sample as much of the $`aϵ`$ plane as possible. The B-S code was used to sample fewer values of $`(a,ϵ)`$, and fewer values of the periastron, but each point in the $`aϵ`$ plane was studied with more realizations of the various phase angles in the problem. The ejection time, as a function of periastron, follows the same general trend for both ensembles of numerical experiments, as shown in the bottom panels of Figures 3 – 6 (where results from the B-S code are plotted as triangular symbols). The main difference between the two explorations of parameter space is that the latter displays a somewhat wider distribution of ejection times for a given periastron value. However, the expectation values of the distributions are in good agreement.

To package these results in a useful format, we fit the numerically determined mean ejection times with an empirical relation of the form

$$\tau _{\mathrm{ej}}=\tau _{\mathrm{ej}}^{}{}_{0}{}^{}\mathrm{exp}\left[\alpha (p1)\right],$$
(4)
where $`\tau _{\mathrm{ej}}^{}{}_{0}{}^{}`$ is a fiducial time scale, $`p`$ is the dimensionless periastron distance $`p(R_{\mathrm{min}}`$/1 AU), and $`\alpha `$ is the dimensionless fitting parameter (the slope of the lines in Figures 3 – 6). The variables for the functional fits are presented in Table 1 below, for the four companion masses $`M_C`$ considered here. The Table also lists the range in $`R_{\mathrm{min}}`$ over which the fitting formulae have been numerically determined. For the stellar mass companions, the ejection time is extremely short (hundreds of years) for a small range of periastron near 1 AU; the fit to equation  thus begins at larger periastron values.

Table 1: Parameters for ejection time scaling laws

| mass $`M_C`$ ($`M_{}`$) | $`\tau _{\mathrm{ej}}^{}{}_{0}{}^{}`$ (yr) | $`\alpha `$ | periastron range |
| --- | --- | --- | --- |
| 0.001 | 1800 | 9.8 $`\pm `$ 3.0 | 1 AU $`<R_{\mathrm{min}}<`$ 1.8 AU |
| 0.01 | 400 | 6.8 $`\pm `$ 1.1 | 1 AU $`<R_{\mathrm{min}}<`$ 2.1 AU |
| 0.10 | 110 | 4.1 $`\pm `$ 0.67 | 1.6 AU $`<R_{\mathrm{min}}<`$ 3.5 AU |
| 0.50 | 0.64 | 4.7 $`\pm `$ 0.32 | 2.4 AU $`<R_{\mathrm{min}}<`$ 4.4 AU |

Although the ejection time depends most sensitively on the variable $`R_{\mathrm{min}}`$, this time scale $`\tau _{\mathrm{ej}}`$ is basically a function of $`a`$ and $`ϵ`$. Close inspection of Figures 3 – 6 shows that nearby points in the $`aϵ`$ plane can lead to rather different ejection times. In addition, systems with the same values of $`a`$ and $`ϵ`$ can display differing ejection times for varying starting phases of the orbiting bodies – as illustrated in Figure 1. These systems are thus chaotic, in the technical sense, and display sensitive dependence on their initial conditions. This variation leads to a range of values for $`\tau _{\mathrm{ej}}`$, as quantified by the error bars in the lower panels of Figures 3 – 6. For a given value of $`R_{\mathrm{min}}`$, the ejection time $`\tau _{\mathrm{ej}}`$ displays a full distribution of values. As mentioned earlier, the width of this distribution arises from two sources. First, the distribution of ejection times has an intrinsic spread for a given pair $`(a,ϵ)`$, as shown in Figure 1. Second, we are averaging over many pairs $`(a,ϵ)`$ with the same periastron; since the ejection time can depend on both $`a`$ and $`ϵ`$ (e.g., HW99), this averaging increases the width of the distribution as shown by the error bars in the lower panels of Figures 3 – 6. The fitting formulae found above (see Table 1) provide an estimate for the expectation value of this distribution as a function of $`R_{\mathrm{min}}`$ (where the expectation value is the mean value, averaged over the underlying probability distribution of ejection times).

The results shown in Figures 3 – 6 are in basic agreement with previous work. In a study of Alpha Centauri, test particles were found to be stable for initially circular orbits within about 3 AU of either star (Wiegert & Holman 1997), where the integration time was 32,000 binary periods. In this binary, the mass ratio is close to 0.5, the periastron distance is 11.2 AU, and the semi-major axis is 23 AU. Scaled to the units of this paper, the boundary of stability corresponds to a time scale of $`\tau 4.8\times 10^5`$ yr for $`p3.7`$. This point falls just above the best fit line shown in the bottom panel of Figure 6, but well within the allowed range. In a more general study, HW99 found the critical radius $`a_{cr}`$ for which a planet will remain stable for $`10^4`$ binary periods. For a given binary mass ratio, the critical radius is a function of eccentricity (see Table 3 of HW99). When scaled to the Sun-Earth-companion systems considered here (with one year planetary orbits), the function $`a_{cr}(ϵ)`$ can be converted to a function $`\tau _{\mathrm{ej}}(p)`$, where $`\tau _{\mathrm{ej}}`$ corresponds to $`10^4`$ binary orbits. The resulting function $`\tau _{\mathrm{ej}}(p)`$ includes only one eccentricity value for a given value of $`p`$, rather than an average over a path in the $`aϵ`$ plane as we use here. Nonetheless, the RMS value of the relative difference444This RMS difference is calculated using equation (1) of HW99 and averaging over the range of eccentricity $`0ϵ0.8`$, i.e., the range listed in their Table 3. (in $`\mathrm{log}\tau `$) between the HW99 result and that predicted by equation  is only about 7% for a companion mass $`M_C`$ = 0.1 $`M_{}`$ and 12% for $`M_C`$ = 0.5 $`M_{}`$.

The results from our lowest companion mass simulations were confirmed for a limited number of cases using the MERCURY integration package (Chambers 1999). For systems with a Jupiter-mass companion, integrations were performed for periastron $`R_{\mathrm{min}}`$ ranging from 1.04 – 2.60 AU. For some of these simulations, the companion began in the same orbital plane as the Sun-Earth system (as before). We also ran an alternate series of simulations in which Jupiter began with an inclination of 1.3 relative to the Sun-Earth orbital plane; for comparison purposes, identical sets of simulations were also performed treating the Earth as a mass-less test particle. In all of the simulations, the other orbital elements were chosen at random. Each system’s evolution was then followed for 100 Myrs or until the Earth was lost from the system.

For simulations with a starting inclination angle for Jupiter, the results for periastron $`R_{\mathrm{min}}`$ = 1.04 AU show an interesting departure from the co-planar simulations. These simulations (where Jupiter starts with $`i`$ = 1.3) display a wider range of ejection times, but the average value remains essentially the same. Although the orbit of the companion brings it alarmingly close to Earth’s orbit, some of the simulations show relatively long ejection times. For example, in a simulation with Jupiter’s $`R_{\mathrm{min}}`$ = 1.04 AU and $`ϵ`$=0.83, after 20,000 years, the (massive) Earth has been perturbed into an orbit with $`a`$ = 9.2 AU, $`ϵ`$ = 0.6, and an inclination of $`i`$ = 139. The Earth’s eccentricity and inclination continue to oscillate for several hundred thousand years until the Earth is ejected from the system.

The MERCURY results verify that the range 1.6 AU $`<R_{\mathrm{min}}<`$ 1.8 AU marks the transition between instability and stability. Most of the systems with $`R_{\mathrm{min}}<`$ 1.7 AU lose the Earth in less than 10 Myr, whereas systems with $`R_{\mathrm{min}}>`$ 1.9 AU remain stable for 100 Myr. The stability of the high $`R_{\mathrm{min}}`$ systems was examined further through a series of 500 Myr runs using $`R_{\mathrm{min}}`$ = 1.87 – 2.01 AU and Jupiter’s eccentricity $`ϵ_J`$ = 0.1 – 0.8. Many of the systems remained stable for the duration of the integration, although there were several exceptions: With $`R_{\mathrm{min}}`$ = 1.87 and $`ϵ`$=0.1, for example, the Earth was ejected from the system in only 0.3 Myr. However, the other systems with $`ϵ`$ = 0.1 and $`R_{\mathrm{min}}>1.87`$ remained stable for the full 500 Myr, as did the systems with 0.1 $`ϵ`$ 0.7. On the other hand, half of the systems with $`ϵ`$ = 0.8 became unstable, with Earth being ejected from the system or accreted by the star. Taken together, these results are consistent with the ejection times predicted by the fit of equation , although the ejection time may be even longer (at large values of periastron) than that predicted by the fitting function. A larger ensemble of long term integrations is needed to clarify this issue. One should also keep in mind that equation  represents the expectation value of the distribution; for a given value of $`R_{\mathrm{min}}`$, the distribution has a width of nearly an order of magnitude in ejection time $`\tau _{\mathrm{ej}}`$.

## 4 Fraction of binary systems that allow Earths

This set of numerical experiemnts can be used to address the question of habitable planets. Although habitability includes many aspects, a key requirement is for the planet to remain stable over long time intervals. The majority of stars reside in binary systems and these companions can preclude the possibility of an Earth-like planet, i.e., a small planet in a 1 AU (nearly circular and stable) orbit about a solar type star. The results of the previous section indicate that the survival time of such an Earth increases with the periastron distance $`R_{\mathrm{min}}`$. In this section, we extrapolate our numerical results to longer times and find an approximate lower bound on the fraction of binary star systems that allow for habitable planets.

The numerical results of the previous sections are complete out to ejection times of 10 million years and include limited results out to 500 million years. Although the time required for life to develop on a planet is largely unknown, the most familiar habitable planet – our Earth – has an age of 4.6 Gyr. Since we want to find a conservative estimate for the fraction of binary systems that allow for habitability, we assume that an Earth-like planet must remain stable for 4.6 Gyr.

The ejection times for a given value of $`R_{\mathrm{min}}`$ span a wide range, as shown in Figures 3 – 6. Again adopting a conservative approach, we assume that the ejection times take the shortest values within their allowed range, and hence the ejection times are about 10 times shorter than the numerically determined expectation values (as fit by equation ). Our estimate for the minimum periastron distance required for survival over a time $`\tau _{SS}`$ thus takes the form

$$p>1+\alpha ^1\mathrm{ln}[10\tau _{SS}/\tau _{\mathrm{ej}}^{}{}_{0}{}^{}],$$
(5)
where $`p`$ is the periastron $`R_{\mathrm{min}}`$ in units of AU. For the sake of definiteness, we take $`\tau _{SS}`$ = 4.6 Gyr, the current age of the solar system. In adopting the form , we are assuming that equation  continues to hold out to larger periastron values and longer ejection times. Since our (limited) longer term integrations indicate that Earth-like planets may survive even longer than predicted by an extrapolation of equation , the limit implied by equation  represents a conservative bound. The parameters $`\tau _{\mathrm{ej}}^{}{}_{0}{}^{}`$ and $`\alpha `$ are given in Table 1 for companions of various masses.

We note that this constraint can be generalized. If we consider an Earth-like planet with circular orbit of radius $`\stackrel{~}{a}`$ = $`a`$/(1AU) in orbit about a primary star of mass $`m=M_{}/(1M_{})`$, and we want to enforce stability for a time $`\tau `$, the required constraint becomes

$$p>\stackrel{~}{a}\left\{1+\alpha ^1\mathrm{ln}[(10\tau /\tau _{\mathrm{ej}}^{}{}_{0}{}^{})m^{1/2}\stackrel{~}{a}^{3/2}]\right\}.$$
(6)
The values of $`\alpha `$ and $`\tau _{\mathrm{ej}}^{}{}_{0}{}^{}`$, as listed in Table 1, apply to rescaled companion masses $`\stackrel{~}{M}_C`$ = $`mM_C`$ (where $`M_C`$ are the masses given in the table).

Now we specialize back to our standard case of a 1.0 $`M_{}`$ primary and an initial Earth orbit of 1 AU. For companion masses $`M_C`$ = 0.1 and 0.5 $`M_{}`$, the minimum value of the periastron distance required for stability (according to equation ) is about $`p=67`$, as derived from the range of allowed slopes listed in Table 1. To be conservative, once again, we use the high end of this estimated range.555We note that this limit is consistent with a related result (HW99). Assuming that our solar system had a solar mass companion, HW99 found the minimum semi-major axis required for stability of the planets. For $`ϵ`$ = 0.4, they found that $`a>400500`$ AU. Since the outermost planet is Neptune with $`a_{\mathrm{Nep}}`$ = 30 AU, this result corresponds to a dimensionless periastron in the range $`p=5.36.6`$, consistent with the constraint discussed here. To find a lower limit on the fraction of binaries that allow for habitable Earths, we thus need to find the fraction of binary systems with dimensionless periastron $`p>7`$.

The observed population of binaries has a well-defined period distribution (DM91), which takes a log-normal form. Focusing on the case of primary stars with masses of approximately 1.0 $`M_{}`$, we can convert the period distribution into the probability distribution for semi-major axes $`a`$. The resulting normalized distribution thus takes the form

$$dP_a=f(\mathrm{ln}a)d\mathrm{ln}a=\frac{1}{\sqrt{2\pi }\sigma }\mathrm{exp}\left[\frac{(xx_0)^2}{2\sigma ^2}\right]dx,$$
(7)
where the variable $`x`$ is the logarithm of the dimensionless semi-major axis. The parameters $`x_0`$ and $`\sigma `$ are determined from fits to the observed binary period distribution (as reported in DM91). The dimensionless width $`\sigma `$ = 3.53. To convert the period distribution to a distribution of semi-major axis, we must include the ratio $`\mu `$ of the companion mass to the primary mass through the relation $`x_0`$ = 3.44 $`+(1/3)\mathrm{ln}[1+\mu ]`$. As a result, the distribution depends (weakly) on the value of the mass ratio $`\mu `$.

The distribution of eccentricity is also well defined, but the observed distribution takes a different form for close, intermediate, and wider binaries (DM91). In the present context, we are interested in binaries with periastron distances greater than about 7 AU, so the semi-major axes are also greater than 7 AU and we can use the eccentricity distribution for wider binaries. In this regime, the eccentricity distribution has the simple form $`dP_ϵ=2ϵdϵ`$ (DM91). Furthermore, for these wider binaries, the eccentricity is independent of semi-major axis (DM91, Heacox 1998).

For a given mass ratio $`\mu `$, we define $`F_\mu (p)`$ to be the fraction of binary systems with dimensionless periastron greater than $`p`$. To evaluate $`F_\mu (p)`$, we must integrate over the portion of the $`aϵ`$ plane with periastron $`a(1ϵ)>p`$, where the integrand is weighted by both the eccentricity distribution and the distribution of semi-major axes (described above). The eccentricity integration can be done analytically and the remaining integral – which defines the fraction $`F_\mu (p)`$ as a function of $`p`$ – takes the form

$$F_\mu (p)=\frac{1}{\sqrt{2\pi }\sigma }_{\mathrm{ln}p}^{\mathrm{}}𝑑x\left(1p\mathrm{e}^x\right)^2\mathrm{exp}[(xx_0)^2/2\sigma ^2].$$
(8)
The fraction $`F_\mu (p)`$ is a slowly varying function of the mass ratio $`\mu `$. To estimate the fraction $`F(p)`$ of all binary systems for which Earth-like planets can remain stable, we must find the weighted average of $`F_\mu (p)`$, i.e.,

$$F(p)_0^1w(\mu )F_\mu (p)𝑑\mu ,$$
(9)
where the distribution of mass ratios $`dP_\mu =w(\mu )d\mu `$. This distribution has been observed and is presented in DM91 (see Table 7 and Figure 10). The peak of the distribution occurs near $`\mu `$ = 0.3 and the majority of systems have mass ratios in the range $`0.1\mu 0.5`$. We use DM91 to specify the distribution $`w(\mu )`$ of mass ratios and evaluate the integral of equation . The resulting function $`F(p)`$ is shown in Figure 7.

Figure 7 also shows a fit to the numerically determined result. The fitting function is chosen to have the simple form

$$\stackrel{~}{F}(p)=F_1\mathrm{exp}\left[(a\xi +b\xi ^2)\right],$$
(10)
where $`F_1`$ = 0.711, $`a`$ = 0.101, $`b`$ = 0.0287, and $`\xi =\mathrm{ln}[p]`$. The function  provides a good approximation to the numerically determined result, with an absolute error less than about 0.011 and a relative error less than 4 percent. More exact fits are not warranted, as the observed distributions of binary orbital parameters are not known to this accuracy. This function can be used to estimate the fraction of binary systems with periastron greater than any specified value within the allowed range $`1<p<10^5`$.

As argued above, stability of Earth over 4.6 Gyr requires a stellar companion to have $`p>7`$ ($`R_{\mathrm{min}}>7`$ AU), where this value has been estimated from numerical experiments using companion masses in the range $`M_C=0.10.5M_{}`$ (which are typical values – see DM91). The fraction of binary systems that meet this constraint and allow habitable Earth-like planets is 0.5 (or 50 percent). This estimate should be regarded as a lower bound on the fraction. Additional binary systems could allow for habitable Earths if the orbits can be inclined. Notice also that this estimate does not include close binary systems, those with separations $`a1`$ AU. Although sufficiently close binary systems could allow Earth to remain stable (see, e.g., HW99), the companion will necessarily reside within 1 AU of the Earth and the companion’s radiative flux could affect considerations of habitability. For the wider binaries considered here, the requirement of dynamical stability demands that $`R_{\mathrm{min}}>7`$ AU so that the companion is always farther than about 6 AU from Earth. The additional radiative flux of the secondary is always less than about 3% of that of the primary. For a more typical companion mass, say $`M_C=0.4M_{}`$, the flux from the secondary is less than 0.2% of the total.

The numerical results of the previous section show that Sun-Earth-companion systems can remain stable for relatively long times – such as the current 4.6 Gyr age of our solar system – even though they do not meet the analytic criteria for stability (see equations ). To illustrate this point, Figure 8 shows the allowed region of the $`aϵ`$ plane for Earth-like planets in binary systems. The solid curves delineate the region of the plane that allows Earth to remain stable for 4.6 Gyr, as estimated here, whereas the dashed curves delineate the (much smaller) region of the plane for which the system is Hill stable. The allowed region is that below each curve. Figure 8 emphasizes that the requirement of system stability over geological, biological, and even cosmological time scales (equations ) can be less restrictive and more relevant than the requirement of Hill stability (equation ).

## 5 Conclusion

This paper presents numerical simulations of Earth-like planets in binary systems. These binaries have solar mass primaries and companion masses in the range 0.001 $`M_{}`$ (Jovian mass) to 0.5 $`M_{}`$ (K stars), although the results can be scaled to other choices. The first result of this work is an estimate of the ejection times for Earth-like planets in these systems (Figs. 3 – 6). For a given companion mass, the ejection time depends most sensitively on the periastron distance $`R_{\mathrm{min}}`$. Although the ejection time shows a wide range for a given value of $`R_{\mathrm{min}}`$, the overall trend is well-defined and the expectation value of the distribution can be described by an exponential function of $`R_{\mathrm{min}}`$. We have fit our numerical results, for each companion mass, to find the mean ejection time as a function of periastron (see equation  and Table 1).

The second result of this work is an estimate for the fraction of binary systems with Sun-like primaries that allow an Earth-like planet to remain stable for a specified time period. Our numerical experiments suggest that the requirement for stability can be written in terms of a minimum periastron distance for the binary orbit. The resulting constraint is provided by equations  (where the fitting parameters $`\alpha `$ and $`\tau _{\mathrm{ej}}^{}{}_{0}{}^{}`$ have been calculated for companion masses $`M_C`$ = 0.001 – 0.5 $`M_{}`$). For the observed distributions of binary orbital parameters (DM91), the fraction of binaries that have periastron distances greater than a given value is specified by equations  and Figure 7. Taken together, these results imply that at least 50 percent of the binary systems allow an Earth-like planet to remain stable for 4.6 Gyr (the current age of the solar system). Because the ejection time is a much more sensitive function of periastron than the fraction of binary systems (with periastron greater than a given value), this estimate is quite robust; for example, if we were to adopt the overly conservative stability requirement that $`R_{\mathrm{min}}>20`$ AU, the fraction of viable binary systems would still be 40 percent. We also note that the condition of system stability over billions of years is much less restrictive than the requirement of Hill stability (see Figure 8).

This type of stability analysis can be applied to the observed planetary systems now being discovered in association with nearby stars (e.g., Mayor & Queloz 1995; Butler et al. 1999, Marcy et al. 2001, Fischer et al. 2002). For a subset of the known extrasolar planetary systems, those with $`a>1.5`$ AU, we have performed additional sets of numerical simulations using the masses and orbital properties of the observed giant planets as the companions. An Earth-like planet, with the mass of Earth and an orbital radius of 1 AU, is assumed to reside in each system; we then study its intermediate term prospects for stability. Figure 9 shows the result – the estimated ejection time for Earth-like planets plotted as a function of the observed periastron distance of the giant planet. Each system displays a range of ejection times for the hypothetical Earth (the vertical bars in the Figure show the standard deviation of this distribution). The systems with periastron greater than about 1.8 AU allow an Earth-like planet to survive for more than 10 Myr. Extrapolating this result to the age of the galaxy, about 10 Gyr, we estimate that systems with Jupiter mass companions and periastron greater than 2.8 AU would remain stable; for companions with 10 Jupiter masses, the minimum periastron distance is about 4 AU.

The parameter space available to this class of systems is enormous and additional numerical work should be carried out. This paper provides an exploration of the $`aϵ`$ plane using integration times out to 10 Myr, but longer term simulations should be done for $`t10`$ Myr. In the regime studied here, the ejection time varies according to the exponential law of equation . At sufficiently large value of periastron, however, chaotic motion should no longer occur and the system should become stable. In addition, this work focuses on co-planar orbits, although longer ejection times can be realized for varying orbital inclinations. This work also studies only single planets, whereas multiple planets systems can also be considered. Under favorable circumstances, multiple planets can protect each other from ejection (M. Holman, private communication); if the orbital precession induced by the other planets occurs on a shorter time scale than that induced by the binary companion, the perturbations of the companion can be washed out. This work only considers the stability of Earth-like planets with outer binaries, where the Earth lies within the orbit of the secondary; the case of inner binaries, with Earth orbiting the binary pair, should also be investigated. Finally, this paper studies the stability of three-body systems after the Earth-like planet has formed. The presence of a binary companion can affect the planetary formation process, and an investigation of this issue is underway (Quintana 2003).

Acknowledgements

This study began as an REU project for EMD at the University of Michigan. We would like to thank Matt Holman and Greg Laughlin for useful discussions and an anonymous referee for many useful comments that improved the paper. EMD and MF are supported by the Hauck Foundation through Xavier University. EVQ and FCA are supported by NASA through a grant from the Origins of Solar Systems Program and by the University of Michigan through the Michigan Center for Theoretical Physics.


text:
# Loss of superfluidity in a Bose-Einstein condensate via forced resonant oscillations

111E-mail address: adhikari@ift.unesp.br

## Abstract

We predict the loss of superfluidity in a Bose-Einstein condensate in an axially symmetric harmonic trap alone during resonant collective oscillations via a classical dynamical transition. The forced resonant oscillation can be initiated by (a) a periodic modulation of the atomic scattering length with a frequency that equals twice the radial trapping frequency or multiples thereof, or by (b) a periodic modulation of the radial trapping potential with a frequency that equals the radial trapping frequency or multiples thereof. Suggestion for future experiment is made.

The first experimental observation of quantum phase effects on a macroscopic scale such as interference of matter waves  was made in a Bose-Einstein condensate (BEC) in an axially-symmetric harmonic trap. Later on more controlled studies of interference of matter waves were performed with BEC loaded on a harmonic plus an one-dimensional optical lattice trap . More recently, a three-dimensional optical lattice trap has been employed by Greiner et al.  in the study of the formation of an interference pattern. The formation of the interference pattern is a consequence of phase coherence in the BEC across the optical lattice sites.

In the experiments with optical lattice, phase coherence in the BEC is generated by quantum tunneling of atoms from one optical lattice site to another leading to a communication among various sites. As the strength of the optical lattice barrier is larger than the typical energy of the system by about two orders of magnitude, classical movement of uncondensed cold atoms is prohibited across the optical lattice barriers. It is the atoms of the BEC which experience this miraculous flow through the high barriers, which is a manifestation of superfluidity. More surprising is the observation that all superfluids are necessarily condensates, but all condensates are not superfluids. There have been both theoretical and experimental studies of breakdown of superfluidity in BEC via quantum  and classical  transitions.

As superfluidity is necessarily a quantum phenomenon the experimental consideration of Greiner et al.  on the loss of superfluidity in a BEC trapped on a three-dimensional optical lattice potential via a quantum phase transition is worth mentioning. The long-range phase coherence in the condensate along the entire optical lattice is a sign of communication among various sites which is necessary for developing superfluidity in the condensate. Equal phase at all points or a slowly (and orderly) varying phase are the ideal examples of coherent phase. On the other hand, a rapidly (and randomly) varying phase in space is generally incoherent. The phase on an optical lattice site and the number of atoms in that site play the roles of conjugate variables obeying the Heisenberg uncertainty principle of quantum mechanics . In the superfluid state the coherent phase is considered to be known and consequently the number of atoms on each site is unknown thus allowing a free movement of atoms from one site to another . The loss of phase coherence is associated with a Mott insulator state where the phase is entirely arbitrary across the optical lattice sites and the number of atoms at each site is fixed and consequently, their free passage from one site to another is stopped.

Greiner et al.  demonstrated that, as the strength of the optical lattice traps is increased, the quantum tunneling of atoms from one optical site to another as well as the communication among different optical lattice sites are stopped resulting in a superfluid to Mott insulator quantum phase transition in the BEC . As the strength of the optical lattice traps in the Mott insulator state is reduced the superfluidity is restored in a short time . This reversible quantum phase transition may occur at absolute zero (0 K) and is driven by Heisenberg’s uncertainty principle  and not by thermal fluctuations involving energy as in a classical phase transition. As the temperature approaches absolute zero all thermal fluctuations die out and at 0 K classical phase transitions are necessarily excluded.

Following a suggestion by Smerzi et al. , Cataliotti et al.  demonstrated in a novel experiment the loss of superfluidity in a BEC trapped in an one-dimensional optical-lattice and harmonic potentials when the center of the harmonic potential is suddenly displaced along the optical lattice through a distance larger than a critical value. Then a modulational instability of classical nature takes place in the BEC. Consequently, it cannot reorganize itself quickly enough via phase coherent states and the superfluidity of the BEC is lost. Later on, it has been suggested that superfluidity could be lost in a resonant collective oscillation of a BEC on an one-dimensional optical lattice potential arising from a periodic modulation of the scattering length  or a periodic modulation of the radial trapping potential .

All the above studies on the destruction of superfluidity in a BEC were performed with an optical lattice trap and the loss of superfluidity in both classical and quantum cases has been traced to the fixing of a specific atom(s) at a definite lattice site with no mobility . An interesting question is if, in case of the loss of superfluidity via the classical transitions above , the forced tunneling of the atoms through the optical lattice barriers plays a fundamental role. We find that the answer to this question is negative and demonstrate that the superfluidity can also be destroyed in a classical resonant oscillation of a BEC in a harmonic trap alone without any accompanying optical lattice trap. In view of comments in the literature  this is surprising. Effectively, the optical lattice is found to have no effect on the loss of mobility of the atoms associated with quantum fluctuations due to Heisenberg’s uncertainty principle.

In the present study the forced classical resonant oscillation of the BEC is initiated by a periodic modulation of the scattering length or the radial trapping potential. A periodic modulation of the scattering length  or the radial trapping potential  is known to lead to collective resonant oscillation in an axially-symmetric BEC. In previous studies we demonstrated that such forced oscillations in the presence of a joint harmonic and optical traps lead to a breakdown of superfluidity . A similar breakdown of superfluidity of a BEC in a harmonic trap alone reveals that the optical lattice trap does not play a decisive role in the loss of superfluidity via collective resonant oscillation. The optical lattice trap is, however, fundamental in the formation of the interference pattern upon release from the traps, which plays a decisive role in the detection of superfluidity in these studies. The strong optical lattice potential essentially cuts the BEC into phase coherent pieces upon release from the joint traps which is fundamental in the formation of the interference pattern.

Another way of detection of phase coherence and superfluidity of a BEC is by cutting it into phase coherent pieces using a laser as in the experiment by Andrews et al. . Upon free expansion these pieces form an interfence pattern. If the BEC is cut into two pieces the interference pattern consists of a large number of dark (absence of matter) and bright (matter) patches . However, we show that if the initial BEC is cut into many equal coherent pieces as in the experiment with optical lattice, only three prominent bright spots can be generated. In the present numerical simulation with an axially symmetric harmonic trap, the initial wave function $`\psi `$ of the BEC is cut into many equal pieces by setting $`\psi (\rho ,z=\pm jd)=0`$ with $`j`$ = 0, 1, 2, 3,…etc. Here $`d`$ is the spacing in the axial $`z`$ direction and $`\rho `$ corresponds to the radial direction. Subsequently, for a sufficiently small $`d`$ ($`0.11\mu `$m), upon release from the trap such a coherent fragmented BEC will form a interference pattern of three bright spots.

We consider here two ways of initiating the collective resonant oscillation of the BEC. First, it is initiated near a Feshbach resonance  by a periodic modulation of the repulsive atomic scattering length $`a`$ $`(>0)`$ via $`aa+\overline{A}\mathrm{sin}(\mathrm{\Omega }t)`$ where $`t`$ is time, $`\overline{A}`$ an amplitude, and $`\mathrm{\Omega }`$ the frequency of modulation. Such modulation of the scattering length can be realized experimentally near a Feshbach resonance by manipulating an external background magnetic field. When $`\mathrm{\Omega }=2\omega `$ or multiples thereof, resonant collective oscillation can be generated in the BEC, where $`\omega `$ is the radial trapping frequency . This resonant oscillation destroys the superfluidity of the BEC provided that the condensate is allowed to experience this oscillation for a certain interval of time called hold time.

Next we generate the resonant collective oscillation by a periodic modulation of the radial trapping potential $`V_\rho `$ in the axially symmetric BEC via $`V_\rho V_\rho [1+\overline{B}\mathrm{sin}(\mathrm{\Omega }t)]`$ with $`\overline{B}`$ the amplitude of modulation. When $`\mathrm{\Omega }=\omega `$ or multiples thereof, resonant collective oscillation can be generated in the BEC . This phenomenon has also been explored experimentally . We find that this resonant oscillation also destroys superfluidity of the BEC after some hold time.

The transition we consider is classical in nature and can be treated  by the mean-field nonlinear Gross-Pitaevskii (GP) equation . The time-dependent axially symmetric BEC wave function $`\mathrm{\Psi }(𝐫;t)`$ at position $`𝐫`$ and time $`t`$ is described by the following GP equation

$`[i\mathrm{}{\displaystyle \frac{}{t}}{\displaystyle \frac{\mathrm{}^2^2}{2m}}+V(𝐫)+gN|\mathrm{\Psi }(𝐫;t)|^2]\mathrm{\Psi }(𝐫;t)=0,`$ (1)
where $`m`$ is the mass and $`N`$ the number of atoms in the condensate, $`g=4\pi \mathrm{}^2a/m`$ the strength of interatomic interaction, with $`a`$ the atomic scattering length, and $`V(𝐫)=\frac{1}{2}m\omega ^2(\rho ^2+\nu ^2z^2)`$ the trapping potential where $`\omega `$ is the angular frequency of the harmonic trap in the radial direction $`\rho `$, $`\nu \omega `$ that in the axial direction $`z`$, with $`\nu `$ the aspect ratio. The normalization condition is $`𝑑𝐫|\mathrm{\Psi }(𝐫;t)|^2=1.`$

In the axially-symmetric configuration, the wave function can be written as $`\mathrm{\Psi }(𝐫,t)=\psi (\rho ,z,t)`$. Now transforming to dimensionless variables $`\widehat{\rho }=\sqrt{2}\rho /l`$, $`\widehat{z}=\sqrt{2}z/l`$, $`\tau =t\omega ,`$ $`l\sqrt{\mathrm{}/(m\omega )}`$, and $`\phi (\widehat{\rho },\widehat{z};\tau )\widehat{\rho }\sqrt{l^3/\sqrt{8}}\psi (\rho ,z;t),`$ Eq. (1) becomes

$`[`$ $`i{\displaystyle \frac{}{\tau }}{\displaystyle \frac{^2}{\widehat{\rho }^2}}+{\displaystyle \frac{1}{\widehat{\rho }}}{\displaystyle \frac{}{\widehat{\rho }}}{\displaystyle \frac{^2}{\widehat{z}^2}}+{\displaystyle \frac{1}{4}}\left(\widehat{\rho }^2+\nu ^2\widehat{z}^2\right)`$ (2)
$`+`$ $`{\displaystyle \frac{V_{\text{opt}}}{\mathrm{}\omega }}{\displaystyle \frac{1}{\widehat{\rho }^2}}+8\sqrt{2}\pi n\left|{\displaystyle \frac{\phi (\widehat{\rho },\widehat{z};\tau )}{\widehat{\rho }}}|^2\right]\phi (\widehat{\rho },\widehat{z};\tau )=0,`$
where nonlinearity $`n=Na/l`$. In terms of the one-dimensional probability $`P(z,t)2\pi `$ $`_0^{\mathrm{}}𝑑\widehat{\rho }|\phi (\widehat{\rho },\widehat{z},\tau )|^2/\widehat{\rho }`$, the normalization of the wave function is given by $`_{\mathrm{}}^{\mathrm{}}𝑑\widehat{z}P(z,t)=1.`$ The probability $`P(z,t)`$ is useful in the study of the present problem under the action of the optical lattice, specially in the investigation of the formation and evolution of the interference pattern after the removal of the trapping potentials.

We use the parameters of the experiment of Cataliotti et al.  with repulsive 87Rb atoms where the radial trap frequency was $`\omega =2\pi \times 92`$ Hz. For the mass $`m=1.441\times 10^{25}`$ kg of 87Rb the harmonic oscillator length $`l=\sqrt{\mathrm{}/(m\omega )}=1.126`$ $`\mu `$m and and the dimensionless time unit $`\omega ^1=1/(2\pi \times 92)`$ s $`=1.73`$ ms. We solve Eq. (2) numerically using a split-step time-iteration method with the Crank-Nicholson discretization scheme described recently . The time iteration is started with the known harmonic oscillator solution of Eq. (2) with $`n=0`$: $`\phi (\widehat{\rho },\widehat{z})=[\nu /(8\pi ^3)]^{1/4}`$ $`\widehat{\rho }e^{(\widehat{\rho }^2+\widehat{z}^2)/4}`$ . First we calculate the ground-state wave function of the system for $`n=5`$ and $`\nu =0.5`$ which we use in numerical simulation. To investigate the superfluidity of the BEC we set $`\phi (\widehat{\rho },\widehat{z})=0`$ for $`\widehat{z}=\pm 0.3j`$, $`j=0,1,2,3,\mathrm{}`$. This cuts the initial BEC in slices of width 0.3 in the axial direction. The free expansion of a phase coherent initial BEC sliced in this fashion will lead to a prominent interference pattern revealing the superfluidity. Essentially, similar slices result in a BEC trapped in an one-dimensional optical lattice as in the experiment by Cataliotti et al . The atom cloud released from one piece expand, and overlap and interfere with atom clouds from neighboring piece to form the robust interference pattern due to phase coherence. The pattern consists of a central peak and two symmetrically spaced peaks moving in opposite directions . Although we employ the dimensionless space units $`\widehat{\rho }`$ and $`\widehat{z}`$ and time unit $`\tau `$ in numerical calculation, the results are reported in actual units: $`r`$ $`\mu `$m, $`\rho `$ $`\mu `$m, $`z`$ $`\mu `$m and $`t`$ ms. In the conversion we used the parameters of the experiment of Cataliotti et al. , e.g., $`\rho =0.8\widehat{\rho }`$ $`\mu `$m, $`z=0.8\widehat{z}`$ $`\mu `$m, and $`t=1.73\tau `$ ms.

First we study the destruction of superfluidity in the condensate after the application of a periodic modulation of the scattering length resulting in a similar modulation of nonlinearity $`n`$ in Eq. (2) via

$$nn+A\mathrm{sin}(\mathrm{\Omega }\tau ),$$
(3)
where $`A`$ is an amplitude. In the present study we employ amplitude $`A=3`$. For present $`n=5`$, this will restrict the modulated nonlinearity to have always positive values corresponding to atomic repulsion. Negative values of nonlinearity corresponding to atomic attraction may lead to collapse and instability  and will not be considered here. The BEC has been found to execute resonant collective oscillation when the modulation frequency $`\mathrm{\Omega }`$ is an even multiple of the radial trap frequency $`\omega `$ . In the presence of an one-dimensional optical lattice trap such resonant nonlinear oscillation destroys superfluidity . Here we investigate if the forced quantum tunneling of the atoms across the optical lattice trap during classical oscillation of the BEC is responsible for the loss of superfluidity.

We explicitly study the destruction of superfluidity in the condensate upon the application of modulation (3) of the scattering length leading to a resonant oscillation at $`\mathrm{\Omega }=2\omega `$. The loss of superfluidity only takes place if the BEC is allowed to experience the resonant oscillation for a certain hold time. For numerical simulation we allow the BEC to evolve on a lattice with $`\rho 20`$ $`\mu `$m and 20 $`\mu `$m $`z`$ $`20`$ $`\mu `$m after the modulation (3) is applied and study the system after different hold times. The one-dimensional probability $`P(z,t)`$ is plotted in Figs. 1 (a), (b), (c) and (d) for hold times 0, 17 ms, 35 ms and 69 ms, respectively. For hold time 17 ms, prominent interference pattern is formed upon free expansion. In Fig. 1 (a) three separate piece in the interference pattern corresponding to three distinct trails can be identified. The interference pattern is slowly destroyed at increased hold times as we can see in Figs. 1 (b), (c) and (d). As the hold time increases the maxima of the interference pattern mix up upon free expansion and finally for the hold time of 69 ms the interference pattern is completely destroyed as we find in Fig. 1 (d). As the BEC is allowed to evolve for a substantial interval of time after the application of the periodic modulation in the scattering length, a dynamical instability of classical nature sets in which destroys the superfluidity .

Next we investigate how the phase $`\delta `$ over the condensate behaves as the hold time is increased and what happens to the modulus $`|\psi (\rho ,z)|`$ of the wave function. For this purpose we plot the phase and the modulus of the wave function of the BEC’s of Figs. 1 (a), (b), (c), and (d), respectively, in Figs. 2 (a), (b), (c), and (d) after hold times 0, 17 ms, 35 ms, and 69 ms. For hold time 0 the wave function is smooth and the phase is slowly varying over the full condensate corresponding to perfect phase coherence. But as the hold time is increased the phase varies rapidly over the condensate. The total variation of phase over the BEC for hold times 0, 17 ms, 35 ms, and 69 ms are 1.5 rad, 15 rad, 20 rad, and 45 rad, respectively. For larger hold times the modulus $`|\psi (\rho ,z)|`$ of the wave function also develops a nonsmooth irregular behavior via the dynamical instability. The rapid variation of $`\delta `$ and the nonsmooth irregular wave function over the BEC for larger hold times are responsible for the destruction of superfluidity. Hence the superfluidity in the BEC could be destroyed in the absence of an optical trap due a classical collective resonant oscillation in the condensate initiated by a periodic modulation of the scattering length with a frequency $`\mathrm{\Omega }=2\omega `$. Such resonances appear for $`\mathrm{\Omega }=2\omega j`$ , where $`j=1,2,3,\mathrm{}`$ and we have verified that similar breakdown of superfluidity also occurs for $`\mathrm{\Omega }=4\omega `$. The resonance becomes more vigorous as $`A`$ is increased and so is the destruction of superfluidity.

However, resonant oscillation arising from modulation (3) is not the only classical dynamical mechanism for the destruction of superfluidity. It can also be destroyed via a periodic modulation of the radial trapping potential . To illustrate this we consider the following modulation of the radial trapping potential in (2) via

$$\widehat{\rho }^2/4(\widehat{\rho }^2/4)[1+B\mathrm{sin}(\mathrm{\Omega }\tau )],$$
(4)
while the axial trapping potential is left unchanged. In this case parametric resonance appear when the modulation frequency $`\mathrm{\Omega }`$ is a multiple of the trap frequency $`\omega `$ . These resonances appear even for a small value of the amplitude $`B`$. However, they appear for a band of frequency around the multiples of $`\omega `$ for a larger value of $`B`$.

We find that the superfluidity is destroyed if the modulation (4) of the radial trapping potential is continued for a certain hold time with $`\mathrm{\Omega }`$ fixed at a resonant frequency. To illustrate the loss of superfluidity in this case we plot the one-dimensional probability distribution $`P(z,t)`$ vs. $`z`$ and $`t`$ in Figs. 3 (a), (b), and (c) for hold times 35 ms, 69 ms, and 104 ms, respectively, for $`\mathrm{\Omega }=\omega `$ and $`B=0.5`$. In Fig. 3 (a) we find that the superfluidity is maintained for a hold time of 35 ms and three prominent peaks are formed in this case. However, as the hold time is increased the clean interference pattern is gradually destroyed as can be seen from Figs. 3 (b) and (c). The same phenomenon is found to occur at a larger resonant frequency of modulation, where $`\mathrm{\Omega }`$ is a larger multiple of $`\omega `$. The destruction of superfluidity is favored for a larger amplitude of modulation $`B`$.

If the frequencies $`\mathrm{\Omega }`$ are off their resonance values and the amplitudes $`A`$ and $`B`$ in Eqs. (3) and (4), respectively, have moderate values, there is no loss of superfluidity after large hold times for modulations of scattering length or radial trapping potential. To illustrate this for modulation (4) with $`B=0.5`$, $`\mathrm{\Omega }=1.2\omega `$ and a hold time of 104 ms we plot the density $`P(z,t)`$ vs. $`z`$ and $`t`$ in Fig. 4, where the clean interference pattern reappears with a small change in the modulation frequency $`\mathrm{\Omega }`$ from the resonance value $`\omega `$ to a nonresonant value $`1.2\omega `$. The maintenance of superfluidity in Fig. 4 off the resonance should be contrasted with its loss in Fig. 3 (c) at resonance for the same hold time.

In conclusion, we studied the destruction of superfluidity in a cigar-shaped trapped BEC upon the application of a periodic modulation of the scattering length or of the radial trapping potential leading to a resonant collective excitation. In the absence of modulation, the formation of the interference pattern upon the removal of the trap clearly demonstrates the phase coherence. At the resonance frequencies a dynamical instability in the BEC leads to the destruction of superfluidity provided that the BEC is kept in the modulated trap for a certain hold time. Consequently, after release from the trap no interference pattern is formed. The superfluidity in the BEC is not destroyed when the frequency of modulation is changed to a nonresonant value. It is possible to study this novel way of the destruction of superfluidity experimentally and a comparison of those results with mean-field models will enhance our understanding of matter wave BEC.

The work was supported in part by the CNPq and FAPESP of Brazil.