text:
# Two-dimensional hydrodynamic simulation of an accretion flow with radiative cooling in a close binary system

## 1 Introduction

The angular momentum loss of the gas in the accretion disc in a close binary system is one of the very important phenomena in the astrophysics. The most accepted model that explains the mechanism of angular momentum transport is the $`\alpha `$ disc model (Shakura & Sunyaev, 1973; Pringle, 1981). In this model, viscosity originating from turbulence, magnetism or whatever, is supposed to transport the angular momentum. However, in this model the details of how the viscosity is generated are still unknown.

An alternative model that needs no viscosity for angular momentum transfer is the spiral shock model in which spiral shaped shock waves are formed in the accretion disc due to the tidal force of the mass losing star. This model was proposed by Sawada, Matsuda & Hachisu (1986a). They conducted two-dimensional hydrodynamic calculations of a flow of invisicid adiabatic gas, and found the spiral shocks numerically. Since then a number of two-dimensional simulations have been carried out, and it has been claimed that the spiral shocks really appear in the accretion discs (Sawada, Matsuda & Hachisu, 1986b; Sawada et al., 1987; Spruit et al., 1987; Larson, 1988; Spruit, 1989; Rozyczka & Spruit, 1989; Matsuda et al., 1990; Savonije, Papaloizou & Lin, 1994). Besides, Spruit (1987) obtained self-similar solutions of spiral shocks for idealized accretion discs.

For this spiral shock model, two major points have been questioned as regard to the validity of the model:

The three-dimensional effect might suppress the formation of spiral shocks (Lin, 1989; Lin, Papaloizou & Savonije, 1990a, b; Lubow & Pringle, 1993).

The cooling effect due to radiation is not taken into account explicitly. Though such effect has been considered by assuming a lower specific heat ratio for adiabatic gas, the obtained disc was still too hot.

As regard to the question (i), three-dimensional simulations have been carried out mostly by particle methods. Such simulations tend to give rather diffusive results if a sufficient number of particles are not available. Spiral shocks are likely disappeared particularly for those cases with a lower specific heat ratio for which the pitch angle of the spiral shock decreases. For example, Yukawa, Boffin & Matsuda (1997) carried out three-dimensional simulations using smoothed particle hydrodynamics method. As many as $`(56)\times 10^4`$ particles were employed in the calculations. They showed that an accretion disc was formed and spiral shocks appeared for the case of $`\gamma =1.2`$, but spiral shocks disappeared for the case of $`\gamma =1.1`$ and $`\gamma =1.01`$.

Three-dimensional simulation using either the finite difference or finite volume method was first carried out by Sawada & Matsuda (1992). They calculated the case of $`\gamma =1.2`$ with a mass ratio of unity using the upwind TVD Roe scheme and a generalized curvilinear coordinate system. Although the calculation was continued up to only a half revolution period because of the limitation of CPU time, they successfully showed the existence of spiral shocks on the accretion disc.

Bisikalo et al. (1997a), Bisikalo et al. (1997b), Bisikalo et al. (1998a), Bisikalo et al. (1998b) and Bisikalo et al. (1998c) carried out three-dimensional numerical simulations of accretion discs by means of the finite difference method. They employed a TVD Roe scheme with a monotonic flux limiter of the Osher’s form. A Cartesian coordinate system was used in the calculations. Their results showed that no hot spot was formed. Furthermore, disc formation was inhibited for higher $`\gamma `$ and no spiral shock was seen.

Makita, Miyawaki & Matsuda (1987) carried out two- and three-dimensional numerical simulations of an accretion disc in a close binary system with a higher resolution using Simplified Flux vector Splitting finite volume method (Jyounouchi et al., 1993; Shima & Jyounouchi, 1994). Their computational region only covered the vicinity of the mass accreting star. The gas from the mass losing star was assumed to flow into the computational domain through a rectangular hole located at the L1 Lagrangian point. They obtained quasi steady solutions, and confirmed the existence of both an accretion disc and spiral shocks for all cases with $`\gamma =1.2`$, $`1.1`$, $`1.05`$ and $`1.01`$. It was shown that a smaller $`\gamma `$ value resulted in a more tightly wound spiral shock in two-dimensional calculations, but such tendency was not so obvious in three-dimensional cases. It was also shown that the spiral shock waves disappeared when tidal force due to the mass losing star was artificially cut off.

Fujiwara et al. (2001) carried out three-dimensional simulations of a close binary system containing both the mass accreting star and the mass losing star, and investigated the interaction between the L1 stream and the accretion disc. In the calculations, the ratio of constant specific heat was chosen as $`\gamma =1.2`$ and $`1.01`$. They confirmed the existence of accretion disc as well as spiral shocks. No hot spot was found in the accretion disc. Instead, they found that a bow shock wave was formed due to the collision of the L1 stream and the rotating disc. It was shown that this bow shock wave heated the outer part of the accretion disc, and also enhanced the density perturbation in the disc resulting in a more effective transfer of angular momentum by the tidal torque.

From above numerical results, one can assume the existence of spiral shocks rather firmly even for three-dimensional accretion discs. This view is further supported by a recent observational evidence. Steeghs, Harlaftis & Horne (1997) found the first convincing evidence for spiral structure in the accretion disc of the eclipsing dwarf nova binary IP Pegasi using the technique known as Doppler tomography.

Next, let us examine the question (ii). In the past numerical simulations of accretion discs, it has been customary to assume an adiabatic gas with a lower specific heat ratio than that of a mono-atomic gas ($`=5/3`$) to account for non-adiabatic process of radiative cooling. So far, no simulation that accounts for radiative cooling effect with $`\gamma =5/3`$ has been made to show whether a disc structure is really formed and spiral shocks appear. It is also yet to know how good is the approximation to employ an adiabatic gas with a lower $`\gamma `$ for simulating actual flowfield with radiative cooling effect.

Recently, Stehle (1999) carried out two-dimensional hydrodynamic calculations of accretion discs including both the $`\alpha `$-type viscosity and the effect of energy loss from the surface by radiation. They followed the evolution of the local disc thickness in the one-zone model of Stehle & Spruit (1999). The computational domain, however, only covered the vicinity of the mass accreting star. In their calculation, spiral shocks were clearly observed which dominated the disc evolution in a hydrodynamical time-scale.

The purpose of this work is first to obtain the two-dimensional flowfield in a close binary system with radiative cooling effect, and then to explore whether a disc structure is really formed and spiral shocks appear in the accretion disc. Moreover, we would like to examine how good is the approximation of using a lower $`\gamma `$ value for simulating non-adiabatic processes of radiative cooling in the hydrodynamic calculation of accretion discs. We therefore solve the two-dimensional Euler equations with a constant specific heat ratio of mono-atomic gas ($`\gamma =5/3`$) for a thin disc on the equatorial plane, and include radiative cooling effect by discharging energy in the vertical directions from the top and bottom surfaces of the disc. We employ the flux-limited diffusion (FLD) approximation to obtain radiative heat flux values. In the present study, we consider a realistic case of the observed close binary system.

This paper is organized as follows. In §2 we briefly describe the numerical method, physical assumptions and numerical treatment of radiative transfer. In §3 we show the obtained results for accretion disc with radiative cooling. A detailed comparisons are made with the results of assuming an adiabatic gas with lower $`\gamma `$ values. In §4 concluding remarks will be given.

## 2 ASSUMPTIONS AND METHOD OF CALCULATIONS

We consider a SU UMa type of CV: OY Car as a realistic close binary system (Ritter & Kolb, 1995). The separation between the centers of these two stars is $`4.22\times 10^{10}[\text{cm}]`$. The orbital period is $`1.55[\text{hour}]`$. The mass of these two stars are $`M_1=0.69M_{}`$ for the mass accreting star and $`M_2=0.07M_{}`$ for the mass losing star, respectively. In the present study, we assume the surface density of the mass losing star to be $`\rho _0=0.5\times 10^9[\text{g cm}^3]`$.

### 2.1 BASIC EQUATIONS AND PHYSICAL CONDITIONS

The basic equations describing a two-dimensional flow of perfect gas in a rotating frame of reference can be described as (Sawada, Matsuda & Hachisu, 1984)

$`{\displaystyle \frac{\stackrel{~}{Q}}{t}}+{\displaystyle \frac{\stackrel{~}{F}}{x}}+{\displaystyle \frac{\stackrel{~}{G}}{y}}+\stackrel{~}{H}=0,`$ (1)
$`\stackrel{~}{Q}=\left(\begin{array}{c}\rho \\ \rho u\\ \rho v\\ e\end{array}\right),\stackrel{~}{F}=\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ (e+p)u\end{array}\right),`$ (10)
$`\stackrel{~}{G}=\left(\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ (e+p)v\end{array}\right),`$ (15)
$`\stackrel{~}{H}=\left(\begin{array}{c}0\\ \rho f_x\\ \rho f_y\\ \rho (uf_x+vf_y)+\frac{F_z}{z}\end{array}\right),`$ (20)
where $`\rho `$, $`u`$, $`v`$, $`e`$, $`p`$, $`f_x`$, $`f_y`$ and $`F_z`$ are the density, the $`x`$ and $`y`$ components of the velocity, the total energy per unit volume, the gas pressure, the $`x`$ and $`y`$ components of the net force that includes the gravity force, the Coriolis force and the centrifugal force, and the $`z`$ component of the radiative energy flux, respectively. We ignore the radiation pressure because it is negligibly small if compared with the gas pressure.

The equation of state is given by

$$p=\left(\gamma 1\right)\left\{e\frac{\rho }{2}\left(u^2+v^2\right)\right\},$$
(21)
which is that of an ideal gas characterized by a constant ratio of specific heats $`\gamma `$.

A finite volume approach is used to discretize the governing equations. The AUSM-DV scheme is used to obtain the numerical convective flux (Wada & Liou, 1994). In the time integration, a matrix free LU-SGS method is employed (Nakahashi et al., 1999). Inner iterations are made to assure the second order temporal accuracy in the implicit integration (Yamamoto & Kano, 1996). A CFL number of $`1000`$ is assumed in the calculations.

The parameters that characterize the gas ejected from the mass losing component are the sound velocity, which is chosen to be $`c_0=73[\text{km s}^1]`$, and the velocity of the gas inside of the star, which is $`u_0=0`$. Although we assume a stationary gas, the higher pressure inside of the mass losing star drives the gas flown into the computational domain. Note that if we employ $`A\mathrm{\Omega }`$ as a typical velocity scale as was in the previous works, where $`A`$ is the separation of the two stars and $`\mathrm{\Omega }`$ is the angular velocity of the system, the dimensionless sound velocity of the ejected gas becomes $`0.15`$. As a boundary condition at the surface of the mass accreting star as well as at the outer boundary of the computational domain, we assume a vacuum condition in which any gas reaching there is simply absorbed in. As an initial condition, we assume a gas with very low density, $`0.5\times 10^{14}[\text{g cm}^3]`$, and high temperature, $`3.22\times 10^6/\gamma [\text{K}]`$, that fills up the entire computational domain.

In order to avoid the numerical instability at the initial stage, we first calculate the flow of the gas with $`\gamma =5/3`$ without radiative cooling up to $`10^5`$ time steps using first-order-accurate scheme. After that we switch to second-order-accurate scheme. At the same time, the $`\gamma `$ value is changed for relevant cases. The radiative cooling effect is also included for non-adiabatic cases.

Our computational domain includes both the mass losing star and the mass accreting star on the equatorial plane (Fig. 1). The computational domain that contains the mass losing star has $`20\times 61`$ grid points, while that contains the mass accreting star has $`44\times 61`$ grid points. These two domains are connected by a narrow strip that has $`3\times 19`$ grid points.

### 2.2 NUMERICAL TREATMENT OF RADIATIVE TRANSFER

It is assumed that energy discharge due to radiative cooling occurs in the vertical directions from the top and bottom surfaces of the disc. Because the present flowfield is confined in a two-dimensional equatorial plan, it is necessary to assume the distribution of physical quantities in the vertical direction in the disc for calculation of radiation. We assume a hydrostatic balance in the vertical direction to determine the thickness of the disc, and that the physical quantities in the vertical direction are constant with the local values in the equatorial plan. The thickness of the disc $`H`$ is given by

$$H=\frac{c_s}{\mathrm{\Omega }_\text{K}}=c_s\sqrt{\frac{r_1^3}{GM_1}},$$
(22)
where $`c_s`$, $`\mathrm{\Omega }_\text{K}`$, $`r_1`$ and $`G`$ are the sound speed of the gas on the disc, the Keplerian frequency, the radial distance from the mass accreting star and the gravitational constant, respectively.

The coupling of the flowfield and the radiation field is occurred through the term of $`F_z/z`$ that appears in the energy equation. $`F_z/z`$ is obtained by solving the following equation for radiation energy

$$\frac{F_z}{z}=4\pi \chi Bc\chi E,$$
(23)
where $`B`$, $`E`$, $`c`$ and $`\chi `$ are the Planck function, the radiation energy density per unit volume, the speed of light and opacity, respectively. Because the speed of light is far larger than the sound velocity, and we only seek for a steady state solution, we ignore the time differential term that appears in the radiative transfer equation.

In the present work we employ the Kramer’s opacity which has no dependence on frequency. This Kramer’s opacity is suitable for the calculation of pure hydrogen gas above $`10^4[\text{K}]`$, and is written in the following form (Cannizzo, Shafter & Wheeler, 1988; Stehle, 1999)

$$\chi =2.8\times 10^{24}\rho ^2T^{3.5}[\text{cm}^1],$$
(24)
where $`T`$ is the temperature of the gas.

As mentioned earlier, we consider the radiative transfer only in the vertical direction of the disc, i.e. the $`z`$ direction perpendicular to the $`x\text{-}y`$ plane. The radiative transfer within the disc is ignored. Moreover, the radiative transfer is taken into account in the calculations only in the computational domain that contains the mass accreting star and hence the accretion disc.

The radiative heat flux in the accretion disc is calculated by the FLD approximation (Alme & Wilson, 1974; Levermore & Pomraning, 1981). We assume the gas is in local thermodynamic equilibrium at a temperature that needs not correspond to that of the radiation field. This is because we solve the radiation energy equation (23) to determine the radiation energy density $`E`$. With the FLD approximation, the radiative flux can be written in the form of Fick’s law of diffusion (Levermore & Pomraning, 1981) as

$$F_z=D\frac{E}{z},$$
(25)
with a diffusion coefficient, $`D`$, given by

$$D=\frac{c\lambda }{\chi }.$$
(26)
The dimensionless function $`\lambda =\lambda (E)`$ is called as the flux limiter. We employ Minerbo’s model (Minerbo, 1978) as a choice of $`\lambda `$, i.e. a constraint on the anisotropy of the radiation field. Minerbo assumed a piecewise linear variation of the specific intensity with angle, and found the functional form of the flux limiter as

$$\lambda (R)=\{\begin{array}{cc}2/(3+\sqrt{9+12R^2})\hfill & \text{if}0R3/2,\hfill \\ (1+R+\sqrt{1+2R})^1\hfill & \text{if}3/2<R<\mathrm{},\hfill \end{array}$$
(27)
where $`R`$ is a dimensionless quantity $`R=\left|E/z\right|/(\chi E)`$. In the optically thin limit $`R\mathrm{}`$, the flux limiter gives

$$\underset{R\mathrm{}}{lim}\lambda (R)=\frac{1}{R},$$
(28)
to the first order in $`R^1`$. The magnitude of the flux therefore approaches $`\left|F_z\right|=c\left|E/z\right|/(\chi R)=cE`$, which obeys the causality constraint. In the optically thick or diffusion limit $`R0`$, the flux limiter gives

$$\underset{R0}{lim}\lambda (R)=\frac{1}{3},$$
(29)
to the first order in $`R`$, so that the flux takes the value given by equation (25).

The method of the solving the radiation energy equation (23) is described in Appendix A.

## 3 RESULTS

Figs. 2 (a), (b) and (c) show the density distribution in a (quasi) steady state (at $`10^6`$ step) for the cases of adiabatic gas with $`\gamma =5/3`$, $`1.2`$ and $`1.01`$, respectively. In these calculations radiative cooling effect is not included. It is shown that the disc structure is formed and the spiral shocks are appeared in each case. For the case of $`\gamma =5/3`$, the spiral arms are less clear as shown in Fig. 2 (a). This is because the contour lines are plotted at equal interval while the density in the disc for this case is relatively lower than that in the L1 stream.

It should be noted that the steady state is not steady in a strict sense, but slightly oscillates. Especially, for the case of $`\gamma =5/3`$, the density pattern and the shock position periodically change with time, although the general pattern in the inner and outer region remains unchanged. The temporal fluctuation of the disc structure, however, is not significant and the density distribution shown in Fig. 2 (a) gives a typical one for the case of $`\gamma =5/3`$. This is suggested by the fact that the present density distribution virtually coincides with that of the ensemble average taken from $`2.0\times 10^5`$ to $`1.0\times 10^6`$ time steps.

The pitch angle of spiral arms has a clear correlation with $`\gamma `$. The lower $`\gamma `$ thus leads to a tightly wound spiral shocks. These results are consistent with the previous works (Sawada, Matsuda & Hachisu, 1986a; Makita, Miyawaki & Matsuda, 1987).

The maximum Mach number in the disc is about $`4`$, $`5`$ and $`23`$ for the cases of $`\gamma =5/3,1.2`$ and $`1.01`$, respectively.

Fig. 2 (d) shows the density distribution for the case in which the specific heat ratio is fixed as $`5/3`$ and the effect of radiative cooing is accounted for. Obviously, one can see that the disc structure is formed and spiral shocks are appeared.

Spiral shocks disappear near the mass accreting star and nearly axisymmetric disc structure is formed. Those flow patterns are quite similar to that obtained in the calculation of adiabatic case with $`\gamma =1.01`$, although the absolute value of density in the disc becomes smaller.

The maximum Mach number in the disc is found to be $`23`$ for this case which is the same as that for the case of $`\gamma =1.01`$.

Figs. 3 (a), (b) and (c) show the temperature distribution in a (quasi) steady state with logarithmic scale for the cases of adiabatic gas with $`\gamma =5/3`$, $`1.2`$ and $`1.01`$, respectively. As can be seen in Fig. 3 (a) and (b), the temperature of the gas in the disc is very high and becomes above $`10^7[\text{K}]`$ near the mass accreting star. The smaller $`\gamma `$ value for this case results in nearly isothermal distribution in the disc that gives substantially lower temperature near the mass accreting star. As shown in Fig. 3 (c), compared with the other adiabatic cases, the temperature distribution for the case of $`\gamma =1.01`$ is found to be close to flat.

Fig. 3 (d) shows the temperature distribution with logarithmic scale in a (quasi) steady state for the case of $`\gamma =5/3`$ with radiative cooling. Compared with the other adiabatic cases without radiative cooling, the temperature in the accretion disc is significantly lowered to less than $`10^5[\text{K}]`$, except for the small disc region near the mass accreting star where the temperature exceeds $`10^7[\text{K}]`$. In other words, the accretion disc exhibits a dual structure, i.e. the significantly cooled outer region and the inner core region where temperature is very high. The precise structure in the inner core region, however, is difficult to resolve with the present mesh system, because the pitch angle of the spiral shocks becomes smaller as approaches to the inner region.

The evolution of the total mass in the mass accreting region is shown for each case in Fig. 4. The total mass is defined as a sum of all the masses in the mass accreting region. The horizontal axis denotes the time step. $`1.0\times 10^6`$ time steps correspond to more than $`700`$ orbital periods. In all the cases, the total mass in the mass accreting region becomes almost constant after $`3.0\times 10^5`$ time steps but begin to oscillate slightly around each mean value. Therefore, we can say that each flowfield in the accretion disc reaches a (quasi) steady state. The total mass for the case accounting for radiative cooling is found to be larger than that for the case of $`\gamma =1.2`$, and smaller than that for the case of $`\gamma =1.01`$. From the results, we can say that the overall behavior of the flowfield as well as the evolution of the total mass of the adiabatic case with a smaller $`\gamma `$ value is quite similar to the one with radiative cooling shown in the present study. Although the precise $`\gamma `$ value can change for case-to-case, one can conclude that the use of a lower $`\gamma `$ value for simulating radiative cooling effect is almost a valid assumption.

Finally, let us examine the influence of gas properties from the mass losing star. Because the process of discharge and absorption of radiation depends on density and temperature of the gas in the accretion disc that is influenced primarily by the gas from the mass losing star, we need to know what are the consequences of changing properties of the gas. Fig. 5 (a) shows the density contours in the accretion disc for the lower temperature case in which $`\rho _0=0.5\times 10^9[\text{g cm}^3]`$ and $`c_0=48[\text{km s}^1]`$. In the calculation, the radiative cooling effect is considered. In this case, the accretion disc slightly shrinks, but the overall feature is unchanged if compared with the standard case shown in Fig. 5 (d). On the other hand, if we increase the surface density of the mass losing star, some changes are seen in the results. Fig. 5 (b) shows the density contours for the case with $`\rho _0=0.5\times 10^8[\text{g cm}^3]`$ and $`c_0=73[\text{km s}^1]`$. Note that the density scale is changed to show the overall distribution. In this case, the inner disc region again appears where density becomes very low, but its radius becomes smaller. Surrounding this inner region, there appear several spiral arms in the intermediate region. This intermediate region has a ring like structure, which is connected to an outer spiral arm and also to the incoming flow from the L1point. In Fig. 5 (c), the density contours for the limiting case of $`\rho _0=1.0\times 10^7[\text{g cm}^3]`$ and $`c_0=24[\text{km s}^1]`$ are shown. One can see that the overall feature is quite similar to that shown in Fig. 5 (b), though the inner disc region further shrinks.

## 4 CONCLUDING REMARKS

Two-dimensional hydrodynamic calculations of an invisicid flowfield in a close binary system are carried out by solving the Euler equations with radiative transfer. In the present study, the specific heat ratio is assumed to be constant while radiative cooling effect is included as a non-adiabatic process. The cooling effect of the disc is considered by discharging energy in the vertical directions from the top and bottom surfaces of the disc. We use the flux-limited diffusion approximation to calculate the radiative heat flux values. The obtained results are summarized as follows:

Around the mass accreting star, there appears an accretion disc and spiral shocks are formed in the disc, even for the case with radiative cooling.

The total mass of the disc reaches a (quasi) steady value. It becomes larger when radiative cooling effect is accounted for than that assuming $`\gamma =5/3`$, and is close to the value for the case of lower $`\gamma `$.

The lower $`\gamma `$ value for an adiabatic gas which has been employed in the simulation of accretion discs in order to account for radiative cooling effect is found to be almost a valid assumption in terms of the overall flow features and also the amount of the total mass.

What are shown in this work are the results of two-dimensional hydrodynamic calculation of inviscid flowfield accounting for radiative cooling effect. In order to obtain more rigorous results for the accreting flowfield, we need to explore three-dimensional flowfield that is coupled with radiative transfer. In such case, one needs to account for radiative transfer that is parallel to equatorial plane. This is because the accretion disc is not necessarily a flat thin disc. It is also needed to employ a more detailed radiation model that can consider spectral dependence of opacity. For such calculation, a flux-limited diffusion model is no more applicable and we need further to employ a detailed simulation method such as to use ray-tracing technique to solve radiative transfer equations. Our future studies will focus on these aspects.

## Appendix A Discretizing of Radiation Energy Equation

Let us consider the method of solving the radiation energy equation (23). Using the flux-limited diffusion approximation (25), equation (23) can be rewritten in the following form for radiation energy density $`E`$ as

$$\frac{}{z}\left(D\frac{E}{z}\right)=4\pi \chi Bc\chi E.$$
(30)
Because the distribution of physical quantities in the vertical direction is assumed to be constant, the diffusion coefficient $`D`$ is independent of the height $`z`$. Therefore, we replace the partial differential in $`z`$ by multiplying $`1/H`$. As a result, equation (30) is discretized as

$$2D^n\frac{E^{n+1}}{(H^n)^2}=4\pi \chi ^nB^nc\chi ^nE^{n+1},$$
(31)
where the superscript $`n`$ denotes the time step. The factor $`2`$ in the left hand side of the equation (31) appears because discharging energy by radiation occurs from the top and bottom surfaces of the disc. Except for the radiation energy density $`E`$, the dependent variables are evaluated at time step $`n`$. Equation (31) is solved for $`E^{n+1}`$. This $`E^{n+1}`$ is substituted into the discrete form of $`F_z/z`$, i.e. $`2D^nE^{n+1}/(H^n)^2`$ which is further substituted into the energy equation of the Euler equations. This completes the coupling procedure of the flowfield with radiation.


format: text

text:
# Discrete Symmetries Underlying Some Continuous Ones: Two Examples From Gravity And Particle Physics

## 1 Introduction

In the standard model of particle physics , we see a very complecated primary postulate about the kind of continuous internal symmetry group. In addition, multiplets in different generations should be transformed by the interwoven representations of the group. It is obvious that such the complecated postulates are not good for a fundamental theory.

In this paper, the author describes two own examples from particle physics and gravity in which only some kind of discrete symmetry takes place for the initial model equations but one has a possibility to introduce a continuous symmetry on this base. There is a common feature for the both examples: one must have some additional coordinates or parameters to assume the new simple demands on solutions with a discrete symmetry obeying transition to a continuous one.

As the first example, we consider here a model of composite fermions . The second example is an embedding the general relativity $`4`$space into a flat $`12`$space that results in an alternative model of gravity .

## 2 Symmetries of the composite fermions <br>model

By linearization of an equation for an energy $`E`$ of a two-component system

$$E=(m_1^2+\overline{p}_1^2)^{1/2}+(m_2^2+\overline{p}_2^2)^{1/2}$$

(and assuming after it that constituents’ masses $`m_1=m_2=0`$; $`\overline{p}_1`$ and $`\overline{p}_2`$ are their momentums), one can get the system of linear quantum equations (13)-(16) from  for a wave function $`\psi (x_1,x_2),`$ where $`x_1`$ and $`x_2`$ are the coordinates of both constituents. We deal with two $`8`$ spaces in the model: $`(x_1,x_2)`$ and $`(x,y),`$ where $`x`$ belonging to the Minkowski space $`R_{1,3}^4`$ are coordinates of a centre of inertia, and $`y`$ are internal coordinates. The author has assumed in  that $`y`$ are transformed independently $`x.`$

The equations of motion of such the composite system have the following algebraic property: they permit eight different solutions of the kind $`G_A\psi (x,y)`$ if $`\psi (x,y)`$ is some solution. Matrices $`G_A,A=1,\mathrm{},8,`$ set up a representation of the discrete group $`D_4.`$ If one assumes that transitions between these solutions are induced by transformations of a space $`(y),`$ it leads that an algebraic structure of field equations puts hard restrictions on this space: the space $`(y)`$ should be discrete.

There are two possibilities: to have $`y0`$ for all $`A`$ or to have two isolated sets of solutions when $`y0.`$ To get the global symmetry group $`SU(3)_c\times SU(2)_l`$ for the model in the latter case, one must introduce an additional postulate about conservation of the norm of a set of solutions (see , section 7). The internal coordinates $`y`$ are not observable; namely it leads to the apparent continuous system’s symmetry for any observer in the $`(x)`$ space. Some other features of the standard model on global level are reflected by this model automatically. A minimal set of generations contains $`4`$ generations. A multiplet of any generation contains two $`SU(2)`$singlets that gives a possibility to introduce non-zero neutrino masses that is important after their observation by the Super-Kamiokande collaboration .

## 3 Symmetries of the model of gravity in a flat $`12`$space

To embed the general relativity curved $`4`$space into a flat $`12`$space $`(x,A,B)`$ with flat $`4`$sectors $`(x),(A),(B),`$ one may map $`(x)(A)(B)`$ and get the connection $`\mathrm{\Gamma }=\stackrel{~}{\mathrm{\Gamma }}+\overline{\mathrm{\Gamma }}`$ with $`\stackrel{~}{\mathrm{\Gamma }}`$ being its tensor part. This connection has a trivial curvature. After it, we may linear deform the connection: $`\mathrm{\Gamma }\gamma =f\stackrel{~}{\mathrm{\Gamma }}+\overline{\mathrm{\Gamma }}`$ where $`f`$ is a scalar parameter. A non-trivial curvature of the connection $`\gamma `$ is proportional to $`(f^2f).`$ The last step is to oblige the according Ricci tensor to Einstein’s equations to get an alternative model of gravitation in the flat $`12`$space .

In this case, $`f`$ is $`U(1)`$symmetry’s parameter; the parameter $`Ff^2f`$ gives a curvature scale. Two-valueness of the mapping $`Ff`$ gives $`D_1`$symmetry which may be transformed into the $`SU(2)`$one as in the previous case. But the parameters of the $`U(1)`$ and $`SU(2)`$transformations will depend from each other in a general case. We have here a very interesting analogy with the standard model: to provide its independence, one should do a rotation on some angle $`\theta `$ in the parameter plane $`(f,F)`$ . If we take in the mind a necessity to unify the two described models then the minimum value $`\theta _{min}`$ is determined by the demand that one component of the $`SU(2)`$doublet should stay massless or near to it. We will have in such the case: $`\mathrm{sin}^2\theta _{min}=0,20`$ from the pure geometrical reasons. $`\theta _{min}`$ is an analog of the Weinberg angle $`\theta _w`$ for which $`\mathrm{sin}^2\theta _w=0,215`$ from an experiment. It is exiting that these values of angles $`\theta _{min}`$ and $`\theta _w`$ are approximate enough. After the rotation, the $`SU(2)`$symmetry will be broken if we take gravitation into account. I would like to note here that as it was shown in my paper , quantum gravity would be super-strong on small distances of the order $`10^{11}m.`$ It would mean that one cannot consider it alone on this scale.

## 4 Conclusion

In the two cases, which are described here, one can use some constructive (algebraic in both situations) features of the model yielding discrete symmetries to get some continuous ones. Today, these two models are not connected between themselves. I think that to do the first step to unify them, one should linearize Eqs. (9) or (11) from  which are uniform relative to the connection $`\stackrel{~}{\mathrm{\Gamma }}.`$ As the second step, I consider now a possibility to assume that super-strong interacting gravitons  are constituents of the composite fermions.


text:
# Reply to ‘Comment on ‘Gravitating Magnetic Monopole in the Global Monopole Spacetime’ ’

In a recent paper A. Achúcarro and J. Urrestilla  have pointed the stability problem for composite, global and magnetic monopoles defect, presented in a previous publication . In the former it is claimed that this topological defect is not classically stable against axially symmetric angular deformation on the boson field associated with the global sector ***The stability problem of the global monopole against angular perturbation has been observed by A. S. Goldhaber , and also by the authors in .. By this deformation an extra tension is created at the north-pole which drags the core of the global monopole upwards with no cost of energy. Once the cores of both monopoles are separated, they are repelled by an induced self-interaction consequence of the distortion on the magnetic fields due to the solid angle deficit $`\delta \mathrm{\Omega }=4\pi \mathrm{\Delta }`$, being $`\mathrm{\Delta }`$ the parameter associated with the energy scale where the gauge symmetry is spontaneously broken. However, in the composite defect system this analysis must take into account the induced self-energy more carefully. In order to infer the behavior of this self-energy I briefly comment previous analysis for the ideal case. In  it was calculated the formal expression to the electrostatic self-energy associated with a test charged particle in the pointlike global monopole spacetime. It was observed that it is positive, and by numerical evaluation can be seen that the induced electrostatic strength increases for larger values of the parameter $`\mathrm{\Delta }`$. Although in  the authors said that the effect of the core has been considered, this is not true. In  the metric tensor considered was associated with the global monopole with no internal structure. Considering the effect of the core of the global monopole, the calculation of the self-energy becomes much more complicated. Although the exact expression for this self-energy becomes impossible to be obtained, we can infer its behavior as following: at the global monopole center, $`r=0`$, the metric tensor is Minkowski. Its component $`g_{11}`$ changes from $`1`$ to $`1/(1\mathrm{\Delta })`$ for distance $`rr_{gm}`$, increasing the value of an effective solid angle deficit. On the other hand, the distortion on the electric field decreases for larger separation. So for the case where increasing the distance the augment of the solid angle deficit becomes dominant, the distortion increases with the separation up to $`r_{gm}`$; consequently we can infer that this self-energy goes to zero for very small separation between the charged particle and the core of the monopole, peaks around the distance of the order of the monopole’s size, $`r_{gm}`$, and decreases approximately with $`1/r`$ for larger value of separation. Assuming that this effect is presented to magnetic charged particle as the authors also admit, to separate the cores of the global and magnetic monopoles needs some amount of energy to overcome this barrier An estimation for the peak of the magnetostatic self-energy can be provided considering the expression for self-energy given in  and that the monopole’s size is of order $`\eta ^1\lambda ^{1/2}`$ .. So the question about the classical stability condition of the composite defect against axially symmetric perturbation must take into account this fact.

Another point raised by Achúcarro and Urrestilla is about the profile of the Higgs field. In fact for the model analysed in  there are three mass parameters to scale the radial distance. In our numerical analysis we adopted the mass of the boson vector: $`x=e\eta r`$. This choice of mass scale seemed convenient for us to analyse the behavior of the fields varying the self-coupling $`\lambda `$, keeping the electric coupling $`e`$ fixed.

Acknowledgment
This work was partially supported by CNPq.


text:
# THE DYNAMIC BEHAVIOR OF SOFT GAMMA REPEATERS

## 1  INTRODUCTION

Soft Gamma Repeaters (SGRs) are a small class (4 confirmed and one candidate) of high-energy transient discovered through their emission of bright X-ray/$`\gamma `$-ray bursts which repeat on timescales of seconds to years. Originally confused with Gamma-Ray Bursts (GRBs see e.g. Fishman and Meegan 1995), SGRs get their name from their burst properties that distinguish them from classical GRBs; namely their softer spectral energy distribution and repetition of bursts.

In addition to being prolific burst emitters, SGRs are also persistent X-ray sources. For three objects, the X-ray radiation is modulated by the spin of the underlying neutron star. In each case, the star is spinning down rapidly, indicative of a strong magnetic field (Kouveliotou et al. 1998a). In fact, it was this observation that clinched the now widely accepted magnetar model for SGRs (Thompson and Duncan 1995, 1996). A magnetar is a strongly magnetized neutron star ($`B_{\mathrm{dip}}`$ $`10^{14}10^{15}`$ G) where the magnetic energy exceeds all other sources of free energy in the system. Unlike ordinary radio pulsars, the rotational energy loss in SGRs is insignificant compared to their overall energy output. It is the decay of the strong field of the SGR that powers both their burst and persistent emissions.

In addition to the four confirmed SGRs, there is another small class of neutron stars dubbed the Anomalous X-ray Pulsars (AXPs) that were proposed by Thompson and Duncan (1996) as magnetar candidates. The AXPs share a number of characteristics with the SGRs in terms of their persistent X-ray emission properties, and two members were recently discovered to emit SGR-like bursts (Gavriil et al. 2002; Kaspi et al. 2002). The recent burst detection results solidify the common nature of AXPs and SGRs. A review of AXPs will be presented within this volume (Gavriil et al. 2003).

Here, I will review the salient properties of SGRs including burst characteristics, persistent and pulsed X-ray emission, and the effects of burst activity on SGR persistent emission properties. Associations with supernova remnants and massive star clusters will be discussed elsewhere in this volume (Gaensler 2003).

## 2  BURST PROPERTIES

The defining behavior of SGRs is their repetitive emission of luminous soft $`\gamma `$-ray bursts. Typically, the bursts last $``$0.1 sec and have energy spectra (E $`>`$25 keV) that can be modeled with an Optically Thin Thermal Bremsstrahlung (OTTB) having temperatures 20$``$40 keV (e.g. Göğüş et al. 2001, Aptekar et al. 2001). At lower energies, however, this model fails to fit the spectrum (e.g. Olive et al. 2003). The spectral temperatures show little or no variation with intensity, between bursts, at different epochs, or between sources. The burst energies follow a power-law number distribution (dN/dE $``$ E-5/3 \[Cheng et al. 1996, Göğüş et al. 1999) from $``$1035 ergs up to at least $``$1042 ergs, consistent with a so-called self-organized critical system (e.g. earthquakes, Solar flares, etc.) where the burst energy reservoir greatly exceeds the energy output within any given burst.

Burst active states of SGRs, often referred to as outbursts, range in duration from a few days to several months. During these episodes, they emit anywhere from a handful to several hundred bursts. The burst rate histories of the four SGRs spanning three decades are shown in Figure 1. These lightcurves are comprised of burst detections from several different detectors having different sensitivity limits. The non-uniformity of the burst sampling becomes significant when comparing separate active episodes that were viewed with different detectors. For example, given a stationary burst energy (and peak flux) distribution with a power-law slope of $`5/3`$, an increase in the detector sensitivity by a factor 10 yields a factor $``$5 increase in the overall rate.

On two occasions, more energetic bursts or giant flares were recorded one each from SGR 0526$``$66 on 1979 March 5 (Mazets et al. 1979) and SGR 1900$`+`$14 on 1998 August 27 (Hurley et al. 1999a). The lightcurve of the August 27 flare is shown in Figure 2. Each of these extraordinary events had a bright ($``$1044 ergs s-1), spectrally hard initial spike followed by a softer, several minute long tail showing coherent pulsations at 8 and 5 s, respectively. More recently, an intermediate flare ($``$1043 ergs) lasting 40 s was recorded from SGR 1900$`+`$14 on 2001 April 18 (Guidorzi et al. 2001), suggesting a continuum of bursts energies rather than a dichotomy of bursts and flares (Kouveliotou et al. 2001).

Recently, discrete features in the spectra of bursts from two SGRs, 1900$`+`$14 and 1806$``$20 have been reported. In the case of SGR 1900$`+`$14, an emission line at $``$6.5 keV was discovered (Strohmayer and Ibrahim 2000) during the onset of a single, bright burst. No candidate line features in other bursts have been reported from this source since. For SGR 1806$``$20, an absorption line feature was found near $``$5 keV in a small subset of spectra taken from selected bursts (Ibrahim et al. 2002). The presence of less significant features at higher energies consistent with a harmonic relationship is intriguing, but their significance (for individual bursts) is not convincing.

## 3  PERSISTENT EMISSION

All SGRs are associated with persistent X-ray counterparts; three of them have quiescent luminosities $``$1035 ergs s-1 (Murakami et al. 1994; Rothschild et al. 1993; Hurley et al. 1999b), while the quiescent flux level of SGR 1627$``$41 has not yet been reached (C. Kouveliotou, private communication). The spectra of SGR 1627$``$41 and SGR 1806$``$20 can be modeled with a power-law (Woods et al. 1999a; Mereghetti et al. 2000) while SGR 1900$`+`$14 requires an additional blackbody component ($`kT`$ $``$0.5 keV \[Woods et al. 1999b; Kouveliotou et al. 2001\]). The presence of a blackbody component in SGR 0526$``$66 is marginal (Kulkarni et al. 2003). Currently, we are investigating the stability of the power-law photon index in SGR 1900$`+`$14. Our preliminary finding is that the photon index evolves from $``$1.0 to $``$2.2 within a few years.

Three SGRs show low-amplitude pulsations in their persistent emission within a narrow period range (Kouveliotou et al. 1998a; Hurley et al. 1999b; Kulkarni et al. 2003). The frequency of these pulsations is increasing rapidly (Kouveliotou et al. 1998a, 1999; Kulkarni et al. 2003), consistent with the interpretation of an underlying strongly magnetized neutron star. Due to the strong timing noise present in both SGR 1806$``$20 and SGR 1900$`+`$14 (Woods et al. 2002), it is best to give ranges for their measured period derivatives (Table 1).

## 4  BURST-INDUCED EFFECTS ON THE PERSISTENT EMISSION

During the last few years, we have noted changes in the X-ray emission properties of SGRs during episodes of burst activity (Woods et al. 2001). Through studying the transient effects imparted upon SGRs (or the lack thereof) during times of burst activity, we have gained deeper insight into the nature of the burst mechanism and the SGR systems in general. The changes in SGR pulsed properties and persistent X-ray emission and their relationship to burst activity are presented below.

### 4.1 Torque Changes

Coherent pulsations from the persistent emission of SGR 1806$``$20 were discovered within an RXTE PCA observation from 1996 November (Kouveliotou et al. 1998a). Including both archival and subsequent monitoring observations, the spin frequency history of SGR 1806$``$20 has now been extended from 1993 through 2001 (Figure 3, from Woods et al. 2002). We have found that at all times, the SGR has been spinning down, but the rate of spindown shows substantial variability. In fact, the measured spin-down torque on this SGR has been found to vary by up to a factor $``$6. The torque variations seen in SGR 1806$``$20 do not correlate directly with the burst activity (Woods et al. 2002).

Pulsations from the X-ray counterpart of SGR 1900$`+`$14 were discovered during an ASCA observation in 1998 April (Hurley et al. 1999b), shortly before the SGR entered an intense, sustained burst active interval (Figure 1). Similar to SGR 1806$``$20, a compilation of data preceeding and following the discovery observation showed that this SGR was spinning down rapidly and irregularly (Kouveliotou et al. 1999; Woods et al. 2002) . As with SGR 1806$``$20, the variations in torque do not directly correlate with the burst activity from this SGR with one notable exception, the giant flare of August 27 (Woods et al. 1999c; Palmer 2002).

We note that although a rapid spin-down event most likely occurred within the hours following the August 27 flare, its impact on the overall spin history of the SGR was very small relative to the much larger variations observed during $``$5 years of monitoring. So, in general, the direct effects of burst activity are insignificant to the overall torque noise in each of these SGRs. For a more complete discussion of the torque variability in these two SGRs, as well as a quantitative analysis of the torque noise, see Woods et al. (2002).

Recently, the 8 sec pulsations observed in the tail of the March 5th flare from SGR 0526$``$66 (Mazets et al. 1979) were confirmed in the persistent X-ray emission of the source (Kulkarni et al. 2003). As of yet, there are only two period measurements which provides a measure of the spindown. Further data are required to determine if the torque varies in this system as well.

### 4.2 Pulse Profiles

Currently, the pulse profiles of both SGR 1900$`+`$14 and SGR 1806$``$20 are very nearly sinusoidal (i.e. they show very little power at the higher harmonics). This has not always been the case, however, as both SGRs have shown significant changes in their pulse profiles during the last several years. The most notable of these changes occurred in the pulse profile of SGR 1900$`+`$14 during the tail of the giant flare of August 27 (e.g., Feroci et al. 2001).

Forty seconds after the onset of the August 27 flare, 5.16 s coherent gamma-ray pulsations at high amplitude emerged. Initially, the pulse profile was complex, having four distinct maxima per rotation cycle. A few minutes later toward the end of the flare, the pulse profile was significantly more sinusoidal (Figure 4 – middle row). The same qualitative behavior was observed in the persistent X-ray pulse profile of SGR 1900$`+`$14. In all observations prior to the August 27th flare, the pulse profile was complex having significant power at higher harmonics (Figure 4 – top row). For all observations after 1998 August 27 through the present, the pulse profile remained relatively simple (Figure 4 – bottom row). Hence, the pulse profile change observed at gamma-ray energies during the tail of the August 27th flare translated to the persistent emission pulse profile of this SGR in a sustained manner (i.e. for years after the August 27 X-ray tail had disappeared \[Woods et al. 2001\]).

A systematic study of the temporal and spectral evolution of the pulse profiles of SGR 1900$`+`$14 and SGR 1806$``$20 is presented in Göğüş et al. (2002). This study confirms the temporal evolution of the pulse profile of SGR 1900$`+`$14 mentioned above and shows that there is a similar time-dependence in the pulse profile of SGR 1806$``$20. From 1996 November to 1999 January, the pulse profile of SGR 1806$``$20 becomes more sinusoidal. Due to the sparseness of the observations, however, we cannot determine the exact time of this change, nor the timescale over which it progressed to better than 2.3 years.

### 4.3 Flux Changes

Changes in the flux of SGRs were first noted in SGR 1900$`+`$14 following the giant flare of August 27 (Remillard et al. 1998; Kouveliotou et al. 1999; Murakami et al. 1999; Woods et al. 1999b). Following this discovery, a compilation of persistent and pulsed flux measurements over several years (Figure 5 \[Woods et al. 2001\]) revealed that, in general, there is an excellent correlation between burst activity (top) and enhancements in the persistent/pulsed flux from this SGR (middle). We have found that the pulse fraction is consistent with remaining constant at most epochs despite changes in the persistent flux. It is by assuming that this fraction remains constant at all times that we can plot both the pulsed flux (RXTE PCA) and the persistent flux (BeppoSAX and ASCA) on the same scale. We note, however, that there are exceptions to this rule when the pulse fraction has increased for short periods of time (see below).

#### 4.3.1 Four Bursts from SGR 1900$`+`$14

The brightest pulsed/persistent flux excess seen in Figure 5 is directly linked with the August 27 flare. The flux decays approximately as a power-law in time ($`Ft^{0.7}`$) following the giant flare (Figure 6 \[Woods et al. 2001\]). The X-ray spectrum as measured with the PCA one day after the flare can be fit with the sum of a blackbody and a power-law. The blackbody temperature is hotter ($`kT=0.94\pm 0.03`$) than the quiescent level and the power-law photon index is steeper ($`\mathrm{\Gamma }=2.76\pm 0.07`$). The spectrum (0.1$``$10 keV) of the X-ray tail at $``$19 days after the flare was found to be exclusively non-thermal (Woods et al. 1999b). These observations are consistent with the thermal excess fading more rapidly to its quiescent value than the power-law component.

For SGR 1900$`+`$14, there are now four X-ray tails that can be linked with specific bursts or flares. The second of these events was recorded on 1998 August 29 (Ibrahim et al. 2001; Lenters et al. 2003), just two days after the flare. This burst had a high gamma-ray fluence and an X-ray tail whose flux decayed more rapidly than a power-law in time (Lenters et al. 2003). We define the onset of the tail emission when the intensity quickly drops by $``$4 orders of magnitude within a few tens of milliseconds. Simultaneously, the X-ray spectrum suddenly hardens. Over the next several minutes as the flux decays through the tail, the energy spectrum softens. Formally, the energy spectrum of the tail is equally well fit by a power-law plus a blackbody or an optically thin thermal bremsstrahlung, each with interstellar attenuation (Ibrahim et al. 2001). However, the bremsstralung model yields a column density $``$5 times larger than the measured column from the persistent emission whereas the two component model fit yields a column consistent with the persistent emission value. The pulsed fraction increases above the quiescent level (11% RMS) up to $``$20% during this tail (Lenters et al. 2003), and the phase of the pulsations does not shift during the tail relative to the pre-burst pulse phase.

On 2001 April 18, a burst was detected from SGR 1900$`+`$14 with a high fluence ($``$1043 ergs) and long duration $``$40 s (Guidorzi et al. 2001). The burst energy of April 18 was intermediate between the giant flare of 1998 August 27 and typical SGR events. Like the August 27 flare, this event showed an extended X-ray tail following the burst which lasted for several days (Kouveliotou et al. 2001; Feroci et al. 2003). The properties of this burst tail (Table 2) are similar in many respects to the two previously mentioned.

A second burst from SGR 1900$`+`$14 also recorded during the 2001 April activation also possessed an X-ray tail. This event was detected on April 28 and like previous tails, showed a transient increase in pulse fraction and a cooling blackbody spectral component (Lenters et al. 2003). In fact, this is the only tail thus far to show an exclusively thermal X-ray tail.

Within the group of four X-ray tails following bright bursts detected from SGR 1900$`+`$14 (Table 2), we find a number of similarities and some slight differences indicating possible trends (within the restrictions of small number statistics). In all cases, the thermal component becomes hotter while the power-law component brightens following three of the four bursts. The ratio of tail energy to burst energy is consistent with being constant among the group (Figure 7 \[Lenters et al. 2003\]). Timing anomalies and pulse profile changes are seen following the two most energetic bursts and the pulse fraction increases for the three weakest. An interesting trend which arises from this small set of X-ray tails is that the pulse fraction enhancement in the separate X-ray tails appears to correlate with the magnitude of the thermal contribution to the X-ray flux. That is, tails with the highest relative blackbody flux show the largest increase in pulse fraction. Given the small numbers, however, more examples are required to establish this trend.

#### 4.3.2 SGR 1627$``$41

SGR 1627$``$41 was discovered in 1998 June when it was observed to burst more than 100 times in a short span of about a month (Kouveliotou et al. 1998b; Woods et al. 1999a). No bursts from this SGR have been seen since 1998 August. Shortly after the onset of burst activity in this SGR, an X-ray counterpart was identified and observed to fade in brightness over the next two years. Dissimilar to the observed flux variations in SGR 1900$`+`$14, the decay in the lightcurve of SGR 1627$``$41 proceeds on a much longer timescale of order years and has yet to “bottom out” at a baseline flux level. This object shows that the response of SGRs to burst activity varies between sources.

## 5  DISCUSSION

We have summarized the recent observations of dynamic behavior in the persistent and pulsed emission from SGR 1900$`+`$14, SGR 1806$``$20 and SGR 1627$``$41. Now, we will discuss what constraints these observations place on the models for the SGRs, in particular the magnetar model.

The torque enhancements in SGR 1806$``$20 and SGR 1900$`+`$14 do not directly correlate with the burst activity. In the context of the magnetar model, the abscence of a direct correlation between these two parameters has strong implications for the underlying physics behind each phenomenon. The magnetar model postulates that the bursting activity in SGRs is a result of fracturing of the outer crust of a highly magnetized neutron star. Furthermore, the majority of models proposed to explain the torque variability in magnetars invoke crustal motion and/or low-level seismic activity (Thompson and Blaes 1998; Harding et al. 1999; Thompson et al. 2000). Since there is no direct correlation between the burst activity and torque variability, then either ($`i`$) the seismic activities leading to each observable are decoupled from one another, or ($`ii`$) at least one of these phenomena is not related to seismic activity (Woods et al. 2002).

The dramatic change in the pulse profile of SGR 1900$`+`$14 in conjunction with the giant flare requires a substantial change in the magnetic field of the neutron star (Woods et al. 2001; Thompson et al. 2002). In the magnetar model, there are at least two possible ways that this can happen. One possibility proposed by Thompson et al. (2002) is that a twist in the magnetosphere is generated following the flare, driving a persistent current which produces an optically thick scattering screen at some substantial distance ($``$10 $`R_{}`$) from the stellar surface. In this model, the surface field geometry remains complex at all times. The pulse profile, however, simplifies when the scattering screen is present (i.e. after the flare). The scattering screen must have the properties of redistributing the radiation in phase, but not in energy in order to account for the reemergence of the blackbody component after the August 27 tail fades away (Thompson et al. 2002). The decay of this magnetospheric twist is believed to be several years. An alternative scenario involves restructuring of the surface magnetic field geometry. In this picture, the field geometry is complex prior to the flare and relaxes to a more dipolar structure following the event giving rise to the observed change in pulse profile.

Thompson et al. (2002) recently investigated each of these scenarios in detail, noting advantages and disadvantages for each model. In this work, they have identified further observational tests involving the energy spectrum of the emission before and after the flare. Simulations of the expected behavior and an analysis of the spectral evolution of SGR 1900$`+`$14 are currently underway.

Currently, only four clear X-ray tails have been detected from SGR 1900$`+`$14. As mentioned earlier, there is a potential correlation between the relative abundance of thermal emission in the tail and the enhancement of the pulse fraction. This correlation, if proven correct with the detection and analysis of several more SGR tails, would provide a strong argument for heating of a localized region on the neutron star during bursts (Thompson et al. 2000; Lyubarsky et al. 2002).

Since the pulse fraction increases during three of the four X-ray tails detected from SGR 1900$`+`$14, the flux enhancement must be anisotropic about the star in these cases. For the August 29 and April 28 burst tails, the location of the heating is also constrained (Lenters et al. 2003). In each of these events, we have precise pulse phase information prior to, during, and after the tail. For both bursts, the phase of the pulsations during the tail does not shift relative to the pulse phase prior to the burst. This requires that the localized region on the neutron star with the largest relative flux enhancement is within or nearby the region giving rise to the persistent X-ray pulse peak (e.g. the polar cap). In the magnetar model, the burst emission is due to the build up of stress in the stellar crust from the evolving magnetic field and the eventual release of this stress when the crust fractures (Thompson and Duncan 1995). Localized heating of the polar cap requires that the fracture site of the burst be at the same location. Since these two events do not occur at all near the same pulse phase (Lenters et al. 2003), the burst emission arising from near the surface must be significantly scattered. An alternative model (Lyutikov 2002) invokes a magnetospheric rather than crustal origin for the bursts. This model is attractive in that it negates the need for large angle scattering of the burst emission and provides a natural explanation for the preferential brightening of the polar cap regions.

Contrary to the short-lived X-ray tails seen following SGR 1900$`+`$14 bursts, the much longer timescale flux decay in SGR 1627$``$41 following its 1998 outburst shows that the effects of burst activity on the persistent emission of SGRs is not uniform among the class. As another example, the recent outburst of the AXP 1E 2259$`+`$586 (Kaspi et al. 2002) shows several source parameters that changed as a result. Some of the observed changes (e.g. persistent and pulsed X-ray flux increases) are similar to those seen in SGR 1900$`+`$14 whereas others are not. As discussed above, episodes where the properties of these objects are rapidly evolving are paramount in determining the nature of these sources and the physical mechanism driving the dynamics of the source spectral and temporal properties. The diversity in SGR/AXP outbursts observed thus far necessitates detailed study of each outburst in order to obtain a comprehensive picture of the class as a whole. This in turn requires diligent monitoring and prompt follow-up of each source when they become active.

## 6 Acknowledgements

I thank the many collaborators who have contributed to the results discussed here. I acknowledge financial support from the Long Term Space Astrophysics program (NAG 5-9350). I would like to thank Chryssa Kouveliotou for useful discussions and a careful reading of the manuscript.

E-mail address of P.M. Woods peter.woods@nsstc.nasa.gov


text:
# Solitons, solitonic vortices, and vortex rings in a cylindrical Bose-Einstein condensate

## I Introduction

The possible occurrence of solitons in a Bose-Einstein condensate (BEC) was theoretically predicted sometime ago, within a one-dimensional (1D) Gross-Pitaevskii model tsuzuki ; zakharov , but experimental observation was delayed because of the absence of a physical realization of a strictly 1D Bose gas. Nonetheless, similar coherent structures were recently observed in BECs confined in traps of varying geometry burger ; denschlag  by a method (phase imprinting) that was directly inspired by the analytical structure of 1D solitons. Yet it has become clear that quasi-1D solitons are susceptible to various instabilities within the three-dimensional (3D) environment of realistic traps.

Indeed, a stability analysis based on the linear Bogoliubov-de Gennes (BdG) equations muryshev1  revealed that an axisymmetric soliton is stable only at sufficiently low particle number or high aspect ratio in an elongated trap. The primary mode of instability was shown to result from non-axisymmetric perturbations, in the sense that a purely imaginary eigenvalue first appears in the spectrum of the BdG equations with azimuthal angular momentum $`m=1`$. This ”snake instability” was further analyzed in Ref. feder  and argued to be responsible for a possible decay of a soliton into vortex rings and/or vortices. In fact, based on the above theoretical work, non-stationary vortex rings were experimentally detected in a spherical trap anderson . On the other hand, Brand and Reinhardt brand1 ; brand2  suggested that the primitive mode associated with the snake instability is a stable structure that may be called solitonic vortex.

The emerging picture is sufficiently perplexing to deserve closer attention. In some recent work komineas1 ; komineas2  we have carried out a detailed calculation of axisymmetric solitons in a cylindrical trap, which were shown to become unstable at high particle density even if we restrict ourselves to radial ($`m=0`$) perturbations. The soliton was then shown to bifurcate to an axisymmetric vortex ring with lower energy. Here, we complete this work by allowing azimuthal deformations of the soliton, in order to eventually account for the primary ($`m=1`$) instability discussed in Ref. muryshev1 ; feder . We thus obtain a reasonably complete description of the two basic modes of instability of the soliton and the corresponding bifurcations to axisymmetric vortex rings and non-axisymmetric solitonic vortices.

The simplest theoretical picture is obtained in the ideal limit of an infinitely elongated cylindrical trap, because static solitary waves may then be treated on equal footing with those propagating at nonzero velocity. Hence, all calculations presented in the main text pertain to an infinite cylindrical trap. In Sec. II we formulate the problem and briefly describe the numerical methods employed. The three basic types of static solitary waves (solitons, solitonic vortices, and vortex rings) and their relative stability are discussed in Sec. III, whereas solitary waves moving along the cylindrical trap at constant velocity are calculated in Sec. IV. The main conclusions are summarized in Sec. V, and some important technical details are relegated to two appendices. In particular, the case of a finite axisymmetric trap is discussed in Appendix B.

## II formulation

An axisymmetric harmonic trap is characterized by a transverse confinement frequency $`\omega _{}`$ and a longitudinal frequency $`\omega _{}`$. We first consider the special limit of a cylindrical trap, with $`\omega _{}=0`$, a restriction that is relaxed only in Appendix B. In this limit, physical units are rationalized as in Ref. komineas2 . Thus time is measured in units of $`1/\omega _{}`$ and distance in units of the transverse oscillator length $`a_{}=(\mathrm{}/M\omega _{})^{1/2}`$ where $`M`$ is the mass of each atom. Complete specification of the system also requires as input the linear density $`\nu `$ which is equal to the average number of particles per unit length of the cylindrical trap. An important dimensionless parameter is then given by $`\gamma =\nu a`$ where $`a`$ is the scattering length. A field rescaled according to $`\mathrm{\Psi }\sqrt{\nu }\mathrm{\Psi }/a_{}`$ satisfies the rationalized Gross-Pitaevskii equation

$$i\frac{\mathrm{\Psi }}{t}=\frac{1}{2}\mathrm{\Delta }\mathrm{\Psi }+V_{\mathrm{tr}}\mathrm{\Psi }+g|\mathrm{\Psi }|^2\mathrm{\Psi },$$
(1)
where $`V_{\mathrm{tr}}=(x^2+y^2)/2=\rho ^2/2`$ is the trapping potential and $`g=4\pi \gamma `$ is a dimensionless coupling constant that is the only free parameter. The conserved energy functional associated with Eq. (1) is given by

$$W=\frac{1}{2}\left(|\mathbf{}\mathrm{\Psi }|^2+2V_{\mathrm{tr}}|\mathrm{\Psi }|^2+g|\mathrm{\Psi }|^4\right)𝑑V$$
(2)
and yields energy in units of $`\gamma _{}(\mathrm{}\omega _{})`$ where $`\gamma _{}=\nu a_{}=\gamma (a_{}/a)`$.

We have also considered various forms of dissipative dynamics, mainly as relaxation algorithms for the calculation of stationary states. The simplest possibility is to make the replacement $`i\mathrm{\Psi }/t(i\zeta )\mathrm{\Psi }/t`$ in Eq. (1), where $`\zeta `$ is some dissipation constant. An extreme special case widely employed in numerical simulations feder  is the fully-dissipative dynamics obtained by the replacement $`i\mathrm{\Psi }/t\mathrm{\Psi }/t`$. However, in either case, particle number conservation is violated and special precautions are necessary to enforce a definite number of particles.

An interesting alternative that preserves particle number was originally introduced by Carlson carlson  for the study of the homogeneous BEC and is simply adapted to the case of a confined BEC. Actually, we have employed a special case of Carlson’s dissipative dynamics governed by

$$i\frac{\mathrm{\Psi }}{t}\zeta \mathrm{\Psi }\frac{}{t}|\mathrm{\Psi }|^2=\frac{1}{2}\mathrm{\Delta }\mathrm{\Psi }+V_{\mathrm{tr}}\mathrm{\Psi }+g|\mathrm{\Psi }|^2\mathrm{\Psi }.$$
(3)
An immediate consequence of Eq. (3) is the continuity equation

$$\frac{n}{t}+\mathbf{}𝑱=0,$$
(4)
irrespective of the specific value of the dissipation constant $`\zeta `$. Here, $`n=|\mathrm{\Psi }|^2`$ is the local particle density and

$$𝑱=\frac{1}{2i}(\mathrm{\Psi }^{}\mathbf{}\mathrm{\Psi }\mathbf{}\mathrm{\Psi }^{}\mathrm{\Psi })=n\mathbf{}\varphi $$
(5)
is the usual current density, where $`\varphi `$ is the phase of the wave function $`\mathrm{\Psi }=\sqrt{n}e^{i\varphi }`$. Therefore, the total number of particles is indeed conserved by the evolution equation (3). However, the total energy defined from Eq. (2) dissipates according to

$$\frac{dW}{dt}=\zeta \left(\frac{n}{t}\right)^2𝑑V,$$
(6)
as long as the particle density $`n`$ remains time dependent.

Although both types of dissipation discussed above are employed in our numerical calculations, Carlson’s proposal has some advantages for our current purposes. We have thus developed a discrete version of Eq. (3) on a cubic lattice to test the stability of solitons and eventually obtain rough estimates of other possible stationary states. These estimates are then fed into a Newton-Raphson algorithm of the type considered in Ref. komineas2  to obtain our final results more accurately as well as efficiently. The latter algorithm is currently upgraded to a fully-3D iterative scheme that does not necessarily enforce axial symmetry. Explicit results are discussed in the following sections.

## III static solitary waves

We first consider a special class of stationary states obtained by making the substitution $`\mathrm{\Psi }(𝒓,t)e^{i\mu t}\mathrm{\Psi }(𝒓)`$ in Eq. (1) to derive the static equation

$$\mu \mathrm{\Psi }=\frac{1}{2}\mathrm{\Delta }\mathrm{\Psi }+\frac{1}{2}\rho ^2\mathrm{\Psi }+g|\mathrm{\Psi }|^2\mathrm{\Psi },$$
(7)
where $`\mu `$ is the chemical potential (in units of $`\mathrm{}\omega _{}`$) required to enforce a definite number of atoms. The simplest and most important solution of Eq. (7) is the ground state wave function $`\mathrm{\Psi }=\mathrm{\Psi }_0(\rho )`$ which depends only on the radial distance $`\rho =(x^2+y^2)^{1/2}`$ and may be chosen to be real and positive. It also satisfies the constraint

$$_0^{\mathrm{}}2\pi \rho 𝑑\rho \mathrm{\Psi }_0^2(\rho )=1,$$
(8)
by a specific choice of the chemical potential $`\mu =\mu (\gamma )`$ for each value of the dimensionless coupling constant $`\gamma `$. Explicit illustrations of the ground state may be found in Ref. komineas2 .

A more general class of excited self-consistent states is obtained by solving Eq. (7) with the same chemical potential, as is appropriate on an infinitely elongated trap, and boundary conditions

$$\underset{z\pm \mathrm{}}{lim}|\mathrm{\Psi }(x,y,z)|=|\mathrm{\Psi }_0(\rho )|,\underset{z\pm \mathrm{}}{lim}\frac{\mathrm{\Psi }}{z}=0.$$
(9)
The class of solutions actually discussed in the present paper is further restricted by the parity relation

$$\mathrm{\Psi }(x,y,z)=\mathrm{\Psi }^{}(x,y,z),$$
(10)
modulo an overall constant phase that is set equal to zero. Relation (10) is obviously consistent with Eq. (7).

To conclude the general description, we also consider the excitation energy defined as

$$E=WW_0\mu (|\mathrm{\Psi }|^2|\mathrm{\Psi }_0|^2)𝑑V,$$
(11)
where both $`W`$ and $`W_0`$ are calculated from the energy functional (2) applied for a general self-consistent state $`\mathrm{\Psi }`$ and the ground state $`\mathrm{\Psi }_0`$, respectively. The presence in Eq. (11) of a term that is proportional to the chemical potential $`\mu `$ ensures that we are calculating the energy difference between two states with the same number of particles.

In recent work komineas1 ; komineas2  we restricted ourselves to axisymmetric states described by a wave function of the form $`\mathrm{\Psi }=\mathrm{\Psi }(\rho ,z)`$. The simplest possibility is a static (black) soliton whose wave function is odd under the parity reflection $`zz`$, and purely imaginary in view of the phase convention adopted in Eq. (10). The excitation energy of a static soliton calculated from Eq. (11) is illustrated in Fig. 1 as a function of the dimensionless coupling constant $`\gamma `$.

When $`\gamma `$ exceeds the critical value $`\gamma _1=4`$, a new class of axisymmetric solitary waves emerges with the standard characteristics of vortex rings komineas1 ; komineas2 . Unlike the case of a homogeneous BEC, a vortex ring may now be static and is described by a complex wave function that satisfies the parity relation $`\mathrm{\Psi }(\rho ,z)=\mathrm{\Psi }^{}(\rho ,z)`$. The calculated excitation energy of a static vortex ring is also shown in Fig. 1 and is seen to smoothly bifurcate from the soliton branch at the critical coupling $`\gamma _1`$, while it remains consistently lower than the soliton energy for $`\gamma >\gamma _1`$. Therefore, one may conclude that a radial ($`m=0`$) perturbation of a soliton will render it unstable when $`\gamma >\gamma _1`$ and eventually transform it into an axisymmetric vortex ring.

One is tempted to presume that the soliton is stable for $`\gamma <\gamma _1`$. However, the linear stability analysis of Ref. muryshev1 ; feder  indicates that the primary soliton instability sets in at a lower critical coupling $`\gamma _0<\gamma _1`$ through azimuthal ($`m=1`$) perturbations. To probe for such a possibility, we solve the dissipative Eq. (3) with initial condition supplied by the calculated static soliton which is found to remain stable at all times when $`\gamma <\gamma _0=1.5`$. However, for $`\gamma >\gamma _0`$, the inherent numerical noise introduces an instability that transforms the soliton into a non-axisymmetric structure which is subsequently rectified by feeding it into a 3D Newton-Raphson iterative algorithm that is designed to directly solve the static equation (7).

We thus obtain a static solitary wave that is partially illustrated in Fig. 2 for $`\gamma =3`$ through density plots over two planes that are perpendicular to each other and both contain the $`z`$ axis. We have furthermore examined the phase to establish that the calculated wave function describes a deformed vortex that pierces through the trap in a direction that is perpendicular to its symmetry ($`z`$) axis. In Fig. 2, the vortex axis is identified with the $`x`$ axis as a matter of convention, whereas a degenerate class of equivalent configurations may be obtained by simple azimuthal rotations around the $`z`$ axis. The resulting solitary wave displays the characteristics of a solitonic vortex calculated earlier by Brand and Reinhardt brand1 ; brand2  in traps of varying geometry.

The energy of the solitonic vortex is depicted in Fig. 1 as a function of $`\gamma `$ and is seen to emerge from the soliton branch at the critical coupling $`\gamma _0=1.5`$, while it remains smaller than the energy of both the soliton and the vortex ring for all $`\gamma >\gamma _0`$. The physical picture can be explained by returning to the initial-value problem associated with the dissipative Eq. (3) using the static soliton as initial condition. As mentioned already, the soliton remains stable for $`\gamma <\gamma _0`$ but decays into a solitonic vortex when $`\gamma _0<\gamma <\gamma _1`$. The pattern is more complex for $`\gamma >\gamma _1`$ where the soliton initially deforms radially to become an axisymmetric vortex ring which eventually decays into a stable solitonic vortex.

Therefore, a vortex ring is strictly speaking unstable. But our real-time simulations indicate that vortex rings are sufficiently long lived to be observed, as was actually the case in the experiment of Ref. anderson  – albeit in the special limit of a spherical trap. On the other hand, it is curious that the predicted solitonic vortex has not yet been seen in experiment.

The picture is completed with a few remarks. One should expect that there exists a series of additional critical couplings (densities) at which the soliton bifurcates into multiple vortex rings and/or solitonic vortices. For example, a double vortex ring was calculated in Ref. komineas2  for $`\gamma >\gamma _2=12`$. In other words, we have only described the two basic patterns of instability that correspond to the appearance of unstable modes in the $`m=0`$ and $`m=1`$ sectors of the linear BdG equations. These two modes may be dominant at intermediate densities, but more complicated patterns should be expected to occur at higher densities, as is evident in the experiment of Ref. anderson  which is further discussed in Appendix B.

## IV moving solitary waves

All three types of static solitary waves discussed in the preceding section are special members of a more general class of stationary states that propagate along the cylindrical trap with constant velocity $`v`$. These are obtained by making the substitution $`\mathrm{\Psi }(𝒓,t)e^{i\mu t}\mathrm{\Psi }(x,y,\xi )`$ in Eq. (1), with $`\xi =zvt`$, to derive the stationary differential equation

$`\mu \mathrm{\Psi }iv{\displaystyle \frac{\mathrm{\Psi }}{\xi }}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Delta }\mathrm{\Psi }+{\displaystyle \frac{1}{2}}\rho ^2\mathrm{\Psi }+g|\mathrm{\Psi }|^2\mathrm{\Psi },`$ (12)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{^2}{x^2}}+{\displaystyle \frac{^2}{y^2}}+{\displaystyle \frac{^2}{\xi ^2}},`$
which reduces to the static equation (7) at $`v=0`$. The boundary conditions (9) and the parity relation (10) are also generalized by the simple replacement $`z\xi `$.

Equation (12) is solved by feeding a static ($`v=0`$) solitary wave into the Newton-Raphson algorithm and by upgrading the velocity typically in steps of $`\delta v=\pm 0.1`$. The results for axisymmetric solitons and vortex rings already presented in Ref. komineas1 ; komineas2  have been confirmed using the currently employed fully-3D algorithm. Therefore, in the following, we shall mostly discuss the case of a moving solitonic vortex for which no calculation has been given in the literature. The static solitonic vortex of Fig. 2 is significantly reshaped when it moves, as is illustrated in Fig. 3 for two values of the velocity ($`v=0.4`$ and 0.8) that are both smaller than the speed of sound $`c=1.3`$ calculated for $`\gamma =3`$ komineas2 . The $`v=0.4`$ entry of Fig. 3 makes it clear that a moving solitonic vortex is shifted off center. While the vortex axis remains within the basal $`xy`$ plane, it is clearly displaced away from the center of the trap. The displacement increases with increasing speed until we reach a critical speed $`v_0=0.8`$ where the vortex disappears from the trap, in the sense that the solitary wave is converted into a gray axisymmetric soliton. This implies, in particular, that the density plots over the $`xz`$ and $`yz`$ planes become identical, as is apparent in the $`v=0.8`$ entry of Fig. 3.

For $`v_0<|v|<c`$ the calculated solitary wave is indeed identical to a moving axisymmetric soliton and reduces to a weakly nonlinear soundlike pulse in the limit $`|v|c`$. A corollary of this discussion is that a high-speed soliton can be stable even for couplings $`\gamma `$ that are greater than the critical coupling $`\gamma _0=1.5`$, in apparent agreement with a similar conclusion reached in Ref. muryshev2  that is now made more precise. Specifically, for $`\gamma <\gamma _0`$ a solitonic vortex does not exist and the soliton is stable over the entire range of velocities $`0<|v|<c`$, as expected. For $`\gamma >\gamma _0`$, the soliton is stable over the limited range $`v_0<|v|<c`$, where $`v_0=v_0(\gamma )`$ is the critical speed at which the solitonic vortex is converted into a soliton and $`c=c(\gamma )`$ is the calculated speed of sound komineas2 . Clearly, $`v_0(\gamma =\gamma _0)=0`$, while a numerical investigation of the ratio $`v_0(\gamma )/c(\gamma )`$ shows that it approaches unity with increasing $`\gamma `$. However, we are not in a position to ascertain whether or not there exists a finite (critical) $`\gamma `$ above which $`v_0=c`$ and, hence, the soliton is not stable for any velocity. In any case, the solitonic vortex is the fundamental mode as long as $`\gamma >\gamma _0`$.

Further insight into the nature of a moving solitonic vortex is obtained by considering the behavior of its phase $`\varphi =\varphi (x,y,\xi )`$ in the asymptotic limits $`\xi \pm \mathrm{}`$. Let $`\varphi _\pm =\varphi (x,y,\xi =\pm \mathrm{})`$ be the two asymptotic values of the phase, which are, in principle, some functions of $`x`$ and $`y`$. However, detailed numerical evidence suggests that $`\varphi _+`$ and $`\varphi _{}`$ are, in fact, constants that depend only on the velocity $`v`$. This conclusion may be substantiated in part by an analytical argument briefly described in Appendix A. Therefore, the asymptotic phase difference

$$\delta \varphi =\varphi _+\varphi _{}$$
(13)
depends only on the velocity $`v`$ and characterizes a solitary wave. Incidentally, the above conclusion applies to all three types of solitary waves discussed in the present paper and is pertinent to their experimental production through phase imprinting burger ; denschlag .

As an illustration, we return to the solitonic vortex calculated for $`\gamma =3`$ and plot $`\delta \varphi `$ as a function of velocity by a solid line in Fig. 4, whereas the corresponding result for an axisymmetric soliton at $`\gamma =3`$ is shown by a dashed line. In both cases, $`\delta \varphi =\pi `$ for $`v=0`$, and $`\delta \varphi =0`$ (mod $`2\pi `$) in the extreme limits $`v\pm c`$ which are not actually reached in Fig. 4 due to numerical inaccuracies. The same figure provides clear evidence for the conversion of a solitonic vortex into an axisymmetric soliton at the critical speed $`v_0=0.8`$ or $`v_0/c=0.6`$.

The preceding discussion of the asymptotic phase is also pertinent to the definition of the impulse of a solitary wave komineas2  which is now simplified to

$$Q=𝒫\delta \varphi ,𝒫=n\frac{\varphi }{z}𝑑V,$$
(14)
where $`𝒫`$ is the usual linear momentum and $`\delta \varphi `$ is given by Eq. (13). The impulse $`Q`$ is measured in units of $`\mathrm{}\nu =\gamma _{}(\mathrm{}/a_{})`$.

A solitary wave may then be characterized by the impulse $`Q`$ and the excitation energy $`E`$ defined in Eq. (11). Both $`E`$ and $`Q`$ are definite functions of the velocity $`v`$ that can be eliminated to yield the dispersion $`E=E(Q)`$. A nontrivial check of consistency is provided by the group velocity relation $`v=dE/dQ`$ which was accurately confirmed by our numerical calculation. For comparison purposes, we reproduce in Fig. 5 the dispersion of a stable axisymmetric soliton calculated for $`\gamma =1<\gamma _0`$ in Ref. komineas2 . This dispersion displays the characteristics of a mode originally derived by Lieb lieb  within a full quantum treatment of a homogeneous 1D Bose gas based on the Bethe Ansatz, and was later rederived by a semiclassical calculation based on the solitary wave of the 1D Gross-Pitaevskii model kulish ; ishikawa . A similar mode was also derived within an effective 1D model that is designed to approximately describe a cylindrical trap jackson .

We are now in a position to calculate the dispersion of a moving solitonic vortex. The dispersion is shown by a solid line in Fig. 6 and is seen to retain the basic characteristics of a Lieb-type mode. We further plot by a dashed line the dispersion of an axisymmetric soliton which is expected to be unstable at $`\gamma =3`$. However, the two dispersions shown in Fig. 6 merge at two critical points that correspond to velocities $`v=\pm v_0`$ where $`v_0=0.6c`$ is the critical speed at which the solitonic vortex is converted into an axisymmetric soliton.

A subtle detail not seen in our limited illustrations is that the calculated soliton for $`\gamma =3`$ exhibits some ringlike features, when $`v0`$, which are not present for $`v=0`$ and also disappear before reaching the critical speed $`v_0`$. These features could be viewed as preliminary signals for the formation of fully-fledged vortex rings that occur for stronger couplings in the region $`\gamma >\gamma _1=4`$.

Since a detailed calculation of both static and moving vortex rings was already described in Ref. komineas1 ; komineas2 , we merely complete the picture by illustrating in Fig. 7 the dispersions of a soliton, a vortex ring, as well as a solitonic vortex calculated for $`\gamma =7`$ where all three possibilities occur according to the bifurcation pattern of Fig. 1. The solitonic vortex is still characterized by a Lieb-type dispersion that merges with the dispersion of the vortex ring at the critical speed $`v_0=0.8c`$, where $`c=1.6`$ is the speed of sound calculated at $`\gamma =7`$. For $`v_0<|v|<c`$ both the solitonic vortex and the vortex ring become identical to a gray axisymmetric soliton. In fact, the vortex ring transforms into a soliton at a speed somewhat smaller than $`v_0`$.

Therefore, it is fair to say that the primary mode associated with the snake instability of a soliton is indeed the solitonic vortex brand1  whose calculation has been significantly extended in the present paper to account for steady motion along a cylindrical trap. Actually, one may envisage a further generalization to a spinning solitonic vortex that undergoes rigid-body precession around the $`z`$ axis at constant frequency $`\omega `$, while it moves along the same axis at constant velocity $`v`$. One could then calculate the excitation energy $`E=E(\omega ,v)`$, the impulse $`Q=Q(\omega ,v)`$, and the azimuthal angular momentum $`L_z=L_z(\omega ,v)`$ to eventually derive a dispersion of the form $`E=E(Q,L_z)`$. The Lieb-type dispersion of a solitonic vortex shown by a solid line in Figs. 6 and 7 would then be the lower edge of a continuum of states parametrized by the angular momentum $`L_z`$. This interesting possibility will not be discussed further in the present paper.

## V conclusion

A phase imprint $`\delta \varphi `$ applied between the two ends of an elongated trap leads to the formation of pure quasi-1D solitons only for weak couplings (densities) in the region $`\gamma <\gamma _0=1.5`$. For intermediate couplings in the region $`\gamma _0<\gamma <\gamma _1=4`$, solitons may be produced temporarily but they eventually decay into stable solitonic vortices. A new pattern emerges for $`\gamma >\gamma _1`$ where a phase imprint should initially produce unstable axisymmetric solitons that subsequently deform radially to become axisymmetric vortex rings and finally decay into solitonic vortices. This pattern is further complicated by the fact that high-speed vortex rings or solitonic vortices are indistinguishable from axisymmetric solitons.

It is thus important to reexamine the parameters of the experiment described in Ref. burger . According to the estimates obtained in Appendix B, the length of the trap is significantly larger than its radius. Hence the approximation by an infinite cylindrical trap may be reasonable. Furthermore, a conservative estimate of the effective coupling leads to $`\gamma 7`$ which lies in the region $`\gamma >\gamma _1`$. Therefore, all three types of solitary waves theoretically discussed in the present paper should have been produced in the original experiment of Ref. burger . Also note that non-stationary vortex rings were actually detected in the experiment of Ref. anderson  but on a spherical trap that is significantly different than the cylindrical trap considered here. Nevertheless, our present analysis can be extended to the case of a finite trap, as discussed in Appendix B.

Finally, the calculation of the present paper clearly demonstrates that a Lieb-type mode appears in the excitation spectrum for all densities, either in the form of a quasi-1D soliton (for $`\gamma <\gamma _0`$) or in the form of a solitonic vortex (for $`\gamma >\gamma _0`$). It is thus interesting to examine whether or not such a mode can be measured.

###### Acknowledgements.

We are grateful to J. Brand for a stimulating discussion that brought the work of Ref. brand1  to our attention, to N.R. Cooper for his constant interest in the present effort, and to A.S. Fokas, and N.S. Manton for a number of useful comments. This work was supported by the EPSRC under contract GR/R96026/01 (S.K.).

## Appendix A ASYMPTOTIC PHASE

The essence of the analytical argument given below was suggested to us by Fokas fokas . For a solitary wave traveling with constant velocity $`v`$, the continuity equation (4) may be written as $`vn/\xi +\mathbf{}(n\mathbf{}\varphi )=0.`$ In the limits $`\xi \pm \mathrm{},nn_0(\rho )=|\mathrm{\Psi }_0(\rho )|^2`$ and all $`\xi `$ (or $`z`$) derivatives vanish on account of the boundary conditions (9). Hence, $`\mathbf{}_{}[n_0(\rho )\mathbf{}_{}\varphi _\pm (x,y)]=0`$, where $`\mathbf{}_{}`$ is the 2D gradient operator and $`\varphi _\pm (x,y)=\varphi (x,y,\xi =\pm \mathrm{})`$ are the unknown asymptotic limits of the phase. If we further employ cylindrical coordinates ($`x=\rho \mathrm{cos}\phi ,y=\rho \mathrm{sin}\phi `$) the asymptotic form of the continuity equation reads

$$\rho \frac{}{\rho }\left[\rho n_0(\rho )\frac{\varphi _\pm }{\rho }\right]+n_0(\rho )\frac{^2\varphi _\pm }{\phi ^2}=0.$$
(15)
For axisymmetric solitons and vortex rings $`\varphi _\pm =\varphi _\pm (\rho )`$ and thus the second term in Eq. (15) vanishes. A first integral is then given by $`d\varphi _\pm /d\rho =c_1/\rho n_0(\rho )`$ where $`c_1`$ is some integration constant. Also note that the energy functional (2) contains a term $`n(\varphi /\rho )^2`$ that becomes $`n_0(\rho )(d\varphi _\pm /d\rho )^2=c_1^2/\rho ^2n_0(\rho )`$ at the two ends of the trap. Since the ground-state density $`n_0(\rho )`$ is finite at $`\rho =0`$ and vanishes in the limit $`\rho \mathrm{}`$, the above term would lead to a badly divergent integral over the radial distance $`\rho `$ unless the integration constant $`c_1`$ vanish. Hence, for an axisymmetric solitary wave with finite excitation energy, $`d\varphi _\pm /d\rho =0`$ and thus the asymptotic phases $`\varphi _\pm `$ can depend only on the velocity $`v`$.

On the other hand, we have not been able to extend the above argument to the case of non-axisymmetric solitary waves for which both terms in Eq. (15) are potentially relevant. Nevertheless, our numerical results for the solitonic vortex discussed in the main text again suggest that the asymptotic phases $`\varphi _\pm =\varphi _\pm (x,y)`$ are, in fact, independent of $`x`$ and $`y`$.

## Appendix B FINITE TRAPS

A finite harmonic trap with axial symmetry is characterized by the two confinement frequencies $`\omega _{}`$ and $`\omega _{}`$. We now introduce rationalized quantities that are in some respects different than those employed for the discussion of the infinite cylindrical trap. Time and distance are still measured in units of $`1/\omega _{}`$ and $`a_{}`$ but the wave function is rescaled according to $`\mathrm{\Psi }\sqrt{N}\mathrm{\Psi }/a_{}^{3/2}`$, where $`N`$ is the total number of atoms, and thus satisfies the constraint $`|\mathrm{\Psi }|^2𝑑V=1`$. The Gross-Pitaevskii equation is then formally identical to Eq. (1) but

$$g=4\pi \alpha ,V_{\mathrm{tr}}=\frac{1}{2}(\rho ^2+\beta ^2z^2),$$
(16)
where $`\alpha =Na/a_{}`$ and $`\beta =\omega _{}/\omega _{}`$ are now the two independent dimensionless parameters. To make contact with the infinite trap, we calculate the average linear density from $`\nu (z)=(N/a_{})|\mathrm{\Psi }|^2𝑑x𝑑y`$ and define the dimensionless quantity $`\gamma (z)=a\nu (z)`$ or

$$\gamma (z)=\alpha |\mathrm{\Psi }|^2𝑑x𝑑y,$$
(17)
which is the analog of the dimensionless coupling constant $`\gamma `$ employed in the main text, except that it now depends on $`z`$. Nevertheless, $`\gamma (z)`$ becomes increasingly uniform in the limit $`\alpha \mathrm{}`$ and $`\beta 0`$, holding $`\gamma (0)\gamma `$ fixed, which is the ideal limit of an infinite cylindrical trap studied in the main body of the paper.

As a simple illustration, we consider the Thomas-Fermi (TF) approximation of the ground state baym  to find that

$$\gamma (z)=\frac{1}{16}(2\mu \beta ^2z^2)^2,$$
(18)
where the chemical potential is given by $`2\mu =(15\alpha \beta )^{2/5}`$. Therefore, the estimated maximum value of $`\gamma `$ is given by $`\gamma _{\mathrm{max}}=\gamma (0)=\mu ^2/4`$. A more conservative estimate is the average value defined from $`\gamma _{\mathrm{av}}=\alpha /L`$ where $`L=2\sqrt{2\mu }/\beta `$ is the TF length of the trap and thus $`\gamma _{\mathrm{av}}=(8/15)\gamma _{\mathrm{max}}`$. The TF radius of the trap is accordingly given by $`R=\sqrt{2\mu }`$. For the parameters of the experiment described in Ref. burger  we find $`R=2\mu \text{m}`$ and $`L=115\mu \text{m}`$, while $`\gamma _{\mathrm{max}}=13`$ and $`\gamma _{\mathrm{av}}=7`$. These estimates are further discussed in the concluding section in the main text.

Now, static solitons, solitonic vortices, and vortex rings exist also in a finite trap and we have calculated them explicitly for various values of $`\alpha `$ and $`\beta `$ to confirm the general picture described in Sec. III. Here, we briefly summarize the basic facts. For each value of the aspect ratio $`\beta `$, there exist two critical couplings $`\alpha _0=\alpha _0(\beta )`$ and $`\alpha _1=\alpha _1(\beta )`$ that lead to a bifurcation pattern similar to that shown in Fig. 1. The critical coupling $`\alpha _0`$ coincides with the coupling where a purely imaginary eigenvalue first appears in the $`m=1`$ sector of the linear BdG equations muryshev1 . Similarly, the critical coupling $`\alpha _1`$ is related to an instability that occurs in the $`m=0`$ sector, but the picture is now significantly complicated by the appearance of complex modes feder . Our numerical solution of the BdG equations shows that there exists an intermediate critical coupling $`\alpha _1^{}<\alpha _1`$ where two real eigenvalues with $`m=0`$ coalesce and eventually become complex for $`\alpha >\alpha _1^{}`$. The precise behavior of the complex eigenvalues with increasing $`\alpha `$ depends on the specific value of $`\beta `$. A common feature is that complex eigenvalues as such again disappear before reaching the critical coupling $`\alpha _1`$. The latter is characterized by the fact that only one purely imaginary eigenvalue persists for $`\alpha >\alpha _1`$ and signals a new bifurcation that leads to the appearance of axisymmetric vortex rings.

As an example we consider the popular ratio $`\beta =1/4`$ for which we find $`\alpha _0=10,\alpha _1^{}=32,\alpha _1=39`$. A static soliton exists for all $`\alpha `$ but is stable only for $`\alpha <10`$, while a static solitonic vortex appears as the most stable structure for $`\alpha >10`$. A new threshold occurs at $`\alpha =39`$ where a static vortex ring emerges with energy intermediate between that of a soliton and a solitonic vortex, in analogy with the situation on an infinitely cylindrical trap demonstrated in Fig. 1. Explicit illustrations of the currently calculated profiles are given in Figs. 8,9 for $`\alpha =100`$ and $`\beta =1/4`$ where all three types of static solitary waves are possible. For smaller $`\beta `$ the critical couplings $`\alpha _0`$ and $`\alpha _1`$ increase indefinitely. But the corresponding effective couplings $`\gamma _0`$ and $`\gamma _1`$, numerically calculated from Eq. (17) applied for $`z=0`$, converge to the asymptotic values $`\gamma _0=1.5`$ and $`\gamma _1=4`$ that are in agreement with the critical couplings obtained in the main text for an infinite cylindrical trap. It is also interesting to consider the special limit of a spherical ($`\beta =1`$) trap in view of the experiment described in Ref. anderson . The critical coupling $`\alpha _0`$ vanishes at $`\beta =1`$ and thus the black soliton is never stable on a spherical trap. On the other hand, the solitonic vortex now becomes degenerate with the ordinary vortex and is stable for all $`\alpha >0`$. The vortex-ring threshold is here predicted to occur at $`\alpha _1=5`$. For the parameters of the actual experiment anderson  we find $`\alpha =430`$ which is substantially larger than $`\alpha _1`$ and eventually explains the observation of multiple vortex rings.