text:
# Stabilization of solitons of the multidimensional nonlinear Schrödinger equation: Matter-wave breathers

## 1 Introduction

The nonlinear Schrödinger equation (NLSE) in its many versions is one of the most important models of mathematical physics, with applications to different fields such as plasma physics, nonlinear optics, water waves and biomolecular dynamics, to cite only a few cases. In many of those examples the equation appears as an asymptotic limit for a slowly varying dispersive wave envelope propagating in a nonlinear medium .

A new burst of interest on problems modelized by nonlinear Schrödinger equations has been triggered by the achievement of Bose-Einstein condensation using ultracold neutral bosonic gases .

In particular in Bose-Einstein condensation one of the key ingredients of the achievement of condensation with alkali gases is the trapping of the neutral atoms which are cooled down below the transition temperatures. Although different trapping techniques have been used in practice, the most commonly used are the magnetic traps which are modelized by external confining forces acting on the system.

There are specific types of Bose-Einstein condensates such as those made of Lithium  in which the interactions between the atoms are attractive. Mathematically this implies that for spatial dimensionalities $`n=2,3`$ collapsing solutions are possible as it is well known in the framework of the usual studies of nonlinear Schrödinger equations . Thus, although in principle one could think that nonlinearity and dispersion could balance mutually and a self-sustained solitonic structure might exist this is only true for $`n=1`$. In fact, it was soon found experimentally  and theoretically supported that such a situation is unstable and leads to collapse . Later studies, taking into account more elaborate models than the simpler NLS mean field models, led to the understanding that the occurrence of collapse during the condensation process would limit the size of an attractive condensate .

Thus, the only confirmed way to observe solitonic states in Bose-Einstein condensates (matter-wave solitons) with negative scattering length involves the elimination of the trap in one direction as proposed in  and found experimentally in .

However, the possibility of using Feschbach resonances to control the scattering length  has provided a way to study large negative scattering length condensates and collapse processes in detail . This control allows, with certain experimental limitations, to modulate appropriately the nonlinear term.

But this possibility opens new ways for the generation and observation of different types of matter-wave solitons not yet fully explored. One of the most striking ones is the possibility of generating *trapless trapped Bose-Einstein condensate solitons*. The idea is that (oscillating) bound states might be obtained by combining cicles of positive and negative scattering length values so that after an expansion and contraction regime the condensate would come back to the initial state. In this way some sort of pulsating trapped condensate, i.e., a *breather*, would be obtained.

A rigorous study of this intuitive idea is necessary in order to make precise predictions and find which are the precise parameter values to be used to obtain such breathers. This is the purpouse of this paper in a general framework of nonlinear Schrödinger systems. This idea has been explored in two previous papers  but here we improve the understanding of the phenomenon and correct mistakes contained in the previous analysis.

Another field of applications of these ideas is nonlinear Optics, there a full control of the nonlinear term is not possible but, because of technical limitations, only piecewise constant values for the nonlinear coefficient can be easily generated experimentally.

¿From the mathematical point of view what we would like is to stabilize solutions close to the stationary unstable solitons of the cubic NLSE by choosing appropriate controls of the nonlinear term. To our knowledge this problem has not been considered previously in the mathematical literature.

This paper is organized as follows: First, in Sec. 2 we present the model equations as they appear in one specific application of the model. Next, in Sec. 3 we reduce the nonlinear model to a system of ordinary differential equations by using the so-called moment method. In Sec. 4 we analyze two-dimensional systems and compare the analytical predictions of the moment method with numerical simulations of the full partial differential equations. We obtain precise conditions for the existence of breathers and discuss the importance of the initial data. In Sec. 5 we analyze the three-dimensional case and explain why moment (or related) equations fail to describe the dynamics of the spherically symmetric case. We consider the case of systems in specific external potentials and give some stabilization conditions. Finally, in Sec. 6 we summarize our conclusions and compare our results with the previous findings of Refs. .

## 2 The model

In this work we study systems ruled by the NLSE with a cubic nonlinearity. In the framework of recent problems in Bose-Einstein condensation this equation appears as the model for the mean field dynamics of a boson system described by a single wavefunction $`\mathrm{\Psi }`$ in the zero–temperature limit. The resulting NLSE equation is called sometimes the Gross–Pitaevskii equation (GPE) which is

$$i\mathrm{}\frac{\mathrm{\Psi }}{\tau }=\left[\frac{\mathrm{}^2}{2m}\mathrm{}+V(𝐫)+\frac{4\pi a\mathrm{}^2}{m}|\mathrm{\Psi }|^2\right]\mathrm{\Psi },$$
(1)
where $`\mathrm{}=_{j=1}^3^2/r_j^2`$ is the three-dimensional Laplacian, $`V(𝐫)=\frac{1}{2}m\omega ^2_{j=1}^3\lambda _j^2r_j^2`$ is the external potential (the so called “trap”) which confines the condensate and $`a`$ is the $`s`$–wave scattering length for the binary collisions within the condensate.

It is more convenient to work with a new set of adimensional quantities defined as $`x_j=r_j/a_0,t=\tau /T`$, $`\psi (𝐱,t)\mathrm{\Psi }(𝐫,\tau )\sqrt{a_0^3/N}`$, where $`a_0=\sqrt{\mathrm{}/m\omega }`$, $`T=1/\omega `$ and $`N=|\mathrm{\Psi }|^2d^3r`$ is the number of particles in the condensate. Then Eq. (1) reads as the next NLSE

$$i\frac{\psi }{t}=\left[\frac{1}{2}\mathrm{}+\frac{1}{2}\underset{j=1}{\overset{3}{}}\lambda _j^2x_j^2+g|\psi |^2\right]\psi ,$$
(2)
where $`g=4\pi Na/a_0`$ and $`|\psi |^2d^3x=1`$. In this paper we consider a system evolving without external potential along one or more dimensions. In Bose-Einstein condensation this corresponds to a condensate generated in a trap and then the trapping is eliminated. This leads to different particular cases of Eq. (2). When the potential is removed in all directions we obtain

$$i\frac{\psi }{t}=\left[\frac{1}{2}\mathrm{}+g|\psi |^2\right]\psi .$$
(3)
This equation holds for three-dimensional and quasi two-dimensional systems as will be discussed in detail in Sec. 5. Eq. (3) has the typical form of the NLSE with cubic power nonlinearity  and has been extensively studied. In this work we will analyze the possibilities of controlling the behavior of the stationary solutions of Eq. (3) by letting $`g`$ be a time-dependent real function. Similar possibilities for nonconstant $`g`$ (in that context z $``$ t) but restricted to $`n=2`$ arise in the context of coherent light propagation in nonlinear Kerr media in the paraxial approximation.

## 3 Moment method

In this section we build a theory which will allow us to reduce the dynamics from the complex model given by Eq. (3) to a simpler one whose analysis will provide insight on the basic features of the phenomenon.

Let us first consider radially symmetric solutions of Eq. (3) $`u=u(r,t)`$, $`r=(_{j=1}^nx_j)^{1/2}`$, satisfying

$$i\frac{u}{t}=\frac{1}{2r^{n1}}\frac{}{r}\left(r^{n1}\frac{u}{r}\right)+g(t)|u|^2u.$$
(4)
where $`n`$ is the spatial dimensionality, in our case $`n=2,3`$. To get information on the solutions of Eq. (4) we will use the moment method . This method proceeds by analyzing the evolution of several integral quantities  related with Eq. (4).

$`I_1(t)`$ $`=`$ $`{\displaystyle |u|^2d^nx}`$ (5a)
$`I_2(t)`$ $`=`$ $`{\displaystyle |u|^2r^2d^nx},`$ (5b)
$`I_3(t)`$ $`=`$ $`i{\displaystyle (u\frac{u^{}}{r}u^{}\frac{u}{r})rd^nx},`$ (5c)
$`I_4(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \left(\left|u\right|^2+\frac{n}{2}g(t)|u|^4\right)d^nx},`$ (5d)
$`I_5(t)`$ $`=`$ $`{\displaystyle \frac{n}{4}}{\displaystyle |u|^4d^nx},`$ (5e)
where $`d^nx=r^{n1}drd\mathrm{\Omega }`$. With our scaling for $`\psi `$, the first one satisfies $`I_1=1`$, for all $`t`$. The remaining ones are related physically to the width, radial momentum and energy of the wave packet. In what follows we assume that the initial data are such that all $`I_j`$ are initially well-defined . Some algebra leads to

$`{\displaystyle \frac{dI_2}{dt}}`$ $`=`$ $`I_3,`$ (6a)
$`{\displaystyle \frac{dI_3}{dt}}`$ $`=`$ $`4I_4,`$ (6b)
$`{\displaystyle \frac{dI_4}{dt}}`$ $`=`$ $`g{\displaystyle \frac{n2}{n}}{\displaystyle \frac{dI_5}{dt}}+{\displaystyle \frac{dg}{dt}}I_5,`$ (6c)
$`{\displaystyle \frac{dI_5}{dt}}`$ $`=`$ $`{\displaystyle \frac{n\pi 2^n}{8}}{\displaystyle \frac{|u|^4}{r}\frac{\mathrm{arg}u}{r}r^{n1}𝑑r}{\displaystyle \frac{nI_3I_5}{2I_2}}.`$ (6d)
All the momenta follow closed evolution laws except for $`I_5(t)`$. In specific situations the moment equations provide completely closed equations for the evolution of all the $`I_j`$  but this is not our case. Extending the number of momenta which are included in the calculation does not help to obtain a closed set of equations.

Thus, to close the system we will take

$$\mathrm{arg}u=I_3r^2/4I_2,$$
(7)
as used in . Physically, this corresponds to approximating the phase of $`u`$ by the spherical wavefront which best fits the distribution. A rigorous justification of this choice is possible  for the case $`n=2`$ and when the initial data is the so-called Townes soliton to be presented later.

When solutions with phase given by Eq. (7) are considered Eqs. (6) are closed and have several (positive) invariants under time evolution given by

$`Q_1`$ $`=`$ $`2(I_4gI_5)I_2I_3^2/4,`$ (8a)
$`Q_2`$ $`=`$ $`2I_2^{n/2}I_5.`$ (8b)
With the help of these quantities, Eqs. (6) can be reduced to a single equation for $`I_2(t)`$, which is

$$\frac{d^2I_2}{dt^2}\frac{1}{2I_2}\left(\frac{dI_2}{dt}\right)^2=2\left(\frac{Q_1}{I_2}+g\frac{Q_2}{I_2^{n/2}}\right).$$
(9)
If we were able to solve Eq. (9) then the use of Eqs. (6) would allow us to track the evolution of the other ones. To simplify Eq. (9) we define $`X(t)=\sqrt{I_2(t)}`$, which is the wave packet width, and substituting it into (9) we get

$$\ddot{X}=\frac{Q_1}{X^3}+g(t)\frac{Q_2}{X^{n+1}}.$$
(10)
## 4 Two-dimensional systems

### 4.1 Rigorous analysis of the moment equations

In physical situations both problems with $`n=2`$ and $`n=3`$ arise. The situation when a Bose-Einstein condensate is tightly confined along a particular direction leads to a quasi two-dimensional system such as the ones obtained in Ref. . From now on we consider the two-dimensional model to be valid in that situation although a more detailed analysis will be made in Sec. 6. Thus, when $`n=2`$ we obtain

$$\ddot{X}=\frac{Q_1+Q_2g(t)}{X^3}\frac{p(t)}{X^3}.$$
(11)
This equation appears also in relation with the classical motion of an atom near an infinite straight wire with oscillating charge (the Paul trap)  and also as an approximate model for beam propagation in layered Kerr media .

Let us consider the case where $`p`$ is continuous and $`T`$-periodic and look for positive periodic solutions and bound states, that is, bounded (nonperiodic) positive solutions without collapse. If $`p(t)0`$ all solutions escape to infinity whereas if $`p(t)<0`$ all solutions collapse. In fact, theorem 2.1 of  implies that a necessary condition for existence of bound states is

$$\overline{p}=\frac{1}{T}_0^Tp<0.$$
(12)
Let $`p(t)`$ be parametrized as $`p(t)=\alpha +\beta a(t)`$ with $`\overline{a}=0`$. We can fix $`a_M=\mathrm{max}a(t)=1`$ without loss of generality. A direct consequence of the latter result is that if $`\alpha 0`$ bound states cannot exist, therefore $`\alpha `$ must be negative. Also, $`p(t)`$ must change sign, so for existence of bound states it is necessary that

$$\alpha +|\beta |>0.$$
(13)
Let us introduce the positive small parameter $`ϵ=2\alpha /3|\beta |`$. We will now prove that *there exists $`ϵ_0>0`$ (depending only on $`a`$) such that Eq. $`(\text{11})`$ has a Lyapunov-stable periodic solution if $`0<ϵ<ϵ_0`$*. Moreover, *there is an infinite number of quasiperiodic solutions and a sequence of subharmonics with minimal periods tending to $`+\mathrm{}`$*. To prove this affirmation let us note that when $`a(t)=\mathrm{cos}\mathrm{\Omega }t`$, Theorem 1 in  asserts that if $`ϵ`$ is small enough there is a periodic solution of twist type . However, the proof is still valid for a general continuous and $`T`$-periodic function with $`\overline{a}=0`$. Moser Twist Theorem  implies that a solution of twist type is Lyapunov-stable as a consequence of the existence of invariant curves (quasiperiodic solutions) around it. Finally, the existence of subharmonics of any order is a consequence of Poincare-Birkhoff Theorem.

This result has interesting physical consequences. The existence of periodic solutions is directly related to the existence of pulsating breathers. The bounded solutions of the theorem would correspond to solutions with quasiperiodic or chaotic oscillations in the width. In fact, Mather sets and Smale’s horseshoes appear in the Poincaré map. The latter correspond to chaotic oscillations whereas Mather sets are Cantorian structures that are more complicated than the quasiperiodic ones.

We have just seen that a necessary condition for the existence of bounded solutions is that $`\alpha +|\beta |>0`$ with $`\alpha <0`$. Moreover the previous result indicates that the amplitude of the oscillating term must be of the same order as the amplitude of the non-oscillating term.

To find particular types of these breathers let us take $`g(t)`$ as $`g(t)=g_0+g_1\mathrm{cos}(\mathrm{\Omega }t)`$ so that $`p(t)=\alpha +\beta \mathrm{cos}(\mathrm{\Omega }t)`$ with $`\alpha =Q_1+g_0Q_2`$ and $`\beta =g_1Q_2`$. An useful tool to get periodic solutions is the stability diagram of Eq. (11). We look for the region of parameters for which there exist stable periodic solutions. To do this we reduce the number of parameters by defining a new time as $`\stackrel{~}{t}=\mathrm{\Omega }t`$. Then Eq. (11) is

$$X^{\prime \prime }=\frac{\stackrel{~}{\alpha }+\stackrel{~}{\beta }\mathrm{cos}(\stackrel{~}{t})}{X^3},$$
(14)
where the prime denotes derivative with respect to $`\stackrel{~}{t}`$, $`\stackrel{~}{\alpha }=\alpha /\mathrm{\Omega }^2`$ and $`\stackrel{~}{\beta }=\beta /\mathrm{\Omega }^2`$. Note that we use the word stable with two different meanings: the more mathematical one to indicate Lyapunov-stability and the more physical one to indicate existence of periodic solutions of the nonlinear Schrödinger equation.

Our objective is to find the region of parameters ($`\stackrel{~}{\alpha }`$,$`\stackrel{~}{\beta }`$) in which there exist initial conditions $`X_0=X(\stackrel{~}{t}=0)`$, $`X_0^{}=X^{}(\stackrel{~}{t}=0)=0`$ for which Eq. (14) has a stable periodic solution. We already know that our search can be restricted to the region given by the conditions $`\stackrel{~}{\alpha }<0`$ and $`\stackrel{~}{\alpha }>|\stackrel{~}{\beta }|`$. To calculate the stabilization region we first solve numerically Eq. (14). Then we take into account the results that finding a periodic solution of Eq. (14) is equivalent to finding a solution which verifies the Neumann conditions $`X^{}(0)=0=X^{}(2\pi )`$. The reason is that by extending in an even way this solution, we obtain a 4$`\pi `$-periodic solution that, of course, is bounded and periodic and this implies that there is a $`2\pi `$-periodic solution . To obtain solutions that verify the Neumann conditions, given a pair of parameters ($`\stackrel{~}{\alpha }`$,$`\stackrel{~}{\beta }`$), we look for two initial conditions $`X_{1,0}`$ and $`X_{2,0}`$ for which $`X_1^{}(2\pi )X_2^{}(2\pi )<0`$. Then, by continuity with respect to the initial conditions it is known that there exists one initial condition $`X_0`$, between $`X_{1,0}`$ and $`X_{2,0}`$, for which the solution of the equation satisfies $`X^{}(2\pi )=0`$. We will use two more observations. First of them is that if $`X(\stackrel{~}{t})`$ is a solution of the equation with parameters ($`\stackrel{~}{\alpha }`$,$`\stackrel{~}{\beta }`$), then $`Y(\stackrel{~}{t})=cX(\stackrel{~}{t})`$ with $`c`$ any arbitrary positive constant is also a solution with parameters ($`\stackrel{~}{\alpha }c^4`$,$`\stackrel{~}{\beta }c^4`$). Therefore, by moving $`c`$, all the points of the line $`y=(\stackrel{~}{\beta }/\stackrel{~}{\alpha })x`$ are obtained and this is nothing but the line which passes through the points ($`\stackrel{~}{\alpha }`$,$`\stackrel{~}{\beta }`$) and $`(0,0)`$. So, to scan the region of parameters we can, for example, to pay attention to an specific value of $`\stackrel{~}{\alpha }`$ and then change continuosly the parameter $`\stackrel{~}{\beta }`$. The second observation is that if $`X(\stackrel{~}{t})`$ is a periodic solution corresponding to a specific choice of ($`\stackrel{~}{\alpha }`$,$`\stackrel{~}{\beta }`$) then $`Y(\stackrel{~}{t})=X(\stackrel{~}{t}+\pi )`$ is also a solution for parameters ($`\stackrel{~}{\alpha }`$,-$`\stackrel{~}{\beta }`$) since $`\mathrm{cos}(\stackrel{~}{t}+\pi )=\mathrm{cos}(\stackrel{~}{t})`$. Therefore, the stability region is symmetric and we can restrict the search to $`\stackrel{~}{\beta }>0`$. With all these considerations the results obtained show that the region of parameters for which there exist initial conditions leading to stable periodic solutions is given by Eq. (13) so this condition seems to be not only necessary but also a sufficient one. In Fig. 1 this stability region is plotted.

### 4.2 Comparison with full numerical simulations of the NLSE

In this subsection we will compare the predictions based on the ordinary differential equation (11) with the simulations of the full nonlinear Schrödinger equation (3).

First of all we have to choose appropriate initial data $`u_0=u(r,t=0)`$ in order to solve numerically the NLSE. Since we are going to try to stabilize the solutions it is reasonable to start trying to stabilize the so-called stationary solutions. So, let us consider stationary solutions of Eq. (3) when $`g=g_s`$ is constant, given by $`u(\text{x},t)=e^{i\mu t}\mathrm{\Phi }_\mu (\text{x})`$. Here $`\mathrm{\Phi }_\mu (\text{x})`$ verifies

$$\mathrm{}\mathrm{\Phi }_\mu 2\mu \mathrm{\Phi }_\mu 2g_s|\mathrm{\Phi }_\mu |^2\mathrm{\Phi }_\mu =0.$$
(15)
As it is precisely stated in , when $`g_s`$ is negative, for each positive $`\mu `$ there exists only one solution of Eq. (1) which is real, positive and radially symmetric and for which $`(|\mathrm{\Phi }_\mu |^2d^nx)^{1/2}`$ has the minimum value between all of the possible solutions of Eq. (15). Moreover, the positivity of $`\mu `$ ensures that this solution decays exponentially at infinity. This solution is called the ground state and, in two dimensions, is known as the Townes soliton. We will denote this solution as $`R_\mu (r)`$ and satisfies the equation

$$\mathrm{}R_\mu 2\mu R_\mu 2g_sR_\mu ^3=0$$
(16)
and the boundary conditions

$$\underset{r\mathrm{}}{lim}R_\mu (r)=0,R_\mu ^{}(0)=0.$$
(17)
Fixed a value of $`\mu `$, the norm and width of $`R_\mu `$ are given by $`\eta _\mu =(|R_\mu |^2d^2x)^{1/2}`$ and $`X_\mu =(|R_\mu |^2r^2d^2x)^{1/2}`$ respectively. Applying scaling trasformations to $`R_\mu (r)`$ it is possible to build stationary solutions having the same shape and norm as the Townes soliton but different widths

$$R_\mu (r)=\mu ^{1/2}R_1(\mu ^{1/2}r),X_\mu =X_1/\mu ^{1/2}.$$
(18)
The equation which verifies the normalized soliton $`R_{\mu ,N}(r)=R_\mu (r)/\eta _\mu `$ is

$$\mathrm{}R_{\mu ,N}2\mu R_{\mu ,N}2g_s\eta _\mu ^2R_{\mu ,N}^3=0.$$
(19)
We have taken as initial condition $`u(r,t=0)=R_{\mu ,N}(r)`$ for some $`\mu `$ (therefore the initial width is fixed) and $`g_s`$ values. To obtain the shape $`R_\mu (r)`$ we have used a shooting method  to solve Eq. (16). The idea is to rewrite Eq. (16) as a dynamical system and impose that the solution of such a system also verifies the boundary conditions (17). In order to avoid the singularity which appears at the point $`r=0`$ we have solved the indetermination by means of a series expansion of $`R_\mu `$ around $`r=0`$.

¿From the theory of nonlinear Schrödinger equations it is known that the Townes soliton is unstable, i.e. small perturbations of this solution lead to either expansion of the initial data or blow-up in finite time. Thus, following the analogy with the stabilization of unstable fixed points in finite dimensional dynamical systems by fast variation of a parameter, we will try to stabilize this unstable but anyway stationary solution of Eq. (3).

Using the predictions of Eq. (11), in order to get stable periodic solutions, we have to choose $`g_0`$ such that $`\alpha =Q_1+g_0Q_2`$ is negative, that is, the following relation must hold

$$g_0<\frac{|u_0|^2d^2x}{|u_0|^4d^2x}.$$
(20)
When only such a $`g_0`$-term is present the solutions collapse as predicted by Eq. (11). So, the $`g_1`$-term is necessary in order to arrest the collapse and achieve the stabilization of the soliton given by $`u_0`$.

Thus, we have first taken as initial data the Townes soliton obtained for $`g_s=0.5`$ and $`\mu =0.5`$ whose width is $`X_0=1.09`$ and its norm satisfies $`g_s\eta _{\mu =0.5}^2=5.85`$. This value is precisely the critical nonlinearity for collapse below which the collapse occurs. Let us note that only for the Townes soliton it is satisfied that the collapse threshold in Eq. (3), when $`g`$ is a constant, is just the value of the nonlinearity in Eq. (19) .

In all the simulations of Eq. (3) we have made there appear two different types of behaviors. First, there is a small-amplitude fast oscillation due to the term $`\mathrm{cos}(\mathrm{\Omega }t)`$ which is present in both the full model Eq. (3) and the reduced ODE system. Secondly, there is a “slow” dynamics due to the internal dynamics of the system which sometimes does not appear in the ODE model. When the results from the ODE (11) and the full PDE are compared qualitatively similar dynamics are observed in many situations. For instance, in Fig. 2 we plot the results of the simulation for parameter values $`g_0=2\pi `$, $`g_1=8\pi `$ and $`\mathrm{\Omega }=40`$. To solve numerically Eq. (3) we have used a pseudo-spectral method combined with a second order split-step method to advance in time . We have also compared this simulation with the one obtained with a fourth order split-step method and found both results in full agreement.

It is important to indicate that it is absolutely necessary to incorporate an absorbing potential at the frontier of the simulation region to avoid possible interferences between the fraction of the wave packet moving outside (i.e. a possible non-trapped fraction of initial data) and the fraction of the initial data which remains trapped. Typical grid sizes were $`256\times 256`$ and the time step was $`dt=0.001`$. All results were tested on different grid sizes and changing time steps.

Clearly, the system is trapped with a fast modulation in the nonlinearity ($`g_1`$-term) and the collapse process is inhibited. Also, there is a fast oscillation with the same frequency $`\mathrm{\Omega }`$ as the one of the oscillating term, (beyond the resolution of the plot), and a slow oscillation due to the proper dynamics of the system.

In Fig. 2(e) it is observed that the norm of the solution decreases in time. This is caused by the action of the absorbing potential on the outgoing way.

The next step in our study is to compare this stabilization in the NLSE with the differential equation for the width. In Fig. 3 we see that the quantitative predictions of Eq. (11) are not very precise as what concerns the amplitude and frequency of the slow oscillation. In any case, at least the simple ODE model predicts correctly the trapping of the solution for these parameters. These results imply that the predictions of the ODE system must be taken only as qualitative indications of the possible dynamics.

Also, there exist parameter ($`g_0`$, $`g_1`$ and $`\mathrm{\Omega }`$) choices which lead to stabilization in Eq. (11) but not in Eq. (3). For example, according to Eq. (11) and the diagram of Fig. 1 the parameters $`g_0=6.4`$, $`g_1=33.4`$ and $`\mathrm{\Omega }=20`$ stabilize the initial soliton as shown in Fig. 4. However, this stabilization does not happen when we solve numerically Eq. (3), instead, after a few oscillations their dynamics are drastically different.

### 4.3 Stabilization of other initial data

Although a Townes soliton is a natural object to stabilize in the framework of Eq. (3) we may wonder (i) if the procedure described above does stabilize other type of initial data and (ii) if the stabilization is appropriately described by Eqs. (11).

From the point of view of applications of Bose-Einstein condensation the possibility of stabilizing other initial data is important since it is not clear how a Townes soliton could be generated in real experiments. On the other hand, other initial data such as Thomas-Fermi type solutions or gaussians are much more natural and easier to obtain. This, together with the fact that the usual time dependent variational method  is usually developed for gaussian profile functions has lead to some interest on the possibility of trapping gaussian initial data (which is the case studied in Refs. ).

Therefore, let us take $`u_0`$ as

$$u_0(r)=\frac{e^{r^2/2X_0^2}}{\sqrt{\pi }X_0},$$
(21)
for which

$`\eta `$ $`=`$ $`\left({\displaystyle |u_0|^2d^2x}\right)^{1/2}=1,`$ (22a)
$`X_0`$ $`=`$ $`\left({\displaystyle |u_0|^2r^2d^2x}\right)^{1/2}.`$ (22b)
In this case we get for the invariants the values $`Q_1=1`$ and $`Q_2=1/2\pi `$ independently of the width $`X_0`$. So, the equation for the width is

$$\ddot{X}=\frac{1+g(t)/2\pi }{X^3}\frac{p(t)}{X^3}.$$
(23)
which corresponds to the evolution equation of Refs. . When $`g`$ is a constant this equation predicts a threshold value for collapse of $`g=2\pi `$ which is, indeed, below the real value which corresponds to the Townes soliton ($`g_s\eta _\mu ^2=5.85`$). This means that when a Gaussian function is taken as initial data, Eq. (23) does not describe accurately the region of trapping, because for $`2\pi <g_0<5.85`$ this equation predicts expansion of the initial data while the real dynamics of the partial differential equation is collapsing.

In Fig. 5 the evolution of a Gaussian with $`X_0=1.09`$ and parameter values $`g_0=2\pi `$, $`g_1=8\pi `$ and $`\mathrm{\Omega }=40`$ is shown. Notice that these are the same width and parameters which allow trapping of a Townes soliton.

We can see that with these parameters the stabilization is achieved although Eq. (23) does not predict trapping (in fact for these parameters $`\overline{p}=0`$). What is the difference with the stabilization of the Townes soliton shown in Fig. 2?. A first important observation is that in the case of Gaussian initial data a readjustment is produced soon by ejecting a significant part of the wave packet far from the trapping region. This corresponds to the first 10 time units where the width significantly increases due to the contribution of the outgoing wave. When this wave hits the absorbing region it is dissipated leading to the step in the norm evolution shown in Fig. 5(e). What remains trapped is indeed a Townes soliton.

Thus, the use of initial data different from a Townes soliton leads to a splitting of the solution into the soliton itself, which is the structure which can be stabilized, plus a certain amount of radiation which goes far from the region of interest. Because of this fact one must be careful to eliminate the outgoint part of the radiation in the numerical simulations. We have verified that when the absorbing region is absent the numerical simulations are misleading and the trapping effects are drastically altered (in fact, when zero boundary conditions at a given distance are used, there appear reflections in the boundary and essential changes on the results take place which would lead to spurious destabilization of the system after very short trapping times).

It is known that Eq. (3) with a cubic nonlinearity and two spatial dimensions corresponds to the so-called critical case  in which diffraction and self-focusing are nearly balanced and collapse is extremely sensitive to perturbations and to changes in the initial conditions. For this reason although Gaussians look roughly like the Townes soliton they are not able to capture the delicate balance between diffraction and nonlinear focusing present in this case. When a Townes soliton is taken as a initial condition in Eq. (3) there is only a very small outgoing component which spreads out (due to numerical errors) and the system responds to the parametric perturbation as a whole. This fact implies that the moment equations may describe reasonably well the dynamics of the system only near the Townes soliton, since they only deal the global dynamics of the system. However, when the initial condition is a Gaussian, the component which spreads out cannot be captured neither by a moment-type formalism nor by the usual variational methods.

## 5 Three-dimensional systems

### 5.1 Analysis of the moment equations

Let us now consider the case $`n=3`$, then Eq. (10) reads

$$\ddot{X}=\frac{Q_1}{X^3}+g(t)\frac{Q_2}{X^4}.$$
(24)
Let us first assume that $`g`$ is continuous and $`T`$-periodic, then if $`\overline{g}=\frac{1}{T}_0^Tg0`$ bound states cannot exist. The reason is that, by Massera theorem, a bound state would imply the existence of a periodic solution. But if there is a periodic solution $`X`$, multiplying $`(\text{24})`$ by $`X^4`$ and integrating over a period we get $`0>4_0^TX^3\dot{X}^2Q_1_0^TX=\overline{g}Q_2T0`$, which is a contradiction. Therefore, a necessary condition for the existence of bounded solutions is

$$\overline{g}=\frac{1}{T}_0^Tg<0.$$
(25)
Note that this result is the opposite to the one inferred in Ref.  from direct numerical simulations an equation similar to (but restricted to the class of gaussian initial data) that of Eq. (24).

Our goal is to find stable periodic solutions of Eq. (24) in three dimensions. A first useful observation is that an arbitrary positive $`T`$-periodic function $`X(t)`$ is a solution of $`(\text{24})`$ if $`g(t)=Q_2^1(X^4\ddot{X}Q_1X)`$. In the following result, we use the definition of positive part of a function as $`f(t)^+=\mathrm{max}\{0,f(t)\}`$.

###### Theorem 1

Let $`X(t)`$ be a $`T`$-periodic positive function such that

(i) $`{\displaystyle _0^T}4{\displaystyle \frac{\ddot{X}}{X}}{\displaystyle \frac{Q_1}{X^4}}>0`$

(ii) $`T{\displaystyle _0^T}{\displaystyle \frac{\ddot{X}^+}{X}}<1`$

Then, there exists $`g(t)`$ $`T`$-periodic function such that $`X(t)`$ is a $`T`$-periodic linear stable solution of $`(\text{24})`$.

Proof.  If we choose $`g(t)=Q_2^1(X^4\ddot{X}Q_1X)`$, then we know that $`X(t)`$ is a $`T`$-periodic solution of $`(\text{24})`$. The linearized equation is

$$\ddot{Y}+\left(4\frac{\ddot{X}}{X}\frac{Q_1}{X^4}\right)Y=0.$$

Now, hypotheses $`(i)(ii)`$ correspond to the assumptions of the classical Lyapunov’s criterion for linear stability (see for instance Theorem 1 in ).

In many cases, hypotheses $`(i)(ii)`$ can be verified numerically. For instance taking $`g(t)=Q_2^1(15+\mathrm{sin}t)^4\mathrm{sin}tQ_2^1Q_1(15+\mathrm{sin}t)`$, then equation $`(\text{24})`$ has the solution $`X(t)=15+\mathrm{sin}t`$ which is linearly stable.

A different way to build solutions to Eqs. (24) is to choose $`g(t)=G(t)X(t)`$ then Eq. (24) becomes equivalent to Eq. (11) and thus all the results of existence of breathers in two dimensions can be applied.

In conclusion, from the analysis of Eq. (24) we get that trapping could be possible in a three-dimensional scenario provided the dynamics of the system is well described by the moment equations (24).

### 5.2 Failure of moment equations for symmetric 3D systems

To check the validity of the predictions made on the basis of Eq. (24) we have compared its predictions with the behavior of the solution as obtained from direct numerical simulations of Eq. (3). We have used as initial data solutions of Eq. (19) for some $`\mu `$ and $`g_s`$ values obtained with a shooting method as we described above for the two-dimensional case. When $`n=3`$, the relations between soliton solutions of different widths are

$$R_\mu (r)=\mu ^{1/2}R_1(\mu ^{1/2}r),X_\mu =X_1/\mu ^{3/4}.$$
(26)
The norm is not conserved by the previous change and verifies the relation $`\eta _\mu =\eta _1/\mu ^{1/4}`$.

In all the simulations we have taken $`g(t)=g_0+g_1\mathrm{cos}(\mathrm{\Omega }t)`$ with parameters $`g_0,g_1`$ and $`\mathrm{\Omega }`$ chosen to give stable behavior on the basis of Eq. (24),

Our numerical results show that the stabilization predicted by Eq. (24) does not occur for solutions of Eq. (3). Although it is possible by a suitable choice of $`g_0`$, $`g_1`$ and $`\mathrm{\Omega }`$ to change the behavior from collapsing to expanding or viceversa from the behavior which corresponds to the $`g_1=0`$ case, we were not able to find stable breather solutions for any choice of the parameters. In the region in which the solutions change from collapsing to expanding it is possible to fine-tune the parameters to get a solution which holds for some time, but this structure is unstable and finally the solution loses its profile.

Thus, the moment equations fail to describe the dynamics of the system in three dimensions. The reason is simple: when the nonlinear term in Eq. (3) is cubic and $`n=3`$ we are in the so-called supercritical case in which collapse occurs by formation of a localized spike on a finite-amplitude background. This behavior cannot be captured neither with the quadratic phase approximation \[Eq. (7)\] nor with any type of variational ansatz.

### 5.3 Three dimensional systems with confinement along one direction

Although stabilization of spherical structures seems not possible in three-dimensional scenarios, we could think about the situation in which the three-dimensional system has cylindrical symmetry. It would seem plausible that if the system is coin-shaped stabilization might take place because, in fact, this system is close to a quasi two-dimensional one and, therefore, the results in two dimensions could be applied. In fact, a numerical simulation reported in Ref.  supports this conjecture. Here we will try to make a more systematic analysis of the phenomenon.

To get a qualitative understanding of the phenomenology for this situation let us consider Eq. (2). We will take gaussian functions as initial data and obtain the evolution equations for the parameters by the use of the collective coordinates method  to find the evolution of the widths in each spatial direction. The idea is to restate the problem of solving Eq. (2) as a variational problem, corresponding to a stationary point of the action related to the Lagrangian density $``$. So, the problem is transformed into the problem of finding $`\psi (\text{x},t)`$ such that the action

$$S=d^nx𝑑t$$
(27)
is extreme. This problem is as complicated as solving the original NLSE. The idea of the method is to restrict the analysis to a set of trial functions. One possible choice is to take a Gaussian ansatz of the form

$$\psi (𝐱,t)=\frac{1}{\pi ^{n/4}_{j=1}^nw_j(t)^{1/2}}\underset{j=1}{\overset{n}{}}e^{x_j^2/2w_j^2+ix_j\alpha _j(t)+ix_j^2\beta _j(t)},$$
(28)
which verifies the following relations

$`\eta (t)`$ $`=`$ $`\left({\displaystyle |u|^2d^nx}\right)^{1/2}=1,`$ (29a)
$`w_j(t)`$ $`=`$ $`2{\displaystyle |u|^2x_j^2d^nx}.`$ (29b)
The variational method leads to the following evolution equations for the widths

$$\ddot{w_j}+\lambda _j^2(t)w_j=\frac{1}{w_j^3}+\frac{g(t)/(2\pi )^{n/2}}{w_j}\left(\underset{k=1}{\overset{n}{}}\frac{1}{w_k}\right).$$
(30)
In the three-dimensional case the equations are

$`\ddot{w_1}+\lambda _1^2(t)w_1`$ $`=`$ $`{\displaystyle \frac{1}{w_1^3}}+{\displaystyle \frac{g(t)/(2\pi )^{3/2}}{w_1^2w_2w_3}},`$ (31a)
$`\ddot{w_2}+\lambda _2^2(t)w_2`$ $`=`$ $`{\displaystyle \frac{1}{w_2^3}}+{\displaystyle \frac{g(t)/(2\pi )^{3/2}}{w_1w_2^2w_3}},`$ (31b)
$`\ddot{w_3}+\lambda _3^2(t)w_3`$ $`=`$ $`{\displaystyle \frac{1}{w_3^3}}+{\displaystyle \frac{g(t)/(2\pi )^{3/2}}{w_1w_2w_3^2}}.`$ (31c)
If we consider cylindrically symmetry solutions we have $`w_1=w_2=w`$, this width $`w`$ being equivalent to the width $`X`$ arising in the radially symmetric two-dimensional case.

$$X=\left(|u|^2(x^2+y^2)d^3x\right)^{1/2}=\left(\frac{w_1^2+w_2^2}{2}\right)^{1/2}=w.$$
(32)
Therefore, if we consider the case where the potential is present only along the $`z`$ direction and we suppose cylindrical symmetry the equations for the widths are

$`\ddot{w}`$ $`=`$ $`{\displaystyle \frac{1}{w^3}}+{\displaystyle \frac{g(t)/(2\pi )^{3/2}}{w^3w_3}},`$ (33a)
$`\ddot{w_3}+\lambda _3^2(t)w_3`$ $`=`$ $`{\displaystyle \frac{1}{w_3^3}}+{\displaystyle \frac{g(t)/(2\pi )^{3/2}}{w^2w_3^2}}.`$ (33b)
Our goal is to stabilize some solutions of this model, i.e. to get periodic solutions for $`(w(t),w_3(t))`$. In particular, if we impose that

$$\frac{g(t)}{(2\pi )^{3/2}w_3(t)}=\frac{g_{2D}(t)}{2\pi },$$
(34)
where $`g_{2D}(t)`$ is an specific modulation of the nonlinear term allowing trapping in two dimensions, then Eq. (33a) is the same equation that in the radially symmetric two-dimensional case (Eq. (23)) and, therefore, by taking the three-dimensional modulation of the form of Eq. (34) we will have stabilization of the width $`w`$. Now, Eq. (34) can be written as

$$g(t)=(2\pi )^{1/2}g_{2D}(t)w_3(t),$$
(35)
and the problem arises when we realize that $`w_3(t)`$ is not known, but it evolves according to Eq. (33b), so Eq. (35) is useless in practical cases. Nevertheless, we could think about the possibility of stabilizing $`w_3`$ around the equilibrium point of Eq. (33b) if the nonlinear term were absent. This equilibrium point is $`w_{3,0}=\lambda _3^{1/2}`$. To ensure that the value of $`w_3`$ is as similar to the equilibrium value as possible, we have to impose that

$$\left|\frac{g(t)/(2\pi )^{3/2}}{w(t)^2w_{3,0}^2}\right|\frac{1}{w_{3,0}^3},$$
(36)
for all $`t`$, and this implies that

$$\lambda _3\underset{t}{\mathrm{max}}\left(\left|\frac{g_{2D}(t)}{2\pi w^2(t)}\right|\right).$$
(37)
If Eq. (37) were satisfied we could take $`w_3(t)w_{3,0}=\lambda _3^{1/2}`$ in Eq. (35) and from this equation we get a prediction for the values of the modulation we should take to get stabilization. As we know the values of $`g_{2D}`$ and $`w`$ to stabilize in two dimensions, Eq. (37) is the condition to get stabilization in three dimensions (with the potential along $`z`$). As we said before, this condition means that the more coin-shaped the system is, the better the stabilization will be, leading to a quasi two-dimensional system. Again, the variational method predicts that stabilization can occur only if $`g_0<(2\pi )^{1/2}(2\pi )`$ but we will see later that it is also possible for other $`g_0`$ values.

In Fig. 6 we plot the evolution of the widths after solving numerically Eqs. (33) for different values of $`\lambda _3`$. We have taken $`g_0=(2\pi )^{1/2}(2.2\pi )`$, $`g_1=(2\pi )^{1/2}8\pi `$, $`\mathrm{\Omega }=40`$, $`w(0)=1.09`$ and $`w_3(0)=w_{3,0}=\lambda _3^{1/2}`$. We expect that the values of $`w`$ during the evolution will be about $`w=1`$, as in 2D case, so the condition (37) says that $`\lambda _35`$, approximately. We can see that the greater the value of $`\lambda _3`$ is the better the stabilization is, according to our previous estimates.

We have compared these results with the simulations of the full NLSE for the case when the potential is restricted to the $`z`$-axis and the system is strongly coin-shaped. To do so we have used a pseudo-spectral fully asymmetric 3D evolution method combined with a second order split-step method to advance in time. Typical grid sizes were of $`128\times 128\times 64`$.

In Fig. 7 we present the results for $`g_0=(2\pi )^{1/2}(2\pi )`$, $`g_1=(2\pi )^{1/2}8\pi `$, $`\mathrm{\Omega }=40`$, $`\lambda _3=100`$ and as initial data we take a Gaussian with $`w(0)=1.09`$ and $`w_3(0)=\lambda _3^{1/2}`$. We see that, as in the two-dimensional case, part of the wave packet goes outside and after a readjustment the solution oscillates in the same way that in two dimensions. Moreover, the width $`w_3`$ remains nearly constant during the evolution.

## 6 Conclusions and discussion

In this paper we have provided a deeper understanding of the phenomenon of stabilization of solitons of the nonlinear Schrödinger equation by means of the control of the nonlinear term. We have developed the moment equations for 2D and 3D systems and obtained, on the basis of their rigorous analysis, precise conditions for the stabilization of 2D systems.

Taking as initial data Townes solitons or Gaussians for simulations based on Eq. (3) we have shown that the former is the structure which is stabilized and that other initial data which can be trapped must eject a fraction of the wave packet in the form of radiation to accomodate to this specific solution.

Also we have analyzed the three-dimensional situation. Here we have made an extensive search of stable regions according to our moment equations improving and extending the analysis of Ref. . We have justified why moment-type equations cannot be used to predict the dynamics of the system. Only limited time stabilization is possible when the parameters are fine-tuned to very precise values. Finally, when a strong trapping along one specific direction is kept the system becomes effectively two-dimensional and it can be described again by variational methods whose predictions agree well with the full numerical simulations of the problem. In the latter case three-dimensional confinement is possible although there are only two spatial directions along which the solution is trapped by the nonlinear forces plus the stabilization mechanisms while the other direction is trapped by harmonic forces provided by the potential.

It is remarkable that the nonlinear Schrödinger equations support these stabilized structures which open new fields for applications, in fact these are the first stable structures obtained in the framework of the cubic nonlinear Schrödinger equation.

Many extensions of this work are possible. First, it would be very interesting to study the robustness of these breather-type solutions under different perturbations, e.g. by mutual collision of different structures in single and multicomponent systems. Secondly it could be interesting to try to stabilize other stationary structures different from the two-dimensional Townes soliton.

Acknowledgments

This work has been partially supported by the Ministerio de Ciencia y Tecnología under grants BFM2000-0521, BFM2002-01308 and Consejería de Ciencia y Tecnología de la Junta de Comunidades de Castilla-La Mancha under grant PAC02-002. G. D. Montesinos is supported by Ministerio de Educación, Cultura y Deporte under grant AP2001-0535.


format: text

text:
# Bulk viscosity of strange quark matter in density-dependent quark mass model and dissipation of r-modes in strange stars

## 1 INTRODUCTION

Since Witten conjectured that the strange quark matter(SQM) composed of comparable number of $`u`$, $`d`$ and $`s`$ quarks might be absolutely stable or metastable phase of nuclear matter, many theoretical and observational efforts have been made on the investigation of its properties and potential astrophysical significance. Because the difficulty of quantum chromodynamics (the lattice approach) in the non-perturbation domain , ones used to adopt phenomenological models. One of the famous models is the MIT bag model with which many authors study the equation of state(EOS) of SQM and the mass-radius relation of the assumed strange stars. Another alternative model is density-dependent quark mass model(DDQM). Originally, some pioneer investigations were carried out but the wrong thermodynamics treatment is applied. Latterly, some works dealt with the thermodynamics with density-dependent particle masses self-consistently with which the structural properties of strange stars have been re-discussed.

It had also been recognized the effects of the dynamics of SQM on strange star’s rotations by the time that the EOS of SQM was gradually uncovered. Wang and Lu find that the non-leptonic weak reaction dominates bulk viscosity of SQM. Further investigations of the bulk viscosity also have been carried out by many other authors in MIT bag model and the fruitful results have been obtained. A recent study incorporate the effective quark masses in quasiparticle description into the bulk viscosity and have a significant result for relevant astrophysics. In parallel, a previous work made the study of the bulk viscosity in mass-density-dependent model but the self-consistent thermodynamics are overlooked. Here we will concentrate the problem in the background of the correct thermodynamics and consider the effect of variable quark masses on the bulk viscosity of SQM. As seen latter in this paper, the effect have twofold contributions. The first, the modification of EOS influences the viscosity. The second, the density-dependent masses directly contribute to the dynamical quantity.

We organize this paper as such. In Sec.II, the self-consistent thermodynamics with density-dependent particle masses are recalled and then further the thermodynamical relations are discussed. In Sec.III, the bulk viscosity of SQM with density-dependent particle masses is derived, which arises from the non-leptonic weak interaction. In Sec.IV, the numerical results are given and a application to the recent astrophysical problem is discussed. In Sec. V, a brief summary is made.

## 2 THERMODYNAMICS WITH DENSITY
-DEPENDENT PARTICLE MASSES

The self-consistent thermodynamics in mass-density-dependent model have been discussed in many works. However Peng et al pointed that the thermodynamics treatments given by have serious problem. They re-derived the pressure and energy expressions from the general ensemble theory and found an extra term can add to the expression of pressure but not to the energy formula. Here we don’t repeat the derivation but simply write the results. When particle masses depend on baryon number density or the volume of the system, the pressure and energy have

$$P=\frac{(V\mathrm{\Omega }(V))}{V}|_{T,\{\mu \}},$$
(1)
$$E=\beta \frac{(\beta V\mathrm{\Omega }(V))}{\beta }+\mu N+\beta \frac{(V\mathrm{\Omega })}{\beta }.$$
(2)
Although the thermodynamics potential density depends on the volume, we still obtain the thermodynamic relation

$$E=V\mathrm{\Omega }(V)+\mu N+TS,$$
(3)
where $`\beta `$ is the reverse temperature, $`V`$ is the system volume, $`\mathrm{\Omega }`$ is thermodynamics potential density, $`\mu `$ is the chemical potential and $`N`$ is the particle number. According to Eq.(3), the basic thermodynamic equation is written as at invariable temperature.

$$\mathrm{d}E=\mathrm{d}(V\mathrm{\Omega }(V))+\mathrm{d}(\mu N)=(\mathrm{\Omega }+V\frac{\mathrm{\Omega }}{V})\mathrm{d}V+\mu \mathrm{d}N,$$
(4)
Therefore, the pressure have an extra term

$$P=\mathrm{\Omega }V\frac{\mathrm{\Omega }}{V},$$
(5)
Equally, it is the formula ,$`P=\mathrm{\Omega }+n_b\frac{\mathrm{\Omega }}{n_b}`$, given in . Those previous works have shown that the thermodynamics are self consistent, where an extra term is appended to the pressure but not to the energy. In the framework of self consistent thermodynamics, the EOS of SQM and the structural properties of strange stars have been investigated. At zero temperature, the thermodynamic potential density is

$$\mathrm{\Omega }=\underset{i}{}\frac{1}{8\pi ^2}\left[\mu _i(\mu _i^2m_i^2)^{1/2}(2\mu _i^25m_i^2)+3m_i^4\mathrm{ln}\frac{\mu _i+\sqrt{\mu _i^2m_i^2}}{m_i}\right].$$
(6)
Accordingly, the energy and pressure can be obtained from Eqs (3) and (5), and the number density is

$$n_i=\frac{1}{\pi }(\mu _i^2m_i^2)^{3/2}.$$
(7)
It has been found that the density-dependent quark masses should take the form

$$m_i=m_{i0}+\frac{D}{n_b^{1/3}},i=u,d,s,$$
(8)
where $`D`$ is a parameter to be determined by stability arguments. For given values $`n_b`$,

$$n_b=\frac{n_u+n_d+n_s}{3},$$
(9)
we very easily determine the quark chemical potentials $`\mu _i`$ to evaluate the thermodynamic quantities of the system by the beta equilibrium and the charge neutrality conditions

$$\mu _d=\mu _s\mu ,\mu _s=\mu _u+\mu _e,$$
(10)
$$\frac{2}{3}n_u\frac{1}{3}n_d\frac{1}{3}n_sn_e=0.$$
(11)
In the following, we will consider the dynamics of SQM( bulk viscosity) near the equilibrium state, which is governed by the EOS calculated with the self consistent thermodynamics.

## 3 DERIVATION OF THE BULK VISCOSITY

In a previou study, the formulae in MIT model, which is presented by Wang and Lu and developed by other investigators, was simply introduced to evaluated the bulk viscosity of SQM in density dependent quark mass model. The investigation used an incorrect thermodynamical results without regard to the consequences due to the dependence of the quark masses on baryon number density. In this section, we derive the bulk viscosity of SQM in the self consistent thermodynamics. Since the quark masses are the functions of the baryon number density represented in(6), we for convenience assume the baryon number density instead of the volume as had adopted in previous works oscillates periodically,

$$n_b=n_{b0}+\mathrm{\Delta }n_b\mathrm{sin}(\frac{2\pi t}{\tau })=n_{b0}+\delta n_b,$$
(12)
where $`n_{b0}`$ is the equilibrium baryon number density and $`\mathrm{\Delta }n_b`$ is the amplitude of perturbation, $`\tau `$ is the period. The presence of viscosity results in the energy dissipation of pulsations in bulk matter. On one hand,using the standard definition of the bulk viscosity of matter, energy dissipation rate per unit volume average over the oscillation period $`\tau `$ can be written as

$$\dot{E}_{\mathrm{kin}}=\frac{\zeta }{\tau }_0^\tau dt(\mathrm{div}𝐯)^2.$$
(13)
where $`𝐯`$ is velocity of the fluid. On the other hand, the hydrodynamic matter deviates from the equilibrium accompanied by the time variations of the local pressure, $`P(t)`$. The dissipation of the energy of hydrodynamic motion due to irreversiblity of periodic compression-decompression process can be expressed as

$$\dot{E}_{\mathrm{diss}}=\frac{1}{\tau V}_0^\tau dtP(t)\dot{V}.$$
(14)
According to the continuity equation $`\dot{n}_b+n_{b0}\mathrm{div}𝐯`$ and the specific volume $`V=\frac{1}{n_b}`$, bearing in mind that $`\dot{E}_{\mathrm{kin}}=\dot{E}_{\mathrm{diss}}`$, from Eqs(13)and (14) we have the bulk viscosity

$$\zeta =\frac{1}{\pi }(\frac{n_{b0}}{\mathrm{\Delta }n_b})_0^\tau P(t)\mathrm{cos}\frac{2\pi t}{\tau }.$$
(15)
It is proven that the bulk viscosity of SQM is mainly produced by the nonleptonic interaction as follows

$$u+ds+u$$
(16)
thereby at quasi-equilibrium the pressure can be regarded as a function of the number density $`n_i`$, where $`i`$ takes $`u,d,s`$. From the baryon number density(9), it is suitable to assume that $`P=P(n_b,n_d,n_s)`$. Thus the pressure can expand in small deviations

$$P=P_0+(\frac{P}{n_b})\delta n_b+(\frac{P}{n_d})\delta n_d+(\frac{P}{n_s})\delta n_s,$$
(17)
where all derivatives are taken at equilibrium and the changes in d- and s-quark number density will be composed of two parts

$$\delta n_i=\delta \stackrel{~}{n}_in_i\frac{\delta V}{V},$$
(18)
where $`\delta \stackrel{~}{n}_i`$ denotes the number density at a given volume for s- and d-quark. From the reaction (16), it evidently satisfies the relation below

$$\delta \stackrel{~}{n}_d=\delta \stackrel{~}{n}_s=_0^t\frac{\mathrm{d}\stackrel{~}{n}_d}{\mathrm{d}t}dt.$$
(19)
Due to the baryon number conservation, we also have

$$\frac{\delta n_b}{n_b}=\frac{\delta V}{V}.$$
(20)
Finally, Eq(17) becomes

$$P=P_0+\left[n_b\frac{P}{n_b}+n_d\frac{P}{n_d}+n_s\frac{P}{n_s}\right]\frac{\delta n_b}{n_b}+\left[\frac{P}{n_d}\frac{P}{n_s}\right]_0^t\frac{\mathrm{d}\stackrel{~}{n}_d}{\mathrm{d}t}dt.$$
(21)
Obviously, only the latest term contributes to the integral in Eq(15). In the formula above the net rate per unit volume can be transform as from the rate per mass given by

$$\frac{\mathrm{d}\stackrel{~}{n}_d}{\mathrm{d}t}=\frac{3}{\pi ^3}G_F^2\mathrm{sin}^2\theta \mathrm{cos}^2\theta JT^2\left(1+\frac{\delta \mu ^2}{4\pi ^2T^2}\right)\delta \mu ,$$
(22)
where $`J`$ is given in  and $`\delta \mu =\mu _s\mu _d`$. Analogously to the pressure, the chemical potential difference can be expanded, giving

$$\delta \mu (t)=(\frac{\delta \mu }{n_b})\delta n_b+(\frac{\delta \mu }{n_d})\delta n_d+(\frac{\delta \mu }{n_s})\delta n_s.$$
(23)
To obtain the reaction rate and then calculate the bulk viscosity, the acquire of the derivatives $`\frac{P}{n_i}`$ and $`\frac{\mu _i}{n_i}`$ is necessary with the equation of state in section above. We respectively regard the thermodynamic potential and the chemical potential as the functional form below

$$\mathrm{\Omega }=\mathrm{\Omega }(\mu _j(n_i,n_b(n_i)),n_b(n_i)),\mu _i=\mu _i(n_i,n_b(n_j))$$
(24)
where $`i,j`$ take $`u,d,s`$. From (5) and (7), we immediately obtain

$$\frac{P}{n_j}=\underset{i}{}(n_i3n_b)\frac{\mu _i}{n_j},$$
(25)
$$\frac{\mu _i}{n_i}=\frac{\mu _i^2m_i^2}{3\mu _in_i},\frac{\mu _i}{n_b}=\frac{m_i}{\mu _i}\frac{m_i}{n_b}.$$
(26)
Thus, the part in(21) contributed to the integral (15), $`\delta P`$, and the chemical potential difference at quasi-equilibrium, $`\delta \mu `$ can be expressed as

$$\delta P=\left(\frac{n_d3n_b}{n_d}\frac{\mu ^2m_d^2}{3\mu }\frac{n_s3n_b}{n_s}\frac{\mu ^2m_s^2}{3\mu }\right)_0^t\frac{\mathrm{d}\stackrel{~}{n}_d}{\mathrm{d}t},$$
(27)
$$\delta \mu =\left(\frac{n_b(m_sm_d)}{\mu }\frac{m_d}{n_b}\frac{m_sm_d}{3\mu }\right)\frac{\delta n_b}{n_b}\left(\frac{\mu ^2m_d^2}{3\mu n_d}+\frac{\mu ^2m_s^2}{3\mu n_s}\right)_0^t\frac{\mathrm{d}\stackrel{~}{n}_d}{\mathrm{d}t}.$$
(28)
Substituting Eq(27) into Eq(15), we finally get the viscosity

$$\zeta =\frac{1}{\pi }\frac{n_{b0}}{\mathrm{\Delta }n_b}\left(\frac{n_d3n_b}{n_d}\frac{\mu ^2m_d^2}{3\mu }\frac{n_s3n_b}{n_s}\frac{\mu ^2m_s^2}{3\mu }\right)_0^\tau dt_0^t\frac{\mathrm{d}\stackrel{~}{n}_d}{\mathrm{d}t}.$$
(29)
For the sake of convenience, Eq(28) can be equivalently rewritten as

$$\frac{\delta \mu }{t}=\left(\frac{n_b(m_sm_d)}{\mu }\frac{m_d}{n_b}\frac{m_sm_d}{3\mu }\right)\frac{\delta n_b}{n_b}\left(\frac{\mu ^2m_d^2}{3\mu n_d}+\frac{\mu ^2m_s^2}{3\mu n_s}\right)\frac{\mathrm{d}\stackrel{~}{n}_d}{\mathrm{d}t}.$$
(30)
## 4 RESULTS AND DISCUSSION

Eqs (22) and (30) must be solved numerically to give the rate$`\frac{\mathrm{d}\stackrel{~}{n}_d}{\mathrm{d}t}`$ and the numerical result of Eq(29) can be obtained.

In figure1 we have plotted viscosity vs. relative perturbation $`\frac{\mathrm{\Delta }n_b}{n_b}`$ respectively for (i)$`D^{\frac{1}{2}}=156`$MeV, $`m_{s0}=80`$MeV at $`n_b=0.4`$fm-3 and $`\tau =10^3`$(solid curves) and (ii)$`D^{\frac{1}{2}}=156`$MeV, $`m_{s0}=140`$MeV at $`n_b=1.36`$fm-3 and $`\tau =10^3`$(dashed curves).

In figure2 we gave the results from our formulae(solid curves) and those in ref.(dashed curves) for the same parameters $`n_b=0.4\mathrm{fm}^3,D^{\frac{1}{2}}=156.0\mathrm{Mev},m_{s0}=80\mathrm{M}\mathrm{e}\mathrm{v}`$. Evidently our viscosity can be higher.

In figure3 we show the results in MIT bag model from Zheng et al for $`m_{s0}=80\mathrm{M}\mathrm{e}\mathrm{V},\tau =10^3`$. The solid curves are because of the consideration of screening mass effect and the dashed curves are regardless of the coupling among quarks in MIT bag model. Comparing figure2 with solid curves in figure3, we notice that our results are higher by 2 to 3 order of magnitude than those of bag model. We also find that the solid curves in figure2 is also higher than the dashed curves. In a word, the facts testify to the viscosity increase of SQM due to medium effect.

Recently it has been realized that rotating relativistic stars are generically unstable against the r(otational)-mode instability The viscous dissipation of dense matter has relevance to the onset of r-mode instability, gravitational radiation detection and the evolution of pulsars. A main motivation for the studies of viscosity is to calculate the damping time for compact star vibrations. For a star of constant density, the r-mode energy is estimated as

$$E=\frac{\alpha ^2\pi }{2m(2m+3)}(m+1)^3(2m+1)!\rho R^5\mathrm{\Omega }^2.$$
(31)
The mode dissipation energy can be computed in term of

$$\frac{\mathrm{d}E}{\mathrm{d}t}=\zeta \delta \sigma \delta \sigma ^{}\mathrm{d}^3x.$$
(32)
The damping time scale is thus

$$\tau _D=\frac{2E}{\frac{\mathrm{d}E}{\mathrm{d}t}}=1.31\times 10^{15}\rho ^2\mathrm{\Omega }^4\zeta ^1.$$
(33)
Figure 4-7show the damping time as functions of temperature for some typical spinning strange star with periods $`𝒯=1.5`$ms, 3ms, 15ms and Keplian limit. In comparison with harmonic oscillations, the damping time scale of r-mode with period $`\tau =\frac{2}{3}𝒯=10^3`$s(figure4) is only short to about 1s, but that of harmonic oscillations can be short to $`10^3`$s(see figure 4 in Ref 17). This implies that the simple radial pulsations can be damped more efficiently by the bulk viscosity. For a Keplian rotation star, the minimum time scale is about $`10^1`$s.

## 5 CONCLUSION

Considering the self-consistent thermodynamics in density dependent quark mass model, we correct the error in and re-derive the bulk viscosity of SQM. The viscosity in DDQM is larger than MIT bag model. Consequences is that the r-modes of strange star will be damped more efficiently in DDQM implying the r-mode instability window of strange stars should be modified.


text:
# Dimension Compactification – a Possible Explanation for Superclusters and for Empirical Evidence Usually Interpreted as Dark Matter

## Abstract

According to modern quantum physics, at the microlevel, the dimension of space-time is $`11`$; we only observe 4 dimensions because the others are compactified: the size along each of the other dimensions is much smaller than the macroscale. There is no universally accepted explanation of why exactly 4 dimensions remain at the microscopic level. A natural suggestion is: maybe there is no fundamental reason why exactly 4 dimensions should remain, maybe when we go to even larger scales, some of the 4 dimensions will be compactified as well? In this paper, we explore the consequences of the compactification suggestion, and we show that – on the qualitative level – these consequences have actually been already observed: as superclusters and as evidence that is usually interpreted as justifying dark matter. Thus, we get a new possible explanation of both superclusters and dark matter evidence – via dimension compactification.

Keywords: dimension compactification, superclusters, dark matter

Main idea. According to modern quantum physics, at the microlevel, space-time has at least 11 dimensions; we only observe 4 dimensions in our macroobservations because the rest are compactified: the size along each of the remaining directions is so small that on macrolevel, these dimensions can be safely ignored; see, e.g., \[Green et al. 1988, Polchinski 1998\].

There is no universally accepted explanation of why exactly 4 dimensions remain at the microscopic level. A natural suggestion is: maybe there is no fundamental reason why exactly 4 dimensions should remain. Maybe when we go to even larger scales, some of the 4 dimensions will be compactified as well?

Could one rigorously argue for “compactification” when these effects must occur over the accessible, kiloparsec and larger, scales? In modern physics, indeed, compactification is related to quantum-size distances, but there is nothing inherently quantum in the compactification idea. Indeed, the very first paper that proposed compactification – the 1938 paper by A. Einstein and P. Bergmann \[Einstein and Bergmann 1938\] – described it as simply the extra dimension being a circle, so that the entire space-time looks like a thin cylinder whose width is negligible if we operate at large enough scales.

In this paper, we explore the consequences of the compactification suggestion, and we show that – on the qualitative level – these consequences have actually been already observed: as superclusters and as evidence that is usually interpreted as justifying dark matter. Thus, we get a new possible explanation of both superclusters and dark matter evidence – via dimension compactification.

Geometric consequences of the main idea. If our idea is correct, then, as we increase the scale, we will observe how a 3D picture is replaced by a 2D and then by a 1D one. This is exactly what we observe: while at a macrolevel, we see a uniform distribution of galaxies in a 3D space, at a larger-scale level, galaxies start forming superclusters – long and thin strands of clusters and galaxies; see, e.g., \[Fairall 1998\]. Superclusters are either close to a 2D shape (as the “Great Wall” discovered in the 1980s) or close to 1D.

Towards physical consequences of the main idea. How does the change in dimension affect physics? In non-relativistic approximation, the gravitation potential $`\phi `$ is related to the mass density $`\rho `$ by the Laplace equation $`\mathrm{}\phi =\rho `$. In the 3D space, this leads to Newton’s potential $`1/r`$, with the force

$$F(r)=\frac{G_0m_1m_2}{r^2};$$

in a 2D space, we get potential $`\mathrm{log}(r)`$, with the force

$$F(r)=\frac{G_1m_1m_2}{r}$$

for some constant $`G_1`$. For intermediate scales, it is reasonable to consider a combination of these two forces:

$$F(r)=\frac{G_0m_1m_2}{r^2}+\frac{G_1m_1m_2}{r}.$$
$`(1)`$
In the Appendix, we explain, on the qualitative level, why such combination naturally follows from the above compactification idea.

Let us consider the simplest possible compactification (along the lines of the original Einstein-Bergmann paper), where the space-time is wrapped as a cylinder of circumference $`R`$ along one of the coordinates. What will the gravitational potential look like in this simple model? The easiest way to solve the corresponding Newton’s equation in this cylinder is to “unwrap” the cylinder into a full space. After this “unwrapping”, each particle in a cylindrical space (in particular, each source of gravitation) is represented as infinitely many different bodies at distances $`R`$, $`2R`$, etc., from each other. For further simplicity, let us consider the potential force between the two bodies on the wrapping line at a distance $`r`$ from each other. The effect of the second body on the first one in cylindrical space is equivalent to the joint effect of multiple copies of the second body in the unwrapped space:

The resulting gravitational potential of a unit mass can be described as a sum of potentials corresponding to all these copies, i.e.,

$$\phi (r)=\frac{1}{r}+\frac{1}{r+R}+\frac{1}{r+2R}+\mathrm{}+\frac{1}{r+kR}+\mathrm{}$$
$`(2)`$
From the purely mathematical viewpoint, this sum is infinite. From the physical viewpoint, however, the actual potential is not infinite: due to relativistic effects, at the current moment of time $`t_0`$, the influence of a source at a distance $`d=r+kR`$ is determined by this source’s location at a time $`t_0d/c`$ (where $`c`$ is the speed of light). Thus, we only need to add the terms for which $`d/c`$ is smaller than the age of the Universe. As a result, we can ignore the slowly increasing infiniteness of the sum when $`k\mathrm{}`$.

How can we estimate this potential? The formula (2) has the form

$$f(0)+f(1)+\mathrm{}+f(k)+\mathrm{},$$

where we denoted $`f(x)\stackrel{\mathrm{def}}{=}1/(r+xR)`$. It is difficult to get an analytical expression for the exact sum, but we can use the fact that this sum is an integral sum for an integral $`{\displaystyle _0^{\mathrm{}}}f(x)𝑑x`$; this integral has an analytical expression – it is $`\mathrm{const}(1/R)\mathrm{ln}(r)`$. In this approximation

$$_0^{\mathrm{}}f(x)𝑑xf(0)+f(1)+\mathrm{}+f(k)+\mathrm{},$$
$`(3)`$
we used the fact that

$$_0^{\mathrm{}}f(x)𝑑x=_0^1f(x)𝑑x+_1^2f(x)𝑑x+\mathrm{}+_k^{k+1}f(x)𝑑x+\mathrm{},$$

and we approximated each term $`{\displaystyle _k^{k+1}}f(x)𝑑x`$ by $`f(k)`$. This approximation is equivalent to approximating the function $`f(x)`$ on the interval $`[k,k+1]`$ by its value $`f(k)`$ at the left endpoint of this interval – i.e., by the first term in the Taylor expansion of the function $`f(x)`$ around the point $`k`$. A natural next approximation is when instead of only taking the first term, we consider the first two terms in this Taylor expansion, i.e., when on each interval $`[k,k+1]`$, we approximate the function $`f(x)`$ by a formula $`f(k)+f^{}(k)(xk)`$. Under this approximation,

$$_k^{k+1}f(x)𝑑xf(k)+\frac{1}{2}f^{}(k),$$

and therefore,

$$_0^{\mathrm{}}f(x)𝑑x=_0^1f(x)𝑑x+_1^2f(x)𝑑x+\mathrm{}+_k^{k+1}f(x)𝑑x+\mathrm{}=$$

$$(f(0)+f(1)+\mathrm{}+f(k)+\mathrm{})+\frac{1}{2}(f^{}(0)+f^{}(1)+\mathrm{}+f^{}(k)+\mathrm{}).$$
$`(4)`$
The second term in the right-hand side can be (similarly to the formula (3)) approximated by

$$f^{}(0)+f^{}(1)+\mathrm{}+f^{}(k)+\mathrm{}\mathrm{const}+_0^{\mathrm{}}f^{}(x)𝑑x\frac{1}{2}f(0).$$

Thus, from the formula (4), we can conclude that

$$f(0)+f(1)+\mathrm{}+f(k)+\mathrm{}_0^{\mathrm{}}f(x)𝑑x+f(0).$$

In particular, for our function $`f(x)`$, we get

$$\phi (r)=f(0)+f(1)+\mathrm{}+f(k)+\mathrm{}\mathrm{const}\frac{1}{R}\mathrm{ln}(r)+\frac{1}{2}\frac{1}{r}.$$

Differentiating relative to $`r`$, we get the desired formula (1) for the gravitational force:

$$F(r)=\frac{d\phi (r)}{dr}=\frac{1}{Rr}+\frac{1}{2r^2},$$

with $`R=2G_0/G_1`$. According to estimates from \[Kuhn and Krugliak 1987\], we expect $`R`$ to be between $`10`$ and $`30`$ kpc.

Physical consequences of the main idea. The force described by formula (1) is exactly the force that, according to Kuhn, Milgrom, et al. \[Kuhn and Krugliak 1987, Milgrom 1983, Milgrom 2002, Sanders and McGaugh 2002\], is empirically needed to describe the observations if we want to avoid dark matter. Indeed, in Newtonian mechanics, for any large-scale rotating gravitational system, if we know the rotation speed $`v`$ at a distance $`r`$ from the center, we can find the mass $`M_g(r)`$ inside the sphere of radius $`r`$ by equating the acceleration $`v^2/r`$ with the acceleration $`G_0M_g(r)/r^2`$ provided by the Newton’s law. As a result, we get $`M_g(r)=rv^2/G_0`$. Alternatively, we can also count masses of different observed bodies and get $`M_L(r)`$ – the total mass of luminescent bodies. It turns out that $`M_L(r)M_g(r)`$ – hence the traditional explanation that in addition to luminescent bodies, there is also “dark” (non-luminescent) matter.

An alternative explanation is not to introduce any new unknown type of matter – i.e., assume that $`M_g(r)M_L(r)`$ – but rather change the expression for the force, or, equivalently, assume that the gravitational constant $`G_0`$ is not a constant but it may depend on $`r`$: $`G_0=G(r)`$. Equating the acceleration $`v^2/r`$ with the acceleration $`G(r)M_L(r)/r^2`$ provided by the new gravity law, we can determine $`G(r)`$ as $`G(r)=v^2r/M_L(r)`$. Observation data show that $`G(r)=G_0+G_1r`$ for some constant $`G_1`$ – i.e., that the dependence of the gravity force on distance is described by the formula

$$F(r)=\frac{G(r)m_1m_2}{r^2}=\frac{G_0m_1m_2}{r^2}+\frac{G_1m_1m_2}{r},$$

which is exactly what we deduced from our dimension compactification idea.

This idea has been proposed 20 years ago, and one of the reasons why it has not been universally accepted is that it was difficult to get a natural field theory explanation of this empirical law. We have just shown that such an empirical explanation comes naturally if we consider the possibility of dimension compactification.

This explanation is in line with Milgrom’s own explanation; see, e.g., \[Milgrom 2003\].

As we go further up in scale, one more dimension starts compactifying, so we start getting a 1D space in which Laplace equation leads to potential $`\phi (r)r`$. Thus, at a large-scale level, we should have a term proportional to $`r`$ added to the normal gravity potential formulas. This additional term is exactly what is add when we take a cosmological constant $`\mathrm{\Lambda }`$ into consideration.

Observable predictions of our new idea. A possible observable consequence of the additional term $`\phi (r)r`$ is that it leads to an additional constant term in the gravitational force and therefore, to a formula $`G(r)=G_0+G_1r+G_2r^2`$. Thus, if the empirical dependence of $`G(r)`$ on $`r`$ turns out to be not exactly linear but rather slightly quadratic, it will be a strong argument in favor of our compactification idea.

Natural open questions. In this paper, we simply formulate the idea and explain why we believe this idea to be prospective. Many related questions are still open:

* how to concoct a geometry that accomplishes “compactification” – what would a metric look like that transitions from large to small scales?
* could we solve for the force law in such a geometry with some rigor?
* will such a profound geometrical perturbation as the one we propose have other consequences beyond a force law change? perhaps not, but it is worth investigating.

Acknowledgments. This work was supported in part by NASA under cooperative agreement NCC5-209, by the Future Aerospace Science and Technology Program (FAST) Center for Structural Integrity of Aerospace Systems, effort sponsored by the Air Force Office of Scientific Research, Air Force Materiel Command, USAF, under grant F49620-00-1-0365, and by NSF grants EAR-0112968 and EAR-0225670.

The author is very thankful to Jeffrey R. Kuhn for his ideas and for his valuable advise.