text:
# Homoclinic Chaos in the Dynamics of a General Bianchi IX Model

## Abstract

The dynamics of a general Bianchi IX model with three scale factors is examined. The matter content of the model is assumed to be comoving dust plus a positive cosmological constant. The model presents a critical point of saddle-center-center type in the finite region of phase space. This critical point engenders in the phase space dynamics the topology of stable and unstable four dimensional tubes $`R\times S^3`$, where $`R`$ is a saddle direction and $`S^3`$ is the manifold of unstable periodic orbits in the center-center sector. A general characteristic of the dynamical flow is an oscillatory mode about orbits of an invariant plane of the dynamics which contains the critical point and a Friedmann-Robertson-Walker (FRW) singularity. We show that a pair of tubes (one stable, one unstable) emerging from the neighborhood of the critical point towards the FRW singularity have homoclinic transversal crossings. The homoclinic intersection manifold has topology $`R\times S^2`$ and is constituted of homoclinic orbits which are bi-asymptotic to the $`S^3`$ center-center manifold. This is an invariant signature of chaos in the model, and produces chaotic sets in phase space. The model also presents an asymptotic DeSitter attractor at infinity and initial conditions sets are shown to have fractal basin boundaries connected to the escape into the DeSitter configuration (escape into inflation), characterizing the critical point as a chaotic scatterer.

The longtime debate on the chaotic dynamics of general Bianchi IX models started with the work of Belinskii, Khalatnikov and Lifshitz (BKL) on the oscillatory behaviour of such models in their approach to the singularityBKL . They showed that the approach to the singularity($`t0`$) of a general Bianchi IX cosmological solution is an oscillatory mode, consisting of an infinite sequence of periods (called Kasner eras) during which two of the scale functions oscillate and the third one decreases monotonically; on passing from one era to another the monotonic behaviour is transfered to another of the three scale functions. The length of each era is determined by a sequence of numbers $`X_s,0<X_s<1,s=integer`$, each of which arises from the preceding one by the transformation $`X_{s+1}`$ = fractional part of $`1/X_s`$. From the properties of this map it is obtained that the behaviour of the model becomes chaotic on approaching the singularity($`t=0`$) for arbitrary initial conditions given at $`t_0>0`$. With the advent of powerful numerical resources the interest in the behaviour of these models revived but - as in the BKL work - the procedure has been basically to obtain maps which approximate the dynamics of the model described by Einstein’s equations and which exhibit strong stochastic propertiesbarrow . How well these discrete maps represent the full nonlinear dynamics has been subject of much research, particularly by Bergerberger1 ; berger2  and Rughrugh . The interest in the chaoticity of Bianchi IX models has been mainly focused on the Mixmaster case (vacuum Bianchi IX model with three scale factorsmisner ), but the chaotic dynamics of other Bianchi IX model universes has also been discussed in the literature (cf. burd  and references therein, and cornish ); homoclinic chaos in axisymmetric Bianchi IX universes with matter and cosmological constant has been treated in barguine . The chaoticity in the Mixmaster dynamics has been object of much dispute in the literature (cf. the contributions to the Section Bianchi IX (Mixmaster) dynamics in burd1 , particularly berger2 ). Latifi et al.latifi  and Contopoulos et al.contopoulos  have shown the nonintegrability of the Mixmaster model in the Painlevé sense, although the question of the generic behaviour (chaotic or not) remained unsettled mainly due to the absence of an invariant (or topological) characterization of chaos in the model (standard chaotic indicators as Liapunov exponents being coordinate dependent and therefore questionablematsas burd1 ). More recently Cornish and Levinlevin  proposed to quantify chaos in the Mixmaster universe by calculating the dimensions of fractal basin boundaries in initial conditions sets for the full dynamics, these boundaries being defined by the code associated with one of the three outcomes on which one of the three axes is collapsing most quickly, as established numerically. Their result received a recent critical review that nevertheless endorses itletelier . In the present paper we shall give an invariant characterization of chaos for the Mixmaster model with a positive cosmological constant, but this criterion does not work for the case of zero cosmological constant.

Our purpose in this paper is to examine the dynamics of a general (three scale factors) Bianchi IX cosmological model with dust and cosmological constant and establish the existence of chaos using the criterion of the homoclinic transversal crossing of topological tubes that organizes the dynamical flow in the model. We show that the dynamics of the model is highly complex and chaotic, and that chaos has a definite homoclinic origin. The phase space of the model is noncompact and the presence of the cosmological constant determines two crucial facts in phase space: first, the existence of a critical point of the type saddle-center-center; second, two critical points at infinity corresponding to the DeSitter configuration, one acting as an “attractor” to the dynamics. With respect to the latter point, our system has the characteristic of a chaotic scattering system with two absolute outcomes consisting of (i) escape into inflation (the DeSitter attractor) or (ii) recollapse to the singularity. The presence of the critical point of saddle-center-center type is responsible for a wealthy and complex dynamics, engendering in phase space topological structures which are homoclinic to a manifold of periodic orbits, having the topology of the three sphere $`S^3`$. The dynamics has a 2-dim invariant plane allowing to expand the dynamical equations about it, in the same procedure used to examine the dynamics about a periodic orbit of the system. The mixmaster universe with a cosmological constant will correspond to the limiting case of dust equal to zero but, as we will discuss, the skeleton of its dynamics is framed in the dynamics of the general case.

The line element of the model has the form

$`ds^2`$ $`=`$ $`dt^2(\theta ^1)^2(\theta ^2)^2(\theta ^3)^2,`$ (1)
where $`t`$ is the cosmological time and $`\theta ^1=M(t)\omega ^1`$, $`\theta ^2=N(t)\omega ^2`$ and $`\theta ^3=R(t)\omega ^3`$. Here the $`\omega ^i`$ are Bianchi IX 1-forms satisfying $`d\omega ^i=ϵ^{ijk}\omega ^j\omega ^k`$. The matter content of the models is assumed to be a pressureless perfect fluid, namely dust, with energy density $`\rho `$, as described by the comoving observers with four velocity orthogonal to the homogeneity surfaces of (1), plus a positive cosmological constant $`\mathrm{\Lambda }`$. The dynamics of the three scale factors M(t), N(t) and R(t) will be given by Einstein’s equations, which are equivalent to Hamilton’s equations generated by the Hamiltonian constraint

$`H`$ $`=`$ $`{\displaystyle \frac{1}{8}}({\displaystyle \frac{M}{NR}}p_M^2{\displaystyle \frac{N}{MR}}p_N^2`$ (2)
$``$ $`{\displaystyle \frac{R}{MN}}p_R^2+{\displaystyle \frac{2}{R}}p_Mp_N+{\displaystyle \frac{2}{M}}p_Np_R+{\displaystyle \frac{2}{N}}p_Mp_R)`$
$`+`$ $`{\displaystyle \frac{1}{2MNR}}(R^4+M^4+N^4(R^2N^2)^2(R^2M^2)^2`$
$``$ $`(M^2N^2)^2)2\mathrm{\Lambda }MNR2E_0=0.`$
Here $`p_M`$, $`p_N`$ and $`p_R`$ are the momenta canonically conjugate to $`M`$, $`N`$ and $`R`$, respectively. $`E_0`$ is a constant proportional to the total energy of the model, arising from the first integral of the Bianchi identities, $`\rho MNR=E_0`$. Hamilton’s equations have the form

$`\dot{M}`$ $`=`$ $`{\displaystyle \frac{H}{p_M}}=\left({\displaystyle \frac{M}{4NR}}p_M+{\displaystyle \frac{1}{4R}}p_N+{\displaystyle \frac{1}{4N}}p_R\right),`$
$`\dot{p}_M`$ $`=`$ $`{\displaystyle \frac{H}{M}}=({\displaystyle \frac{1}{8NR}}p_M^2{\displaystyle \frac{N}{8M^2R}}p_N^2`$
$``$ $`{\displaystyle \frac{R}{8M^2N}}p_R^2+{\displaystyle \frac{2}{8M^2}}p_Np_R)`$
$`+`$ $`{\displaystyle \frac{1}{2M^2NR}}(R^4+M^4+N^4(R^2N^2)^2`$
$``$ $`(R^2M^2)^2(M^2N^2)^2)+2\mathrm{\Lambda }NR`$
$``$ $`{\displaystyle \frac{1}{2MNR}}(4M^3+4(R^2M^2)M`$
$``$ $`4(M^2N^2)M),`$
and cyclically in $`M`$, $`N`$, $`R`$, $`p_M`$,$`p_N`$, $`p_R`$. The dynamical system (Homoclinic Chaos in the Dynamics of a General Bianchi IX Model) and (Homoclinic Chaos in the Dynamics of a General Bianchi IX Model) presents one critical point E with coordinates

$$p_M=p_N=p_R=0,M=N=R=M_0,$$
(5)
where $`M_0=\sqrt{\frac{1}{4\mathrm{\Lambda }}}`$. The energy of the critical point E is defined by

$`2E_{crit}=\sqrt{{\displaystyle \frac{1}{4\mathrm{\Lambda }}}}.`$ (6)
Much of our understanding of nonlinear systems derives from linearization about critical points and the determination of existing invariant submanifolds, which are structures that actually organize the dynamics in phase space. The system presents a two dimensional invariant manifold of the dynamics, defined by

$$p_M=p_N=p_R,M=N=R.$$
(7)
This invariant plane is actually the intersection of the two four dimensional invariant submanifolds defined by ($`M=N,p_M=p_N`$) and ($`N=R,p_N=p_R`$). The saddle-center-center critical point $`E`$ belongs to the invariant plane, and from it emerge the separatrices $`S`$. The FRW singularity ($`M=N=R=0`$,$`p_M=p_N=p_R=0`$) is a degenerate critical point of the dynamics on the invariant plane. The phase picture of the motion in the invariant plane is given in Fig. 1, already expressed in canonical coordinates $`(x,p_x)`$, to be introduced in Eqs. (21)-(Homoclinic Chaos in the Dynamics of a General Bianchi IX Model), that are variables defined on the invariant plane. Also, a straightforward analysis of the infinity of the phase space under consideration shows that it has two critical points in this region, corresponding to the DeSitter solution, one acting as an attractor (stable DeSitter configuration) and the other as a repeller (unstable deSitter configuration) for the dynamics at infinity. The scale factors $`M`$, $`N`$ and $`R`$ approach the DeSitter attractor as $`M=N=Rexp[t\sqrt{\mathrm{\Lambda }/3}]`$ and $`p_M=p_N=p_Rexp[2t\sqrt{\mathrm{\Lambda }/3}]`$. It is easy to see that the DeSitter asymptotic configurations also belong to the invariant plane and that two of the separatrices $`S`$ approach them for times going to $`\pm \mathrm{}`$.

To proceed with the study of the phase space dynamics, let us now linearize the dynamical equations (Homoclinic Chaos in the Dynamics of a General Bianchi IX Model), (Homoclinic Chaos in the Dynamics of a General Bianchi IX Model) about the critical point E. We define

$`X`$ $`=`$ $`(MM_0),W=(p_M0),`$
$`Y`$ $`=`$ $`(NM_0),K=(p_N0),`$ (8)
$`Z`$ $`=`$ $`(RR_0),L=(p_R0).`$
We obtain

$$\left(\begin{array}{c}\dot{X}\\ \dot{Y}\\ \dot{Z}\\ \dot{W}\\ \dot{K}\\ \dot{L}\end{array}\right)=\frac{1}{M_0}S\left(\begin{array}{c}X\\ Y\\ Z\\ W\\ K\\ L\end{array}\right)$$

where

$$\left(\begin{array}{cccccc}0& 0& 0& 1/4& 1/4& 1/4\\ 0& 0& 0& 1/4& 1/4& 1/4\\ 0& 0& 0& 1/4& 1/4& 1/4\\ 3& 1& 1& 0& 0& 0\\ 1& 3& 1& 0& 0& 0\\ 1& 1& 3& 0& 0& 0\end{array}\right)$$

In what follows, without loss of generality, we fix $`\mathrm{\Lambda }`$ such that $`M_0=1`$. The nature of the critical point E is determined by the characteristic polynomial associated with the linearization matrix $`S`$. We obtain

$`P(\lambda )=(\lambda ^2{\displaystyle \frac{1}{4}})(\lambda ^2+2)^2,`$ (9)
with roots

$`\lambda =\pm {\displaystyle \frac{1}{2}},\lambda =\pm i\sqrt{2},`$ (10)
where the second pair has multiplicity two. The pair of real eigenvalues generates a saddle structure, while the degenerate pair of imaginary eigenvalues generate a double center structure. The analysis of the center structure will reveal a manifold of linearized unstable periodic orbits with the topology of a $`3`$-sphere, as we will see.

To display the structure of the linearized motion, we start by diagonalizing $`S`$ with the use of a transformation matrix $`\mathrm{}`$ whose columns are composed of six independent eigenvectors of $`S`$ siegel . A judicious choice of $`\mathrm{}`$ yields primed variables defined by the transformation

$`X^{}`$ $`=`$ $`{\displaystyle \frac{1}{3}}(X+Y+Z),`$
$`Y^{}`$ $`=`$ $`(XY),`$
$`Z^{}`$ $`=`$ $`(X+Y2Z),`$ (11)
$`W^{}`$ $`=`$ $`(W+K+L),`$
$`K^{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(WK),`$
$`L^{}`$ $`=`$ $`{\displaystyle \frac{1}{6}}(W+K2L).`$
In these new variables, the quadratic Hamiltonian about $`E`$ is expressed in the form

$`H=(E_{crit}E_0)+{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{1}{6}}W^26X^2\right)`$ (12)
$`\left({\displaystyle \frac{1}{2}}K^2+Y^2\right)\left({\displaystyle \frac{3}{2}}L^2+{\displaystyle \frac{1}{3}}Z^2\right).`$
These variables are conjugated to the pairs according to $`[X^{},W^{}]=1,[K^{},Y^{}]=1,[L^{},Z^{}]=1`$, other Poisson brackets (PB) zero. The Hamiltonian (12) is separable, and we can immediately identify the following constants for the linearized motion

$`E_{hyp}={\displaystyle \frac{1}{4}}\left({\displaystyle \frac{1}{6}}W^26X^2\right),`$ (13)
$`E_{rot_1}={\displaystyle \frac{3}{2}}L^2+{\displaystyle \frac{1}{3}}Z^2,`$ (14)
$`E_{rot_2}={\displaystyle \frac{1}{2}}K^2+Y^2,`$ (15)
$`Q_1={\displaystyle \frac{1}{3}}Y^{}Z^{}+{\displaystyle \frac{1}{2}}K^{}L^{},`$ (16)
$`Q_2=L^{}Y^{}{\displaystyle \frac{1}{3}}Z^{}K^{},`$ (17)
in the sense that they all have zero PB with the Hamiltonian (12). The first three constants appear already as separable pieces in the Hamiltonian (12). $`E_{hyp}`$ corresponds to the energy associated with motion in the saddle sector, while $`E_{rot_1}`$ and $`E_{rot_2}`$ are the rotational energies associated with the motion in the linear part of the center-center manifold, namely, on stable periodic orbits in 2-dim tori. The remaining two constants are additional symmetries that arise as consequence of the multiplicity two of the imaginary eigenvalues, and are associated to the fact that the linearized dynamics in the center-center sector is that of a two-dimensional isotropic harmonic oscillator. They are not all independent but are related by

$`4E_{rot_1}E_{rot_2}=12Q_1^2+6Q_2^2.`$ (18)
We are now ready to describe the topology of the general dynamics in a linear neighborhood of the critical point $`E`$. To describe all possible motions, the following situations must be taken into account in the linearized Hamiltonian (12). If $`E_{hyp}=0`$ two possibilities arise. First, we have $`W^{}=X^{}=0`$, implying that the motions are periodic orbits of the two-dimensional isotropic harmonic oscillator

$`H=\left({\displaystyle \frac{1}{2}}K^2+Y^2\right)+\left({\displaystyle \frac{3}{2}}L^2+{\displaystyle \frac{1}{3}}Z^2\right)`$ (19)
$`=(E_{crit}E_0).`$
By a proper canonical reescaling of the variables in (Homoclinic Chaos in the Dynamics of a General Bianchi IX Model) and (19) it is easy to see that these constant energy surfaces are hyperspheres and that the constants of motion $`Q_1`$, $`Q_2`$ and $`Q_3=(E_{rot_1}E_{rot_2})`$ satisfy the algebra of the three dimensional rotation group under the Poisson bracket operation, namely,

$`[Q_i,Q_j]=ϵ^{ijk}Q_k`$ (20)
The constant of motion $`Q_1`$ considered as a generator of infinitesimal contact transformations has a peculiar significance in characterizing the topology of the underlying group of the algebra (20)mcintosh . While $`Q_2`$ generates infinitesimal rotations of the orbits, $`Q_1`$ generates infinitesimal changes in eccentricity. The action of $`Q_1`$ is to take an orbit - let us say nearly circular - and transforms it into an orbit of higher and higher eccentricity until it collapses into a straight line. Continued application of $`Q_1`$ produces again an elliptic orbit, but now traversed in the opposite sense , so that it takes a $`720^0`$ rotation to bring the orbit back into itself. The two-valuedness of the mapping arises from the fact that the orbits are oriented. Therefore we are led to the conclusion that the group generated by these constants of motion is the unitary unimodular group and not the rotation group. It is therefore compatible to consider that the center-center manifold has indeed the topology of a three sphere $`S^3`$. Now the separate conservation of $`E_{rot_1}`$ and $`E_{rot_2}`$ (cf.Eq.(19)) allows us to show that the center-center manifold in the linear neighborhood of $`E`$ is foliated by Clifford two-dimensional surfaces in $`S^3`$ clifford , namely, two-tori $`\mathrm{}_{E_0}`$ contained in the energy surface $`E_0=const`$. Such surfaces, as well as the $`S^3`$ manifold that contains them depend continuously on the parameter $`E_0`$. We remark that these two tori will have limiting configurations which are periodic orbits, whenever $`E_{rot_1}=0`$ or $`E_{rot_2}=0`$, and correspond to the case of maximum eccentricity (for instance, a straight line in the plane $`(Y^{},Z^{}))`$.

The second possibility will be $`W^{}=\pm 6X^{}`$, that defines the linear stable $`V_S`$ and unstable $`V_U`$ manifolds of the saddle-sector. $`V_S`$ and $`V_U`$ limit regions $`I`$ ($`E_{hyp}<0`$) and regions $`II`$ ($`E_{hyp}>0`$) of motion on hyperbolae which are solutions in the separable saddle-sector (cf. (12)) $`E_{hyp}=\frac{1}{4}\left(\frac{1}{6}W^26X^2\right)`$, (cf. (12)). Note that the saddle-sector depicts the structure of Fig. 1 in the neighborhood of $`E`$, with $`V_U`$ and $`V_S`$ tangent to the separatrices at $`E`$ barguine . The direct product of $`\mathrm{}_{E_0}`$ with $`V_S`$ and $`V_U`$ generates, in the linear neighborhood of $`E`$, the structure of stable $`(\mathrm{}_{E_0}\times V_S)`$ and unstable $`(\mathrm{}_{E_0}\times V_U)`$ 3-dim tubes, which coalesce into the two-dimensional tori $`\mathrm{}_{E_0}`$ for times going to $`+\mathrm{}`$ or $`\mathrm{}`$, respectively. The energy of any orbit on these tubes is the same as that of the orbits on the two tori $`\mathrm{}_{E_0}`$. These 3-dim tubes are confined inside the 4-dim tubes which are the product of $`V_S`$ and $`V_U`$ by the $`S^3`$ center-center manifold. It is obvious that the existence of the $`S^3`$ center-center manifold is restricted to energy surfaces such that $`(E_{crit}E_0)>0`$. The intersection of the center-center manifold with the energy surface $`H=E_0`$ is in general a three sphere parametrized with $`E_0`$, from which two pairs (one stable and one unstable) of 4-dim tubes emanate (cf. Fig. 2 and caption). More particularly, in the linear neighborhood of $`E`$ (namely, for $`(E_{crit}E_0)`$ small), two pairs of 3-dim tubes emanate from the Clifford two-dimensional surface $`\mathrm{}_{E_0}`$ contained in $`S^3`$. For the case $`E_{hyp}0`$, the motion is restricted on infinite tubes resulting from the direct product of the hyperbolae, lying in the regions $`I`$ and $`II`$ of Fig. 1, with $`\mathrm{}_{E_0}`$. A general orbit which visits the neighborhood of $`E`$ belongs to the general case $`E_{hyp}0`$, $`E_{rot_1}0`$ and $`E_{rot_2}0`$. In this region the orbits have an oscillatory approach to the 3-dim linear tubes, the closer as $`E_{hyp}0`$.

We remark that the surface of the tubes constitutes a boundary for the general flow and, in the neighborhood of $`E`$, is defined by $`E_{hyp}=0`$. Depending on the sign of $`E_{hyp}`$ the motion will be confined inside the 4-dim tube (for $`E_{hyp}<0`$) and will correspond to a flow of motion distinct from the one confined outside the tube (for $`E_{hyp}>0`$). This can be seen from Eq.(12), that can be rewritten

$`E_{rot_1}+E_{rot_2}E_{hyp}=(E_{crit}E_0).`$
Let us consider, for instance, a set of initial conditions corresponding to initially expanding universes. Two distinct flows will be associated with these initial conditions, depending whether they are contained inside the tube or outside the tube. A careful examination of Fig. 2 shows us that the flow corresponding to initial conditions inside the stable tube will reach a neighborhood of $`E`$ and will return towards the FRW singularity inside the unstable tube (with $`E_{hyp}<0`$), while the flow of orbits associated to initial conditions that are outside the tube will reach the neighborhood of $`E`$ and escapes towards the DeSitter attractor along the exterior of the unstable tube that is directioned to the right. Also the flow of the stable tube on the bottom right of Fig. 2 will escape towards the DeSitter attractor along the interior of the unstable tube of the top right. We remark that the direction of the separatrices and hyperbolae in the linearized saddle-sector are the topological guide for the flow. Our main interest will rely on the stable and unstable tubes that emerge from the neighborhood of the critical point $`E`$ towards the FRW singularity and their homoclinic transversal crossings, that constitutes an invariant signature of chaos in the models.

By continuity, the nonlinear extension of the center-center manifold will maintain the topology of the three-sphere $`S^3`$, but will not be decomposable into $`E_{rot_1}`$ and $`E_{rot_2}`$ so that now only the four dimensional tubes with topology $`R\times S^3`$ are meaningful for the nonlinear dynamics.

The extension of the structure of the four dimensional tubes away from the neighborhood of $`E`$ are now to be examined, and our basic interest will reside in on the stable and unstable pair, $`S_{E_0}^3\times V_S`$ and $`S_{E_0}^3\times V_U`$, that leave the neighborhood of $`E`$ and proceed towards the Friedmann-Robertson-Walker (FRW) singularity of the model. To this end let us consider the two dimensional invariant manifold of the dynamics, defined by Eqs.(5). The two dimensional invariant plane is contained in a six dimensional phase space and it is obvious that, contrary to examples in lower dimensional systems, it does not separate the phase space in two disjoint parts. In fact the general motion about a neighborhood of the invariant plane is an oscillatory flow confined along the interior or the exterior of the four dimensional tubes $`R\times S^3`$, where the invariant plane (or more properly, one of the curves of the invariant plane) may be thought as a structure in the center of the tube, as we proceed to show. To see this, let us introduce the canonical coordinate transformation with the generating function

$`G=(MNR)^{\frac{1}{3}}p_x+{\displaystyle \frac{M}{N}}p_y+{\displaystyle \frac{MN}{R^2}}p_z`$
where $`p_x`$, $`p_y`$ and $`p_z`$ are the new momenta, resulting in

$`x=(MNR)^{\frac{1}{3}},y={\displaystyle \frac{M}{N}},z={\displaystyle \frac{MN}{R^2}}`$ (21)
and

$`p_M={\displaystyle \frac{1}{3}}{\displaystyle \frac{NR}{(MNR)^{\frac{2}{3}}}}p_x+{\displaystyle \frac{1}{N}}p_y+{\displaystyle \frac{N}{R^2}}p_z`$
$`p_N={\displaystyle \frac{1}{3}}{\displaystyle \frac{NR}{(MNR)^{\frac{2}{3}}}}p_x{\displaystyle \frac{M}{N^2}}p_y+{\displaystyle \frac{M}{R^2}}p_z`$ (22)
$`p_R={\displaystyle \frac{1}{3}}{\displaystyle \frac{NR}{(MNR)^{\frac{2}{3}}}}p_x{\displaystyle \frac{2MN}{R^3}}p_z.`$
It is worth remarking that the linearization of this canonical transformation about the critical point $`E`$ yields exactly the linear transformation (Homoclinic Chaos in the Dynamics of a General Bianchi IX Model), and that the variables $`(y,p_y,z,p_z)`$ correspond to the primed variables $`(K^{},Y^{},L^{},Z^{})`$ defined on the $`S^3`$ center-center manifold about a linear neighborhood of the critical point. The variable $`x`$ is obviously the average scale factor of the model. In these new variables, the full Hamiltonian (2) assumes the form

$`H={\displaystyle \frac{1}{24x}}p_x^2{\displaystyle \frac{y^2}{2x^3}}p_y^2{\displaystyle \frac{3z^2}{2x^3}}p_z^2`$
$`{\displaystyle \frac{x}{2z^{\frac{4}{3}}}}{\displaystyle \frac{1}{2}}xz^{\frac{2}{3}}y^2{\displaystyle \frac{1}{2y^2}}xz^{\frac{2}{3}}`$
$`+{\displaystyle \frac{x}{yz^{\frac{1}{3}}}}+{\displaystyle \frac{yx}{z^{\frac{1}{3}}}}+xz^{\frac{2}{3}}2\mathrm{\Lambda }x^3=2E_0,`$ (23)
and the equations of the invariant plane reduce to

$$y=1,z=1,p_y=0=p_z.$$
(24)
It is clear that $`(x,p_x)`$ are variables defined on the invariant plane, and now the expansion of the Hamiltonian (Homoclinic Chaos in the Dynamics of a General Bianchi IX Model) about the neighborhood of the invariant plane can be easily implemented, producing a linearized Hamiltonian parametrized by the variables $`(x(t),p_x(t))`$ describing the curves in the invariant plane, in analogy with the way we expand a dynamical system about a periodic orbit. We obtain

$`H={\displaystyle \frac{1}{24x}}p_x^2+{\displaystyle \frac{3x}{2}}2\mathrm{\Lambda }x^3{\displaystyle \frac{1}{2x^3}}p_y^2{\displaystyle \frac{3}{2x^3}}p_z^2`$
$`x(y1)^2{\displaystyle \frac{1}{3}}x(z1)^2=2E_0`$ (25)
with dynamical equations

$`(\delta y\dot{)}={\displaystyle \frac{1}{x^3}}\delta p_y`$
$`(\delta p_y\dot{)}=2x\delta y`$
$`(\delta z\dot{)}={\displaystyle \frac{3}{x^3}}\delta p_z`$ (26)
$`(\delta p_z\dot{)}={\displaystyle \frac{2x}{3}}\delta z,`$
where $`\delta y=y1`$, $`\delta z=z1`$, $`\delta p_y=p_y0`$ and $`\delta p_z=p_z0`$. The linearization matrix in (Homoclinic Chaos in the Dynamics of a General Bianchi IX Model) has imaginary eigenvalues, both with multiplicity two, given by $`\lambda =\pm \frac{i\sqrt{2}}{x(t)}`$, so that in the neighborhood of the invariant plane we have only elliptic modes, namely, the motion is oscillatory about the invariant plane with increasing frequency as the orbits approach $`x=0`$. Close to the critical point $`E`$, the motion in the invariant plane is slow, i.e., $`(\dot{x},\dot{p_x})`$ is small. Following an orbit in this plane with diminishing $`x`$, its phase space speed increases but not faster than the frequency of the surrounding motion in the other four variables. We are thus justified in deriving a qualitative picture for the motion near the invariant plane from an adiabatic approximation for the relevant orbit manifolds. That is, we will assume that the overall geometry of the four dimensional stable and unstable tubes that are asymptotic to the $`S^3`$ center-center manifold can be approximated by the adiabatic deformation of this sphere as $`x`$ decreases away from the critical point. In this way, we have in the six dimensional phase space the four dimensional tubes of motion $`R\times S^3`$ with the two dimensional invariant plane as the structure at its center, about which are the oscillatory degrees of freedom $`(y,p_y,z,p_z)`$. We will use this fact to show that the motion in the invariant plane guides the four dimensional tubes of motion inducing inevitably the crossing of the unstable tube with the stable one in the neighborhood of $`(x=0,p_x=0)`$. Also the eigenvectors of the matrix in (Homoclinic Chaos in the Dynamics of a General Bianchi IX Model) can give us an idea of the behaviour of the flow in the tubes as $`x`$ diminishes, although in the neighborhood of $`x=0`$ the linear expansion (Homoclinic Chaos in the Dynamics of a General Bianchi IX Model)-(Homoclinic Chaos in the Dynamics of a General Bianchi IX Model) is no longer valid. Associated with the eigenvalue $`\lambda =\frac{i\sqrt{2}}{x}`$ we may choose the two independent eigenvectors

$`e_1`$ $`=`$ $`[{\displaystyle \frac{1}{x}},{\displaystyle \frac{2x}{i\sqrt{2}}},0,0]`$
$`e_2`$ $`=`$ $`[0,0,{\displaystyle \frac{1}{x}},{\displaystyle \frac{2x}{3i\sqrt{2}}}],`$
and for the eigenvalue $`\lambda =\frac{i\sqrt{2}}{x}`$,

$`e_3`$ $`=`$ $`[{\displaystyle \frac{i\sqrt{2}}{2x}},x,0,0]`$
$`e_4`$ $`=`$ $`[0,0,{\displaystyle \frac{3i\sqrt{2}}{2x}},x].`$
From the above form of the eigenvectors, we can draw a series of important informations about the behaviour of the four dimensional tubes as they approach the singularity $`x=0`$. As $`x`$ decreases, the sections $`x=const.`$ of the four dimensional tubes about the invariant plane stretches in two directions and contracts in other two, while maintaining their symplectic invariants constant. Also as $`x`$ decreases the oscillations have increasing instantaneous frequency $`\frac{\sqrt{2}}{x}`$. This appears as a characteristic of the dynamical approach to the singularity even for the extension of the tubes to a nonlinear neighborhood of the invariant plane. There is a further fundamental fact that results from the $`x`$-dependence of the above eigenvectors: let us consider for instance the sections $`x=1`$ of the tubes which emanate from the two-tori $`\mathrm{}_{E_0}`$ defined in the linear neighborhood of the critical point $`E`$. Their image projected on the planes $`(y,p_y)`$, $`(z,p_z)`$ can be rescaled to small circles that deforms into ellipses, as $`x`$ is infinitesimally diminished; the ellipse is rotated infinitesimally counterclockwise if the section corresponds to the inferior branch of the curves in the invariant plane ($`p_x<0`$) or clockwise if it corresponds to the superior branch ($`p_x>0`$). This will be fundamental to garantee the transversal crossing of the tubes in the neighborhood of $`x=0`$, as we will discuss. We must however remark that this conclusion is taken in the adiabatic approximation, namely, considering that the variation of $`x`$ with time is adiabatic; this certainly does not occur near the singularity but we assume that no drastic change in the dynamics will alter this behaviour for the actual dynamical tubes. We illustrate this numerically in Fig. 3 by depicting the sections $`x=0.5`$ projected on the planes $`(y,p_y)`$. Figs. 4 show a numerical illustration of the unstable and stable 4-dim tubes emanating from the neighborhood of $`E`$ towards the FRW singularity. We note that, for $`E_{crit}E_0`$ small, the motion proceeds towards the singularity along a tube whose projection on the invariant plane $`(x,p_x)`$ shadows the invariant curve corresponding to $`H=E_0`$, along the upper branch and the lower branch (cf.Fig. 4(b)).

We proceed now to show the homoclinic crossing of the stable and unstable tubes at the neighborhood of the singularity $`(x=0,p_x=0)`$, a phenomenon analogous to Poincaré’s homoclinic tangle, and source of chaos in the modelkoiller . The four following points are essential for this. First, in the five dimensional energy surface $`H=E_0`$, the four dimensional tubes of motion $`R\times S^3`$ actually divide the space, separating the dynamics outside from the dynamics inside the tubes. Second, in the canonical variables $`(x,p_x,y,p_y,z,p_z)`$ introduced above, we have seen that $`(x,p_x)`$ are the variables on the invariant plane, whose phase space picture is depicted in Fig. 1, and the remaining degrees of freedom $`(y,p_y,z,p_z)`$ correspond to elliptic modes of motion (cf. analysis in the neighborhood of the invariant plane). Actually, the tubes have the two dimensional invariant plane as the structure at its center (more properly, one of the curves of the invariant plane corresponding to $`H=E_0`$), about which the flow with the oscillatory degrees of freedom $`(y,p_y,z,p_z)`$ proceeds and in fact the motion in the invariant plane guides the flow, inducing inevitably the crossing of one of the four-dimensional unstable tube with the stable one in a neighborhood of $`(x=0,p_x=0)`$. Third, we note that the invariant plane cannot leave the interior of the four dimensional tubes since they are separation surfaces in the five dimensional manifold $`H=E_0`$, and fourth the property of clockwise rotation or counterclockwise rotation of the sections of the stable and unstable tubes, respectively, in their approach to the FRW singularity. Summing up, in the five dimensional energy surface $`H=E_0`$, the model presents the structure of four dimensional tubes of orbits, emanating from a three sphere along the two saddle directions about the critical point saddle-center-center $`E`$, one stable and one unstable (the stable one corresponding to initially expanding universes). These tubes have the invariant plane curve as the structure at its center and the motion in the invariant plane guides the oscillatory flow towards the singularity $`(x=0,p_x=0)`$, inducing inevitably the crossing of the unstable tube with the stable one in a neighborhood of the singularity; that the crossing is not just a point comes from the conservation of sympletic areas in Hamiltonian dynamics, and the transversality of the crossing comes from the clockwise or counterclockwise rotation of the sections $`(y,p_y)`$ and $`(z,p_z)`$ of the stable and unstable tubes, respectively. If we consider the transversal crossing in a section, say $`x=const.`$, it is not difficult to see that the intersection is a $`S^2`$ manifold. Therefore the intersection manifold will be a 3-dim tube of flow (with topology $`R\times S^2`$) which is contained both in the 4-dim stable tube and in the 4-dim unstable tube, and homoclinic to the $`S^3`$ center-center manifold. It is the equivalent of a homoclinic 1-dim orbit in lower dimensional cases.

Typically, if the 4-dim tubes intersect transversally once they will intersect each other an infinite number of times producing an infinite set of homoclinic orbits which are actually topological 3-dim tubes $`R\times S^2`$. The homoclinic 3-dim tube, which is bi-asymptotic to the $`S^3`$ manifold provides the mechanism for stretching and contraction, giving origin to the homoclinic tangle, which is an invariant signature of chaos in the model wiggins . We just mention that the dynamics near homoclinic orbits is very complex, forming well-known horseshoe structures (cf. wiggins ,koiller ,ozorio  and ozorio1 , and the bibliography therein).

Let us consider the first intersection of the tubes; a part of the orbits inside the unstable tube will enter in the interior of the stable tube and the flow will proceed along the stable tube towards the neighborhood of the critical point $`E`$, from where it will re-enter the unstable tube and proceeds towards the FRW singularity and by a new intersection a part of these orbits will again enter the stable tube and proceed back towards the neighborhood of the critical point $`E`$, and so on, producing an infinite recurrence of the motion for a class of orbits, which will be periodic orbits of very long periods and bounded oscillatory orbits. Another part of the orbits inside the unstable tube will flow along the exterior of the stable tube towards the neighborhood of the critical point $`E`$, and afterwards will escape towards the DeSitter attractor at infinity along the exterior of the unstable tube of the second pair (cf. Fig. 2). This for one single intersection. The same pattern will be reproduced for each of the infinite intersections of the unstable tube with the stable one. In this process, the surface of the tubes, which is in fact a boundary for distinct types of flows, will become more and more stretched and folded, resulting in an intricate structure which will have some important physical consequences for the long time behaviour of the orbits. Indeed, basins of initial conditions for initially expanding universes have a fractal boundary associated with recollapse or escape into the DeSitter configuration, as we show soon.

We must comment that the chaotic aspect of the dynamics described here depends crucially on the presence of the cosmological constant that engenders the structure of a saddle-center-center critical point in phase space together with the homoclinic four dimensional tubes, the transversal crossings of which (in a neighborhood of the FRW singularity) characterize chaos definitely and unambiguously in the model. If we restrict ourselves to the energy surface $`E_0=0`$ we have the Mixmaster universe with a cosmological constant which presents the same structure of chaos. For the limiting case ($`E_0=0,\mathrm{\Lambda }=0`$), the absence of the critical point appears to eliminate the set of homoclinic orbits and periodic orbits of arbitrarily large periods, whose recurrence engenders chaos in the models. In this instance, our approach does not allow to invariantly characterize chaos in the Mixmaster model.

Our system is open (noncompact) with two definite asymptotic exits, namely, collapse and escape to the DeSitter configuration. The code escape/collapse defines basin boundaries in the initial conditions set; these boundaries are associated precisely to the surface of the homoclinic tubes in the set, and correspond to the initial conditions for the homoclinic intersection manifold contained simultaneously in the stable and in the unstable 4-dim tubes. Together with the countable set of periodic orbits of arbitrarily large periods that exist in the neighborhood of each homoclinic orbit and that have the homoclinic orbit as an accumulation set ozorio ozorio1 , they constitute the set of orbits which neither escape or collapse. This set of bounded orbits that are in the boundary of collapse and escape constitutes , in this way, the strange repellor. The fractal dimension of the basin boundary sets are calculated below by a box-counting method.

In other words, our system has the characteristic of a chaotic scattering system, having the saddle-center-center as a chaotic scatterer, with two absolute outcomes consisting of (i) escape into inflation (the DeSitter attractor) or (ii) recollapse to the singularity, for initially expanding universes. The DeSitter attractor defines a way out from the initial singularity to the inflationary phase but this exit to inflation is completely chaotic, namely, small fluctuations in initial conditions may cause the universe to change its asymptotic behaviour from recollapse to escape into inflation and vice-versa. This is due to the presence of chaotic sets in the phase space of the system, as a consequence of the structure of homoclinic tubes and their transversal crossings, which provides an invariant characterization of chaos in the model. One of the questions to be examined here is the characterization of sets of initial conditions for which the DeSitter attractor is attained and establish the character of a chaotic scattering system. To see this chaotic exit to inflation, let us consider initial conditions sets for initially expanding universes in the energy surface $`E_0`$ such that the orbits visit a linear neighborhood of the saddle-center-center critical point $`E`$. In this neighborhood, the Hamiltonian is separable according to (12), namely, $`H=(E_{crit}E_0)+E_{hyp}E_{rot_1}E_{rot_2}=0`$. The partition of the energy $`(E_{crit}E_0)`$ into the energies $`E_{hyp}`$ and $`E_{rot_1}+E_{rot_2}`$ of motion about the critical point will determine the outcome of the oscillatory regime about $`E`$ into collapse or escape to inflation (DeSitter attractor) whether $`E_{hyp}<0`$ or $`E_{hyp}>0`$, respectively. However the nonintegrability of the system, with the consequent homoclinic crossing of the tubes, will cause this partition to be chaotic in general, namely, given an arbitrary initial condition of energy $`E_0`$ we are no longer able to foretell in which of regions $`I`$ or $`II`$ about the saddle-center the orbit will land when it approaches $`E`$. Since $`E_{hyp}=0`$ is a limiting case, for initially expanding universes, the motion inside the tubes corresponds to orbits that will recollapse after reaching the neighborhood of $`E`$, while the motion outside the tubes corresponds to orbits that will escape into the inflationary phase (towards the DeSitter attractor at infinity). In other words, the intrincate crossing and merging of the tubes produces chaotic sets in the phase space of the model, in particular establishing fractal basin boundaries associated with the chaotic exit to inflation. To illustrate this, we finally present a numerical experiment where we obtain a measure of the fractal dimension of basin boundaries in sets of initial conditions by using the box-counting method with the uncertainty code collapse/escape into inflation. In this numerical experiment we evaluate the fractal dimension of portions $`𝒟`$ of phase space, about a point of the separatrix $`S`$ on the invariant plane. The box-counting method, due to Ott and collaborators ott , consists of determining the uncertainty exponent $`\alpha `$ appearing in the scaling law $`fϵ^\alpha `$, where $`ϵ`$ is the uncertainty radius about 5,000 initial conditions taken inside the sets $`𝒟`$, and $`f`$ is the fraction of $`𝒟`$ of uncertain initial conditions with the uncertainty code collapse/escape into inflation. The uncertainty exponent is related to the fractal dimension $`d`$ of the fractal basin boundary by $`d=N\alpha `$, where $`N=5`$ is the phase space energy surface dimensionality. In our numerical experiment we selected a set $`𝒟`$ of homogeneously and randomly distributed initial conditions about the point $`(x=0.4,y=1,z=1,p_x=2.036467455281656,p_y=0,p_z=0)`$ on the separatrix $`S`$, such that the orbits visit a neighborhood of the saddle-center-center $`E`$ before collapse or escape into inflation. In Fig. 5 we depict the scaling law $`fϵ^\alpha `$ in which the small linear portion lies between the saturation region for large $`ϵ`$ and noise for very small $`ϵ`$. The following value for the fractal dimension of the set $`𝒟`$ is found, $`d4.627`$, confirming the chaotic nature of these sets.

###### Acknowledgements.

The authors acknowledge CNPQ/Brazil for partial financial support.


format: text

text:
# A distorted radio shell in the young supernova SN 1986J

## 1 Introduction

Type II supernovae (SNe) are associated with massive stars that have expelled slow, dense winds during their supergiant phase. The stellar explosion drives a shock into this wind, at speeds as high as 20000 $`\mathrm{km}\mathrm{s}^1`$ and temperatures of $`\mathrm{\hspace{0.17em}10}^9`$ K. In addition, a reverse shock propagates back into the stellar envelope at speeds of 500-1000 $`\mathrm{km}\mathrm{s}^1`$ relative to the expanding ejecta. This is the so-called standard interaction model (hereafter SIM; Chevalier 1982, Nadyozhin 1985), and radio, optical, and X-ray emission from Type II supernovae have been usually interpreted within this model. The outgoing shock forms a high-energy-density shell that is responsible for the production of synchrotron radio emission, while the reverse shock accounts for the optical and soft X-ray emission.

SN 1986J in NGC 891 is one of the most radio luminous SNe ever discovered. Indeed, at a distance of $`9.6`$ Mpc (Tully 1998), it had a peak luminosity at 5 GHz about 8 and 13 times greater than SN 1979C and SN 1993J, respectively. The precise date of its explosion is not known, but on the basis of the available radio and optical data SN 1986J was estimated to have exploded around the end of 1982, or the beginning of 1983 (Rupen et al. 1987, Chevalier 1987, Weiler, Panagia & Sramek 1990). Based upon its large radio luminosity, Weiler, Panagia, & Sramek  suggested that the progenitor star was probably a red giant with a main-sequence mass of 20 – 30 $`\mathrm{M}_{}`$ that had lost material rapidly ($`\dot{\mathrm{M}}`$$`2\times \mathrm{\hspace{0.17em}10}^4`$$`\mathrm{M}_{}\mathrm{yr}^1`$) in a dense stellar wind. VLBI observations by Bartel et al.  showed that the radio structure of SN1986J had the form of a shell, in agreement with expectations from the SIM model, with the minimum of emission located approximately at its center.

SN 1986J is peculiar in several respects. (i) Bartel et al.  showed the presence of protrusions in the brightness distribution of the supernova. This finding was interpreted as evidence of deviation from spherical symmetry. Another supernova for which considerable asymmetry in its radio structure has been found is SN1987A (Gaensler et al. 1997, and references therein). Those cases contrast with the remarkable spherically symmetric expansion observed for the radio shell of SN 1993J (Marcaide et al. 1995, 1997; Bartel et al. 2000); (ii) the radio light curves for SN 1986J were not well fitted within the SIM for radio supernovae, which led Weiler et al. (1990) to invoke the existence of a mixed medium (thermal absorbers and non-thermal emitters), or significant filamentation in the circumstellar medium; (iii) optical spectra obtained in 1986 (Rupen et al. 1987) and in 1989 (Leibundgut et al. 1991) showed optical narrow lines with FWHM $``$ 500 - 700 $`\mathrm{km}\mathrm{s}^1`$. In the SIM, the optical and X-ray emission arise in the reverse shock, which moves with a speed close to that of the forward shock, i.e., several thousands $`\mathrm{km}\mathrm{s}^1`$. Since the optical velocities are far lower than those determined from VLBI (radio) observations at similar epochs, these narrow lines question the validity of the SIM to successfully explain the observed optical and X-ray emission from SN 1986J. Chugai  proposed a different model which overcomes the problems of the SIM and seems plausible for SN 1986J. In his model, the supernova envelope is not colliding with a smooth distribution of wind material, but with a clumpy one, and the bulk of the observed X-rays originates in the shocked dense wind clumps. Therefore, the narrow-line ($`\mathrm{\Delta }v\mathrm{\hspace{0.17em}500}\mathrm{km}\mathrm{s}^1`$) material is not directly related to the shock wave (which moves at much higher speeds), but to the speed of the shock-excited dense clouds. More recently, Chugai & Belous  have shown that the evolution of the radio light curves of SN1986J can be well interpreted by free-free absorption in a clumpy wind. The possible presence of clumps of radio absorbing plasma is also consistent with the proposal by Weiler et al. (1990) for a filamentary structure in the circumstellar medium.

Unfortunately, VLBI observations of SN 1986J have been so scarce since its explosion that only one high-resolution radio image for SN 1986J exists (Bartel et al. 1991). In this paper, we use very sensitive 5 GHz global VLBI observations taken on 21 February 1999 to obtain the second high-resolution radio image of SN 1986J, 16 yr after its explosion. We present a brief report of the VLBI observations in section 2, present and discuss our results in section 3, and summarize our conclusions in section 4.

## 2 VLBI Observations and Image Processing

We used archival global VLBI data for SN 1986J. The supernova was observed at a frequency of 5 GHz from 17:00 UT on 21 February 1999 to 04:50 UT on 22 February 1999, using a very sensitive VLBI array that included the following antennas (diameter, location): The whole VLBA (25m, 10 identical antennas across the US), phased-VLA (130m-equivalent, NM, US), Effelsberg (100m, Germany), Medicina and Noto (32m each, Italy), and Onsala (20m, Sweden). The telescopes received both left- and right-hand circular polarizations (LCP and RCP) which, after correlation, were combined to obtain the total intensity mage presented in this paper. Effelsberg had technical problems (damaged gear) and thus could not take part in the observations. Onsala only recorded in LCP mode, and therefore its data were not used to obtain our total intensity VLBI image.

The observations were made with a bandwidth of 64 MHz. The data were correlated at the VLBA Correlator of the National Radio Astronomy Observatory (NRAO) in Socorro (NM, US). The correlator used a pre-averaging time of 3 s. Since SN 1986J was expected to be very faint, at a level of a few mJy, the observations were carried out in phase-reference mode. SN 1986J and the nearby ICRF source, 3C 66A, were alternately observed during the 12-hr experiment. The observations consisted of $`120`$ s scans on SN1986J and of $`70`$ s scans on 3C 66A, plus a few additional seconds of antenna slew time to make a duty cycle of 190 s. 3C 66A was also observed as amplitude calibrator for SN 1986J. The source 3C 84 was observed as a fringe finder.

The correlator data were analyzed using the Astronomical Image Processing System (AIPS). The visibility amplitudes were calibrated using the system temperature and gain information provided for each telescope. The instrumental phase and delay offsets among the 8-MHz baseband converters in each antenna were corrected using a phase calibration determined from observations of 3C 84. The data for the calibrator 3C 66A were then fringe-fitted in a standard manner. We exported the 3C 66A data into the Caltech imaging program DIFMAP (Shepherd et al. 1995) for mapping purposes. 3C 66A showed a flux of $`0.95`$ Jy at 5 GHz, and displayed a one-sided core-jet structure at mas scales, with the jet extending southwards up to $`28`$ mas. The final source model obtained for 3C 66A was then included as an input model in a new fringe-fitting search for 3C 66A. In this way, the solutions obtained were structure-free. The phases, delays, and delay-rates determined for 3C 66A were then interpolated and applied to the source SN 1986J. The SN 1986J data were then transferred into the DIFMAP program. Standard self-calibration techniques and purely uniform weighting were used to achieve maximum resolution in the hybrid maps shown in Fig. 1.

We note two powerful aspects of the use of the phase-reference technique: (i) it effectively increases the integration time on a source from minutes to hours, thus allowing the detection and imaging of very faint objects (e.g. Beasley and Conway 1995), and (ii) it retains the positional information of SN 1986J with respect to the phase-reference source, 3C 66A, almost 1 degree apart. Indeed, a Fourier inversion of the SN 1986J data within DIFMAP, before applying any self-calibration, showed that the a priori position of SN 1986J used in the correlator ($`\alpha =02\stackrel{h}{.}22\stackrel{m}{.}31\stackrel{s}{.}320`$, $`\delta =+42\stackrel{}{.}19\stackrel{}{.}57\stackrel{}{.}282`$; J2000.0) was off from the actual one by $``$17 mas in right ascension (westwards), and by $``$21 mas in declination (southwards). Therefore, phase-referencing proved to be also useful in providing further refinements in the precise coordinates of SN 1986J.

## 3 Results and Discussion

The only other available radio image of SN 1986J is that obtained by Bartel et al. (1991; hereafter B91) of 29 September 1988, at a frequency of 8.4 GHz. Our global VLBI image of 21 February 1999 was obtained at 5 GHz. The most remarkable thing in our image (panel (b) in Fig. 1) is the fact that the brightness distribution of SN 1986J still shows the form of a shell, even though more than 10 yr since the last observations have passed. However, the shell, definitely distorted and with strong asymmetries, shows changes in the brightness distribution between the two epochs.

### 3.1 Asymmetric shell morphology and protrusions in SN 1986J

We distinguish at least four bright features in the rim of the shell labeled A,B, C, and D in Fig. 1. The absolute minimum ($``$0.26 mJy beam-1) of the shell brightness distribution is indicated by the X in panel (b). The flux densities and distances of the bright features A,B,C, and D from X are shown in Table 1. At a flux density of 1.13 mJy beam-1, the corresponding brightness temperature of feature A is $`1.8\times \mathrm{\hspace{0.17em}10}^7`$ K. Features B and C are local maxima with brightness temperatures of $``$1.0 and $``$1.2 $`\times \mathrm{\hspace{0.17em}10}^7`$ K, respectively. Component D is an extended feature in the AC direction. Therefore, its flux density and radial distance from X are just indicative.

The image by B91 showed several ‘knots’ along the rim of the shell of the supernova brightness distribution. The strongest knot was located northeast of the center, with a peak brightness temperature of $`3\times \mathrm{\hspace{0.17em}10}^8\mathrm{K}`$. B91 also showed the presence of protrusions in the brightness distribution of the supernova, at a distance of $``$ 1.9 mas from the centre. We note that the radial distances to features A, B, and C in Fig. 1 are comparable to the radial distances to the protrusions observed by B91, and therefore one could be tempted to identify features A, B, and C with stationary clumps in the CSM of SN1986J. However, caution must be taken, since not only the position angles are different, but also our synthesized beam is significantly larger than that of B91. Our image also shows three spots outside the shell structure at contours above 3 times the rms noise off-source. These spots have distances and position angles with respect to X that are significantly different from those obtained by B91, and therefore do not correpond to the same physical features shown in B91. The peaks of these spots range from 0.23 to 0.25 mJy/beam, less than the minimum of emission within the shell. Therefore, they could be mere artifacts of the image reconstruction procedure, implying that the protrusions seen by B91 have disappeared during the period ranging from September 1988 to February 1999, and pointing to a change in the density properties of the presupernova wind.

Our image also shows that the radio shell morphology of SN 1986J is asymmetric. The four bright features in Fig. 1 seem to delineate a highly distorted shell –whose minimum of emission is not at its centre– and are indicative of a significant deformation of the shock front. If the point X in panel (b) corresponds to the centre of the explosion, we obtain the angular radial distances shown in Table 1 for the knots in the shell of the supernova. It would appear then that the supernova has expanded more along the AC axis. In addition, the angular radial distance XA is more than twice the distance XC, indicative of a strong asymmetric expansion. We note, however, that our identification of X with the centre of the supernova explosion is not unique. In this respect, future phase-reference VLBI observations will be most useful, since this technique will allow a proper registration of the maps, thus permitting the identification of the centre of the explosion.

Nevertheless, most of the results presented and discussed below do not depend on the putative centre of the explosion. In particular, the brightness distribution of SN 1986J is clearly asymmetric, independent of the location of the centre of the explosion. How could such an asymmetric structure have arisen? Various models have been proposed to account for asymmetric morphologies in supernovae. Blondin, Lundqvist & Chevalier  have suggested that an axisymmetric density distribution in the wind from a supernova progenitor leads to protrusions emerging along the symmetry axis. These authors find that for a power-law supernova density profile, the flow approaches a self-similar state in which the protrusion length is 2–4 times the radius of the shell, after $`10`$ yr. They cite the supernova remnant 41.95+575 in M 82 as an example where such axisymmetric protrusions are likely to have emerged, but point out that their model is not compatible with the protrusions seen by B91 for SN 1986J. A case where considerable asymmetry in its overall structure has been observed is SN1987A. Gaensler et al.  showed that the eastern and western regions of the SN1987A radio remnant were brighter than the northern and southern regions, and interpreted this as evidence for the shocked wind being axisymmetric in its shape and/or density distribution. From our image, it follows that an axisymmetric circumstellar interaction does not seem to apply for SN1986J. The protrusions outside the shell, on the other hand, could be formed by an asymmetric distribution in the wind similar to that invoked by Blondin et al. , though not restricted to be axisymmetric. Other models that predict supernova aspherical morphologies include that of Khokhlov et al. , which have modeled jet-induced explosions of core collapse supernovae. The end result is a highly aspherical supernova with two high-velocity jets of material moving in polar directions, and a slower moving, highly distorted ejecta containing most of the supernova material. A recent VLBI image of 41.95+975 by McDonald et al.  reveals structure resembling that of a collimated outflow, as expected from the model by Khokhlov et al. . Although in our case one could identify the brightest knots A and C with a collimated outflow, there is also strong emission in the almost perpendicular direction BD. In addition, all four knots seem well confined within the supernova shell. Therefore, the existence of two high-velocity, well collimated jets in SN 1986J seems unlikely. In the case of SN1987A, a directional anisotropy in the supernova explosion might be indicated by the fact that the eastern limb of SN1987A was brighter, and was also increasing in brightness more rapidly, than the western limb (Gaensler et al. 1997). Although such a directional anisotropy in SN1986J could explain the existence of feature A (which roughly corresponds to the strongest knot of emission seen by B91), it has difficulties in explaining the change in position angles for the other peaks, as well as the appearance of features B and D. We therefore favour an scenario where the strongly asymmetric brightness distribution of the SN 1986J shell structure is due to the collision of the supernova ejecta with a clumpy (Chugai 1993, Chugai & Belous 1999), or filamentary wind (Weiler et al. 1990). Our image gives support to this model, and shows that the clumpy, or filamentary wind must probably be inhomogeneous to produce a highly distorted shell.

### 3.2 Deceleration of the expanding shell

At a distance of 9.6 Mpc (Tully 1998), 1 mas corresponds to a linear size of $`1.4\times \mathrm{\hspace{0.17em}10}^{17}\mathrm{cm}0.05\mathrm{pc}`$. Based on 8.4 GHz VLBI observations, B91 found an angular size of $``$3.7 mas ($`5.3\times \mathrm{\hspace{0.17em}10}^{17}\mathrm{cm}0.17\mathrm{pc}`$) for the shell of SN 1986J, 5.74 yr after its explosion (assuming it took place on 1983.0). The corresponding mean linear velocity would then be $``$14700 km s-1 for the first 5.74 yr. However, this velocity applies only to the protrusions found by B91, not to the shell. Indeed, the velocities reported in B91 were calculated for the protrusions, and assuming that these originated in the centre at the time of the explosion. These authors also pointed out that the protrusions extended from the centre to twice the radius of the shell, i.e., the protrusions were outside the shell. Therefore, a value of $``$1.85 mas ($`2.7\times \mathrm{\hspace{0.17em}10}^{17}\mathrm{cm}0.09\mathrm{pc}`$) for the angular size of the shell of SN 1986J at epoch 1988.74 is indicated, and a mean linear velocity of the radio shell of $``$7400 km s-1 is more appropriate to characterize the first 5.74 yr of the expansion of SN 1986J, as has been previously noticed by Chevalier  and Houck et al. . Chevalier  also pointed out that such a velocity at roughly the time of the peak flux is evidence against the synchrotron self-absorption mechanism acting in SN 1986J.

The angular size of the supernova (at the 5-rms noise level) along the AC direction in Fig.1 is $`\theta 4.7\pm 0.3`$ mas, equivalent to $`6.8\times \mathrm{\hspace{0.17em}10}^{17}`$ cm $`\mathrm{\hspace{0.17em}0.22}`$ pc. Combining this angular size measurement with that obtained by Bartel et al. (1991) for epoch 1988.74 ($`\theta \mathrm{\hspace{0.17em}1.85}`$ mas), we then obtain a mean angular expansion velocity of the shell between 29 September 1988 (1988.74) and 21 February 1999 (1999.14) of $`0.14\mathrm{mas}\mathrm{yr}^1`$, which corresponds to a linear velocity of $`6300\mathrm{km}\mathrm{s}^1`$. If we assume that SN 1986J freely expanded for the first 5.74 yr of its life ($`rt^m,m=1`$) and then started to decelerate, the expansion between the two epochs of VLBI observations is characterized by $`m=0.90\pm 0.06`$. This is a mild deceleration, and contrasts with the case of several other supernovae, for which a strong deceleration has been measured (SN1979C: Marcaide et al. 2002; SN1987A: Staveley-Smith et al. 1993, Gaensler et al. 1997; SN1993J: Marcaide et al. 1997, Bartel et al. 2000). Thus, the angular size that we obtain for the shell seems consistent with the conclusion that the bulk of the shell structure expanded for the first 5.74 yr at speeds smaller than $`15000\text{ km s}\text{-1}`$ (see above and Chevalier 1998).

For a standard presupernova wind velocity, $`v_w`$=10  km s-1, the linear size of SN 1986J at epoch $`t`$=16.14 yr implies that we are sampling the progenitor wind about 11000 yr prior to its explosion. Since the mass loss rate of SN 1986J seems to have been $`2\times \mathrm{\hspace{0.17em}10}^4\text{ }\mathrm{M}_{}\mathrm{yr}^1`$ (Weiler, Panagia & Sramek 1990), the swept-up mass must have been $`\mathrm{M}_{\mathrm{sw}}2.2\mathrm{M}_{}`$, and the thermal electron density $`8\times \mathrm{\hspace{0.17em}10}^3\text{ cm}\text{-3}`$ (for a standard density profile of the progenitor wind, $`\rho _{\mathrm{cs}}r^2`$). Since the expansion of the supernova has not decelerated significantly between the two VLBI observations, it follows that the swept-up mass by the shock front must be much less than the mass of the ejected hydrogen-rich envelope, $`\mathrm{M}_{\mathrm{env}}`$, as otherwise we should have observed a much stronger deceleration. In fact, momentum conservation implies that $`\mathrm{M}_{\mathrm{env}}\mathrm{\hspace{0.17em}12}\mathrm{M}_{}`$, significantly larger than $`\mathrm{M}_{\mathrm{sw}}`$. If the hydrogen-rich mass envelope was as high as 12 $`\mathrm{M}_{}`$, this is a hint that the progenitor of SN 1986J was probably a single, massive Red Super Giant (as previously suggested by Weiler et al. 1990), which lost mass rapidly, but managed to keep intact most of its hydrogen-rich envelope by the time of explosion. This result contrasts with the cases of SN 1993J and SN 1979C, whose hydrogen-rich envelopes had masses of 0.2–0.4 $`\mathrm{M}_{}`$(Woosley et al. 1994, Houck & Fransson 1996) and $`\mathrm{\hspace{0.17em}0.9}\mathrm{M}_{}`$ (Marcaide et al. 2002), and whose progenitor stars were part of binary systems.

### 3.3 Total energy and magnetic field of SN 1986J

Since the radio emission is of synchrotron origin, we can estimate a minimum total energy (in magnetic fields, electrons, and heavy particles) and a minimum magnetic field for SN 1986J. If we assume equipartition (magnetic field energy is approximately equal to the total particle energy), then the minimum total energy is (Pacholczyk 1970) $`E_{\mathrm{min}}^{(\mathrm{tot})}=c_{13}(1+k)^{4/7}\varphi ^{3/7}R^{9/7}L^{4/7},`$ where $`L`$ is the radio luminosity of the source, $`R`$ is a characteristic size, $`c_{13}`$ is a slowly-dependent function of the spectral index, $`\alpha `$, $`\varphi `$ is the fraction of the supernova’s volume occupied by the magnetic field and by the relativistic particles (filling factor), and $`k`$ is the ratio of the (total) heavy particle energy to the electron energy. This ratio depends on the mechanism that generates the relativistic electrons, ranging from $`k1`$ up to $`k=m_p/m_e2000`$, where $`m_p`$ and $`m_e`$ are the proton and electron mass, respectively.

Based on snapshot VLA observations of SN 1986J carried out on 13 June 1999 (Weiler, private communication), we determined an spectral index of $`\alpha =0.69\pm 0.06`$ ($`S_\nu \nu ^\alpha `$) between 1.6 and 8.5 GHz. This value is very close to $`\alpha =0.67`$ obtained by B91, and is a typical value for supernovae that are in the optically thin radio regime. With $`\alpha =0.69`$ and $`S_{4.9\mathrm{GHz}}=7.2`$ mJy from our observations, we obtain a radio luminosity $`L_\nu 1.38\times 10^{37}(D/9.6\mathrm{Mpc})`$ erg s-1 for the frequency range $`10^7`$$`10^{10}`$ Hz. The value of the function $`c_{13}`$ is approximately $`2.39\times 10^4`$ (for $`\alpha =0.69`$, and $`\nu _1`$ and $`\nu _2`$ equal to $`10^7`$ and $`10^{10}`$ Hz, respectively). The filling factor, $`\varphi `$, is by definition $``$ 1. In the case of supernova SN 1993J, Marcaide et al. (1995a, 1995b, 1997) showed that the radio emitting shell of SN 1993J had a width, $`\mathrm{\Delta }R0.3R_{\mathrm{out}}`$, where $`R_{\mathrm{out}}`$ is the outer shell radius, and therefore $`\varphi 0.66`$. As the characteristic size for SN 1986J, we will take half the largest diameter of the shell, $`\theta =2.35\mathrm{mas}`$, which corresponds to a linear size of $`R=D\theta 3.4\times 10^{17}`$cm. With these values, the minimum total energy is then

$$E_{\mathrm{min}}^{(\mathrm{tot})}1.13\times 10^{48}(1+k)^{4/7}\varphi _{0.66}^{3/7}\theta _{2.35}^{9/7}D_{9.6}^{14/7}\mathrm{erg}$$

where $`\varphi _{0.66}=(\varphi /0.66)`$ $`\theta _{2.35}=(\theta /2.35\mathrm{mas})`$ , and $`D_{9.6}=(D/9.6\mathrm{Mpc})`$. The value of the magnetic field that yields $`E_{\mathrm{tot}}^{(\mathrm{min})}`$ is then equal to

$$B_{\mathrm{min}}10.6(1+k)^{2/7}\varphi _{0.66}^{2/7}\theta _{2.35}^{6/7}D_{9.6}^{2/7}\mathrm{mG}$$

Since $`1k2000`$, $`E_{\mathrm{min}}^{(\mathrm{tot})}`$ can have values between $`2\times 10^{48}`$ and $`9\times 10^{49}\mathrm{erg}`$, while the corresponding values of the magnetic field can lie between $`13`$ and $`93\mathrm{mG}`$ (for $`\varphi =0.66`$ and $`D=9.6\mathrm{Mpc}`$). These values of the magnetic field for SN 1986J are in agreement with, e.g., those obtained for SN 1993J at similar radii. For example, Pérez-Torres, Alberdi, and Marcaide  showed that $`B\mathrm{\hspace{0.17em}0.58}(r/3\times \mathrm{\hspace{0.17em}10}^{16}\mathrm{cm})^1`$ G, which results in a value of $`\mathrm{\hspace{0.17em}50}`$ mG for a radius of $`3.4\times \mathrm{\hspace{0.17em}10}^{17}`$ cm. Since it is very unlikely that the magnetic field energy density in the wind is larger than its kinetic energy density, i.e., $`B^2/8\pi \rho v_\mathrm{w}^2/2`$, we can then obtain an upper limit for the magnetic field in the circumstellar wind of SN 1986J:

$$B_{\mathrm{cs}}(\text{ }\dot{\mathrm{M}}v_w)^{1/2}r^1\mathrm{\hspace{0.17em}0.78}(\text{ }\dot{\mathrm{M}}_4v_{10})^{1/2}r_{17}^1\mathrm{mG}$$

where $`\text{ }\dot{\mathrm{M}}_4=\text{ }\dot{\mathrm{M}}/10^4\text{ }\mathrm{M}_{}\mathrm{yr}^1`$, $`v_{10}=v_w/10\text{ km s}\text{-1}`$, $`r_{17}=r/10^{17}\mathrm{cm}`$, and we have assumed a standard wind density profile. For SN 1986J, $`\text{ }\dot{\mathrm{M}}=2\times \mathrm{\hspace{0.17em}10}^4\text{ }\mathrm{M}_{}\mathrm{yr}^1`$ and $`r=\mathrm{\hspace{0.17em}3.4}\times \mathrm{\hspace{0.17em}10}^{17}`$ cm, and we obtain $`B_{\mathrm{cs}}\mathrm{\hspace{0.17em}0.32}`$ mG. These arguments suggest that the magnetic field in the shell of SN 1986J is in the range 13–93 mG, or about 40 to 300 times the magnetic field in the circumstellar wind. Since a strong shock yields a fourfold increase in the particle density, the post-shock magnetic field is also fourfold increased. Hence, compression alone of the circumstellar wind magnetic field cannot account for the large magnetic fields existing in SN 1986J, and other field amplification mechanisms are needed to be invoked, e.g., turbulent amplification (Chevalier 1982; Chevalier & Blondin 1995). The same conclusion was reached for the case of SN 1993J (Fransson & Björnsson 1998, Pérez-Torres, Alberdi & Marcaide 2001), where magnetic field amplification factors $`100`$ were found to be necessary.

## 4 summary

We used 5 GHz global VLBI observations of SN 1986J taken on 21 February 1999 to obtain a high-resolution radio image of the supernova about 16 yr after its explosion (assumed to have occurred on 1983.0). Our main results can be summarized as follows:

(i) The radio structure of SN 1986J is shell-like, even after more than 16 yr since its explosion. However, the shell is highly distorted, indicative of a strong deformation of the shock front.

(ii) Assuming that the shell seen in our 1999.14 image corresponds to the shell in the 1988.74 image, the expansion of SN 1986J has decelerated from $``$7,400  km s-1 down to $``$6,300  km s-1. This mild deceleration can be characterized by a power-law ($`rt^m`$) with a deceleration parameter $`m=0.90\pm \mathrm{\hspace{0.17em}0.06}`$. Since the swept-up mass must be $`2.2\mathrm{M}_{}`$ (for $`\text{ }\dot{\mathrm{M}}=2\times \mathrm{\hspace{0.17em}10}^4\text{ }\mathrm{M}_{}\mathrm{yr}^1`$ and $`v_w=10\text{ km s}\text{-1}`$), this mild deceleration seems to indicate a mass of the hydrogen-rich envelope ejected at explosion of $`\mathrm{\hspace{0.17em}12}`$$`\mathrm{M}_{}`$. In this case, the supernova progenitor kept almost intact its hydrogen-rich envelope, in spite of its strong mass-loss wind rate, and is a strong hint that SN 1986J went off in a single, massive star scenario.

(iii) We show that the brightness distribution in the shell is strongly asymmetric, and suggest that the most likely scenario for its origin is the collision of the supernova envelope with a clumpy, or filamentary, anisotropic wind. However, we cannot exclude that other mechanisms, such as a directional anisotropy in the supernova explosion, may also play a role in shaping the radio emitting structure of SN1986J. We detect three emission spots at the 3-rms noise level outside the shell. Though these spots could be protrusions alike those detected by Bartel et al. (1991), we suggest that they are mere artifacts of the image reconstruction procedure.

(iv) The radio flux of SN 1986J is $`7`$ mJy, which results in a luminosity of about $`1.4\times \mathrm{\hspace{0.17em}10}^{37}\mathrm{erg}\mathrm{s}^1`$. By assuming equipartition between fields and particles, we estimate a minimum energy for the supernova shell in the range $`E_{\mathrm{min}}^{(\mathrm{tot})}`$ 0.2–9 $`\times \mathrm{\hspace{0.17em}10}^{49}`$ erg, which corresponds to a minimum magnetic field of $`B_{\mathrm{min}}`$ 13–93 mG. We find that the average magnetic field in the circumstellar wind of SN 1986J is $`B_{\mathrm{cs}}0.32`$ mG. Therefore, powerful amplification mechanisms (e.g. turbulence) should be acting in SN 1986J, since compression alone of the circumstellar field cannot account for the magnetic fields existing in the supernova shell.

The high-resolution radio image of SN 1986J presented here shows that VLBI is a very powerful technique for investigating the interaction of supernovae, because it allows to directly image the emission from the interaction region (supernova ejecta and shocked wind). Supernova SN 1986J belongs to the rare group of radio supernovae that permit their follow-up with VLBI. Further high-resolution, sensitive VLBI observations of SN 1986J are necessary to provide us with a more detailed view of its evolution and interaction with the circumstellar medium. The use of phase-referencing will be mandatory, since this technique will not only ensure the detection of the supernova, but will also allow a proper registration of the images for all epochs. This registration will likely permit the unambiguous identification of the centre of the explosion. Moreover, such observations would also allow to shed further light on the apparent asymmetric presupernova wind of SN1986J, obtain a better determination of the deceleration parameter, and definitively confirm, or reject, the protrusions outside the supernova shell structure we have marginally detected here.

## Acknowledgments

We thank the National Radio Astronomy Observatory (NRAO) for providing us with the archival data used in this paper. We thank Roger Chevalier, Nikolai Chugai, and Alan Roy for comments on the manuscript. We are grateful to an anonymous referee for many substantive and useful comments, which significantly improved our manuscript. This research has been supported by a Marie Curie Fellowship of the European Community (contract IHP-MCFI-99-1), and by the Spanish DGICYT grants AYA2001-2147-C02-01/02. KWW thanks the Office of Naval Research (ONR) for the 6.1 funding supporting this research. The NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Incorporated. The European VLBI Network is a joint facility of European and Chinese radio astronomy institutes funded by their national research councils.