text:
# Dynamics of the Bose-Einstein condensate: quasi-one-dimension and beyond

## 1 Introduction

A family of Bose-Einstein condensates (BECs) has been created in a series of ground-breaking experiments  over the last five years. These condensates, made of trapped, dilute gases of 23Na, 87Rb, 1H, or 7Li, are well-modeled by the nonlinear Schrödinger equation (NLSE) , also known as the Gross-Pitaevskii equation . However, they have yet to be experimentally investigated in quasi-one-dimension (quasi-1D.) We define the BEC to be in the quasi-1D regime when its transverse dimensions are on the order of its healing length and its longitudinal dimension is much longer than its transverse ones. In this case the 1D limit of the 3D NLSE is appropriate, rather than true 1D or 2D mean field theories .

The 1D NLSE is expected to be the correct description of the quasi-1D BEC . It is a ubiquitous equation. Among other natural phenomena, it models Bose-condensed photons , disordered media , helical excitations of a vortex line , and light pulses in optical fibers . In this latter context soliton solutions to the NLSE have been used to greatly increase the efficiency of communications .

Because the BEC is a coherent matter-wave it has a natural analogy with light. This insight has led to light-inspired experiments, as for example four-wave mixing  and both continuous-wave and pulsed atom lasers . The creation of a quasi-1D BEC would allow the exploration of a great body of experimental and theoretical results from fiber optics .

In addition to the analogy with fiber optics there are several other reasons to explore the quasi-1D BEC. Cold atom waveguides have recently been experimentally realized . This work suggests the manipulation of matter waves in quasi-1D and perhaps the eventual creation of matter-wave circuits, just as fiber optics have come to replace wires in certain applications. The quasi-1D regime of the BEC is qualitatively different from its 2D and 3D counterparts: in the attractive case collapse is prevented, while in both the repulsive and attractive cases solitons are stablized.

In a pair of recent articles all stationary states of the 1D NLSE under box and periodic boundary conditions  were analytically elucidated. This provides a starting point for the analysis of the quasi-1D BEC. In the present investigation we impose box or periodic boundary conditions in the longitudinal dimension and box boundary conditions in the transverse dimensions. These boundary conditions model the potential quasi-1D regime of a number of present experiments: the atom waveguide ; the prolate harmonic trap in which is formed a cigar-shaped BEC ; the toroidal trap ; and finally an oblate harmonic trap with a barrier formed in the middle either by a second spin state of the same atom  or by a laser , in which is formed a pancake-shaped BEC with the center removed. Periodic boundary conditions are the best first model for toroidal geometries; box boundary conditions are a good starting model for cigar-shaped geometries.

Our approach is to find the stationary states of the NLSE under simple, analytically solvable boundary conditions. These give physical insight into the corresponding results for potentials which require numerical solution. Here we demonstrate the merits of this approach by the extension of 1D into quasi-1D, 2D, and 3D. Elsewhere we have applied such insight to the harmonic potential , among others.

The paper is organized as follows. In section 2 we review the analytic stationary solutions in 1D. In section 3 we numerically illustrate these solutions and their properties and predict how they should extend into quasi-1D and higher dimensions. In section 4 we numerically verify these predictions and thereby demonstrate that the quasi-1D regime is experimentally accessible for both repulsive and attractive condensates. In section 5 we show how quasi-1D-type stationary states extend into 2D and 3D and the rich phenomena that result. Finally in section 6 we give our conclusions.

## 2 Solitons and stationary states

We consider two cases of the BEC: repulsive and attractive atomic pair interactions. For our purposes they differ by the sign of the nonlinear term in the NLSE. The former is the one that has received the most experimental attention. The latter is generally unstable in 3D  but stable in 1D .

The NLSE in 1D has a number of special properties which are described in the mathematical literature. It is integrable , may be solved exactly by the Inverse Scattering Transform  and has a countably infinite number of conserved quantities . We have taken advantage of its special properties in finding its stationary solutions. We first describe all stationary states of the 1D NLSE subject to box and periodic boundary conditions. We then make clear their connection with solitons and provide physical interpretations of their form.

### 2.1 Mathematical form

We begin with the NLSE that describes a BEC of $`N`$ atoms of mass $`M`$, confined in an external potential $`V(\stackrel{}{r})`$:

$$[\frac{\mathrm{}^2}{2M}^2+g\psi (\stackrel{}{r},t)^2+V(\stackrel{}{r})]\psi (\stackrel{}{r},t)=ı\mathrm{}_t\psi (\stackrel{}{r},t)$$
(1)
where $`\psi (\stackrel{}{r},t)^2`$ is the single particle density such that $`\rho (\stackrel{}{r},t)=N\psi (\stackrel{}{r},t)^2`$, the coupling constant $`g4\pi \mathrm{}^2aN/M`$, and $`a`$ is the scattering length for binary collisions between atoms. The case of repulsive interactions corresponds to $`a>0`$; that of attractive interactions to $`a<0`$.

The characteristic length scale of variations in the condensate wavefunction is the healing length $`\xi `$:

$$\xi (8\pi \overline{\rho }a)^{1/2}$$
(2)
where $`\overline{\rho }`$ is the mean particle density. The BEC is in the quasi-1D regime when the lengths of the tranverse $`y`$ and $`z`$ dimensions of the trap, $`L_y`$ and $`L_z`$, satisfy the following criteria: $`L_y,L_z\xi _{eff}`$ and $`L_y,L_z<<L_x`$. The former ensures that the condensate remains in the ground state in the two transverse dimensions while the latter ensures that longitudinal excitations are much lower in energy than possible transverse excitations. Under these conditions one may make an adiabatic separation of longitudinal and transverse variables . Equation (1) then reduces to:

$$[(\xi _{eff})^2_x^2\pm f(x)^2+\stackrel{~}{V}(x)]f(x)=ı_{\stackrel{~}{t}}f(x)$$
(3)
where all terms have been made dimensionless: $`f(x)`$ is the longitudinal portion of the wavefunction; $`\stackrel{~}{t}(\mathrm{}/2M\xi ^2)t`$ we will call the natural time; the $`\pm `$ refer to the repulsive and attractive cases, respectively; $`\stackrel{~}{V}(x)(2M\xi ^2/\mathrm{}^2)V(x)`$; and $`\xi _{eff}`$ is an effective healing length, discussed in section 4.1. If we furthermore assume a stationary state:

$$[(\xi _{eff})^2_x^2\pm f(x)^2+\stackrel{~}{V}(x)]f(x)=\stackrel{~}{\mu }f(x)$$
(4)
where $`\stackrel{~}{\mu }(2M\xi ^2/\mathrm{}^2)\mu `$ is a dimensionless chemical potential which is now an eigenvalue.

All stationary solutions to equation (4) may be written in terms of Jacobian elliptic functions . The properties of such functions are reviewed elsewhere . There are five normalizable symmetry-breaking solution-types to the stationary NLSE with a constant potential on a finite interval. They are pictured in figure 1. The stationary solutions appropriate to longitudinal box boundary conditions are:

$$f(x)=A\mathrm{sn}(\frac{2jK(m)x}{L_x}+\delta m)$$
(5)
$$f(x)=A\mathrm{cn}(\frac{2jK(m)x}{L_x}+\delta m)$$
(6)
for the repulsive and attractive cases, respectively. Equation (5) is shown in figure 1(a) and equation (6) is shown as the anti-symmetric solution in figure 1(c). $`A`$ is the amplitude, $`j1`$ is the number of nodes, $`K(m)`$ is the complete Jacobian elliptic integral which is the quarter period of these functions, $`\delta `$ is the offset, here equal to zero, and $`0m1`$ is the Jacobian elliptic parameter. The general notation sn$`(um)`$ is standard for Jacobian elliptic functions .

Under periodic boundary conditions equations (5) and (6) are still solutions, but $`j`$ is required to be even and the number of nodes is $`j`$, i.e. $`j=2,4,6,\mathrm{}`$ Translational symmetry is restored, so the offset $`\delta `$ is no longer required to be zero. This leads to a Kosterlitz-Thouless-type entropy in 1D .

In addition to the above two solution-types there are three nodeless, symmetry-breaking solution-types under periodic boundary conditions. The first, a solution to the attractive case, is real:

$$f(x)=A\mathrm{dn}(\frac{2jK(m)x}{L_x}+\delta m)$$
(7)
An example is shown as the symmetric solution in figure 1(c). The other two types are intrinsically complex. $`f(x)r(x)\mathrm{exp}(ı\varphi (x))`$ and for the repulsive and attractive cases, respectively:

$$r(x)^2=A^2(1\gamma \mathrm{dn}^2(\frac{2jK(m)x}{L_x}+\delta m))$$
(8)
$$r(x)^2=A^2(\mathrm{dn}^2(\frac{2jK(m)x}{L_x}m)\gamma (1m))$$
(9)
where in both cases the phase must be found by numerical integration from the equation:

$$\varphi ^{}(x)=\frac{\alpha }{r(x)^2}$$
(10)
The phase and amplitude are shown in figure 1(b) and 1(d). In the repulsive case $`A^2\gamma `$ is the depth of the density minima below the constant background. When $`\gamma =1`$ equation (5) is recovered. In the attractive case $`\gamma `$ interpolates between the real anti-symmetric and symmetric solutions in equations (6) and (7). In both cases $`0\gamma 1`$. $`\alpha `$ is a constant of integration. $`j`$ is the number of density minima or maxima, respectively.

All five solution-types are found by solving equation (4) subject to normalization and boundary conditions, and are described in detail in references  and .

### 2.2 Physical interpretation

The nonlinear mean field term in equation (4) causes a spreading out or a clumping together of the wavefunction between nodes for the repulsive and attractive cases, respectively. Thus the stationary solutions to the NLSE with box boundary conditions are in one-to-one correspondence with the sinusoidal solutions to the particle-in-a-box problem in linear quantum mechanics. The number of nodes is $`j1`$ and the complete set of stationary states is described by equations (5) and (6), where $`j\{1,2,3,\mathrm{}\}`$.

Likewise, if $`j`$ is restricted such that $`j\{2,4,6,\mathrm{}\}`$, equations (5) and (6) are in one-to-one correspondence with the sinusoidal solutions to the particle-on-a-ring problem in linear quantum mechanics. As the ring has rotational symmetry these solutions are symmetry-breaking. We note that, as in the linear Schrödinger equation, there are also complex, constant-amplitude, plane-wave solutions on the ring.

The Jacobian elliptic parameter $`m`$ governs the strength of the nonlinearity. As $`m0^+`$, sn $``$ sin and cn $``$ cos, respectively. This is the linear, sinusoidal limit. As $`m1^{}`$, sn $``$ tanh and cn $``$ sech. These are the dark and bright soliton, stationary solutions to the NLSE on the infinite line, respectively . This shows the connection between these stationary solutions and solitons. In the box the $`j^{\mathrm{th}}`$ stationary state is a $`j1`$ or $`j`$ soliton-train in the repulsive and attractive cases, while on the ring the $`2j^{\mathrm{th}}`$ stationary state is a $`2j`$ soliton-train. As we have shown elsewhere , soliton-trains interact in a manner similiar to single solitons. And as with single solitons their form and phase relationships are quite robust.

The other three solution-types have no analogue in the particle-on-a-ring problem in linear quantum mechanics. They exist only on the ring and are nodeless and symmetry-breaking. The complex solutions have a monotonically increasing phase. Each has a complex-conjugate, degenerate partner. The complex solutions for the attractive case, in particular, are an entirely new stationary solution-type for the NLSE.

The complex solutions for the repulsive case are interpreted as density-notch solitons moving with speed $`c`$ on the ring with an opposing momentum boost of the condensate of speed $`c`$, which results in a stationary state in the lab frame. Density-notch solitons have a speed between zero and the Boguliubov sound speed, ranging from maximal to zero depth, respectively . Those not of maximal depth are called grey solitons, while those which are of maximal depth and therefore form a node are called dark. Figure 1(b) shows the bounded, quantized version of a grey 2-soliton-train.

At typical experimental trap sizes the single density-notch stationary state on the ring is the lowest energy, symmetry-breaking excitation above the real, constant-amplitude ground state. When any of these stationary excited states is perturbed it gives rise to soliton-type motions . Recent reports suggest that such motions can be induced in repulsive BECs by optical phase engineering of the condensate phase .

All attractive symmetry-breaking, longitudinally periodic, stationary solution-types, i.e. the anti-symmetric and symmetric ones shown in figure 1(b) and the complex one shown in figure 1(d), are described by the $`C_j`$ point symmetry group, where $`j`$ is the number of peaks. There are $`j`$ nearly degenerate solutions. For even $`j`$ there is a real symmetric-anti-symmetric pair and $`(j2)/2`$ degenerate complex pairs. For odd $`j`$ there is a real symmetric solution and $`(j1)/2`$ degenerate complex pairs. Thus by group theory these three solution-types form the complete set of stationary states made of evenly spaced peaks.

As the attractive BEC has been made in 3D but not in quasi-1D, these stationary solutions are open to experimental investigation.

### 2.3 Phase engineering of dynamics

There are two ways to impart velocity to a density-notch soliton. One may boost the constant background via a *phase ramp* and the notch will drift with the resulting super-current. One may also apply a *phase jump* and the resulting soliton speed with respect to the background will be $`c=c_{max}\mathrm{sin}(\delta /2)`$, where $`c_{max}`$ is the Boguliubov sound speed and $`\delta `$ is related to the phase difference across the moving density notch . Both types of phase profile are apparent in figure 1(b). The phase ramp may be seen in the constant background slope in the phase. The phase jump is shown over each of the two notches. The velocities cancel and a stationary state results.

To impart a collective velocity to trains of density-notches one applies an equal phase difference across each of the member notches or a phase ramp across the whole train. One may also treat the members of the train individually by applying varying jumps and thereby study solitons interactions. A bright soliton may be set in motion by a phase ramp but not a phase jump. For a phase ramp $`\mathrm{exp}(ıkx)`$ the resulting velocity $`ck`$. There is no upper bound on the velocity of bright solitons. As with density-notches, a train of bright solitons may be manipulated collectively by applying a phase ramp across the whole set or individually by applying different phase ramps across the members of the train.

Soliton-trains are thus collective excitations which can be directly manipulated by simple phase profiles. One may manipulate the members of the train individually or collectively. This makes clear the connection between the above stationary states and solitons.

## 3 One dimension

### 3.1 Stability

We first demonstrate that the five solutions types described in section 2 are indeed stationary. In figure 2 we propagate them in time numerically. In the numerical work presented for one, two, and three dimensional solutions of the NLSE the condensate wave function was expanded as $`\psi (\stackrel{}{r},t)=\mathrm{\Sigma }_ic_i(t)\varphi _i(\stackrel{}{r})`$. The $`\varphi _i(\stackrel{}{r})`$ were chosen as a co-ordinate-discretized, pseudo-spectral basis, and the $`c_i(t)`$ were propagated in time using a fourth order variable-step Runge-Kutta algorithm, made highly efficient by implementation with spatial fast fourier transform techniques. We shall use the same numerical methods throughout our article.

The denumerably infinite number of conserved quantities in the 1D NLSE makes solitons robust. Single solitons, in whatever equations they arise, are in general stable. Here in both the repulsive and attractive cases soliton-trains are stable . In the attractive case the two nodeless solution-types are *cyclically stable*. That is, in response to stochastic perturbation they evolve in a complicated manner but come asymptotically close to their original state in an aperiodic cycle. In figure 3 we illustrate the stability properties of the stationary states by introducing $`1\%`$ stochastic white noise into all five fundamental stationary solution types and observing their evolution.

Stochastic perturbation of density-notch solitons can cause a velocity shift or small-scale emission of phonons . The phase relationships between solitons, however, remain intact. The depth of a density-notch is a simple function of its velocity. Therefore the soliton-train changes in depth but not in spacing or form. It drifts as a unit in a diffusive manner. This is true for both the real and complex solution-types.

Stochastic perturbation of bright solitons can cause a velocity shift, small-scale emission of phonons, and/or a phase shift . Although a single soliton is stable, multiple solitons are cyclically stable. As the phase difference between adjacent solitons changes from zero to $`\pi `$ their interaction changes from attractive to repulsive , and it is this changing interaction that leads to cyclical stability. The anti-symmetric solution is stable. The repulsive interactions between solitons lock them into place. In the other two solution-types the solitons attract and repel as their relative phases are perturbed by the phonons that make up the white noise. We have not shown their recurrence here.

### 3.2 Expectations in higher dimensionality

We expect that in quasi-1D the five solution-types described analytically in section 2 and numerically in section 3.1 will have the same properties as in 1D. We also expect, and will verify in the following section, that excellent quasi-1D stationary states can be found as the direct product of eigenstates of the 1D NLSE in each of the dimensions, even though the NLSE is only approximately separable.

To reiterate our criteria for quasi-1D from section 2.1: $`L_y,L_z\xi _{eff}`$; $`L_y,L_z<<L_x`$. In the repulsive case a non-overlapping 2-soliton solution becomes possible when the length is $`2\pi \xi _{eff}`$ . This provides a criterion for the breakdown of the quasi-1D approximation: if $`L_y`$, $`L_z>2\pi \xi _{eff}`$ a transverse excitation becomes likely.

Solutions to the attractive NLSE become unstable when the negative mean field energy is not sufficiently opposed by the positive kinetic energy. This may be seen by performing a variational calculation. So one criterion we could use is the length for which the transverse solution changes from positive to negative energy. This is about $`2\frac{1}{2}\xi _{eff}`$. In the literature it has been found that the 3D attractive NLSE is more likely to go unstable than the 2D NLSE . So we expect that the quasi-1D approximation will break down at a smaller transverse length in 3D than in 2D. In the attractive case we expect that out of the quasi-1D regime the wavefunction should collapse.

We note that in the repulsive case the criteria for separability and quasi-1D are not identical. In the limit as the transverse lengths become very large the NLSE becomes exactly separable for a constant transverse solution. However, the separable solutions can be unstable, as we shall show in section 5.

In the repulsive case it has been shown in the mathematical literature that band soliton and planar soliton solutions of the NLSE in 2D and 3D can decay into vortices via long-wavelength transverse modulations known as snake instabilites . They are called band or planar solitons both to describe their shape in 2D and 3D and to differentiate them from 1D solitons of the type formally obtained by the Inverse Scattering Transform . Thus as we increase the length of the transverse dimensions out of the quasi-1D regime and into 2D and 3D we expect to find that our soliton-train solutions break up into vortices. The BEC conserves vorticity. When the band or planar solitons decay the circulation quantum numbers of the resulting vortices should add to zero. Since the energy of a vortex is proportional to its quantum number squared we expect the decay to be into pairs of singly-quantized vortices of opposite circulation.

## 4 Quasi-one-dimension

To demonstrate the experimental viability of the 1D stationary states described above we add one or two transverse dimensions of scale size $`\xi _{eff}`$ and study the stability and phase engineering of solitons and stationary states in quasi-1D. The connection between our numerical studies and present BEC experiments is made clear. We also illustrate the phase rigidity of the quasi-1D BEC. We hope the following work will encourage experimentalists to create a quasi-1D attractive BEC and to study density-notches and phase rigidity in the quasi-1D regime of the repulsive BEC.

### 4.1 Connection with experiments

In table 1 we show the length and time scales and the number of atoms for the four atomic species of BEC. The below conversion formula allows one to find the number of atoms given the healing length, trap size, and type of BEC.

$$N=\frac{\xi _{eff}L_xL_yL_z}{8\pi a}$$
(11)
where the healing length $`\xi _{eff}`$ and the scattering length $`a`$ must be in the same units and the trap dimensions $`L_x`$, $`L_y`$, and $`L_z`$ are in units of $`\xi _{eff}`$. From section 2.1 the time scale is:

$$t=\frac{2M\xi _{eff}^2}{\mathrm{}}\stackrel{~}{t}$$
(12)
where $`M`$ is the mass of the atomic species. $`\stackrel{~}{t}`$ are the natural, unitless time units we use in the numerical studies.

In the table we choose $`\xi _{eff}=3.0\mu `$m to make solitons in the condensate easily observable. This is the order of the healing length used in recent soliton experiments . The time scales are on the order of 0.1 to 1 ms per natural time unit. Since we present results on time scales from 25 to 100 natural time units all our numerics take place on the experimental time scales of from 1 to 100 ms. Similiar results persist to 1000 natural time units. As BECs last from 1 to 75 s our studies extend over the lifetime of current BECs .

For computational reasons we have used somewhat smaller traps than current experiments, 10x1x1 or 10x$`\frac{1}{2}`$x$`\frac{1}{2}`$ in the attractive case and 25x1x1 in the repulsive case. Thus the number of atoms is small. In the attractive case the trap can be made longer and remain in the quasi-1D regime as long as the number of peaks in the solutions is increased. In the repulsive case the trap size could be increased in each transverse dimension by a factor of 5 and arbitrarily in the longitudinal dimension, resulting in a factor of 100 or more increase in the volume. Therefore the numbers we have used in the simulations scale into the experimentally accessible regime.

The assumption of separability in quasi-1D is made by projecting the total wavefunction onto the ground states of the transverse dimensions and integrating. This produces an effective healing length which depends on the transverse width:

$$\xi _{eff}^2=(\frac{3(E(m_t)/K(m_t)1)^2}{2(1+m_t)E(m_t)/K(m_t)(2+m_t)})^2\xi ^2$$
(13)
$$\xi _{eff}^2=(\frac{3(E(m_t)/K(m_t)(1m_t))^2}{2(1+2m_t)E(m_t)/K(m_t)(25m_t+3m_t^2)})^2\xi ^2$$
(14)
for the repulsive and attractive cases, respectively. $`m_t`$ is the Jacobian elliptic parameter which determines the width of the transverse wavefunction, which is a sn or cn function described by equations (5) or (6). $`E(m_t)`$ and $`K(m_t)`$ are complete Jacobian elliptic integrals.

In the limit as $`m_t0^+`$, $`\xi _{eff}^2\frac{4}{9}\xi ^2`$ in both cases and the transverse wavefunction is sinusoidal. As $`m_t1^{}`$, $`\xi _{eff}^2\xi ^2`$ and $`\xi _{eff}^20`$ monotonically in the repulsive and attractive cases, respectively. For the purposes of the present quasi-1D numerical investigations we shall use a transverse healing length of $`1`$ or $`\frac{1}{2}`$. For these healing lengths $`m_t0.1`$ and $`\xi _{eff}^2\frac{4}{9}\xi ^2`$. We have used these rescalings in determining the number of atoms in table 1.

Besides the ratio of the box lengths to the healing length, the numerical studies do not depend on the choice of the experimental parameters. In all figures we have set $`\xi _{eff}=1`$ and thus are working in length units scaled to the effective healing length.

### 4.2 Stability

As in section 3.1 we again add $`1\%`$ stochastic white noise to the five solution-types, but in both the 2D and 3D quasi-1D limits. The initial state is created from the product of box transverse wavefunctions of the form described by equations (5) and (6) with a box or periodic longitudinal wavefunction of each of the five types shown in figure 1.

In figure 4 we show results of simulations in the 3D case. Trap parameters which enforce the quasi-1D dynamics are used: a volume of 25x1x1 in the repulsive case and 10x$`\frac{1}{2}`$x$`\frac{1}{2}`$ in the attractive case. The stability properties are identical with those found in 1D. Although we have not shown it here, the 2D case displays essentially the same results.

The repulsive solutions drift but are otherwise unchanged. The complex solution-type moves with a steady velocity to the left because the approximation of separability has failed to give the right momentum boost of the condensate to counteract the motion of the soliton-train. But clearly a larger momentum boost would return it to a stationary state. Note that the phase relationships are quite stable: constant between notches in figure 4(a) and changing linearly between notches in figure 4(c).

The attractive anti-symmetric solution is locked into place and is stable as it was in 1D. The two nodeless solution-types break up and recur. The symmetric, real type shown in figure 4(g)-(h) shows an example of this recurrence. In the third time slice the left-hand peak has spread out and changed its phase relationship to the right-hand peak. In the fourth time slice it has returned to its initial state.

As predicted, there are small fluctuations due to the assumption of separability in our choice of initial state. In the repulsive case we found that beyond six healing lengths the solutions became unstable in both 2D and 3D. In the attractive case we found that the requirements for the quasi-1D regime were more complicated. A transverse length of $`\frac{1}{2}\xi _{eff}`$ was required in 3D while $`\xi _{eff}`$ was sufficient in 2D. In general the stability depended both on the number of peaks and the solution-type. The anti-symmetric solution-type with a large number of peaks was the most stable configuration. This kept the density spread out and the contribution of the mean field term in the NLSE small. Other researchers have found similiar results .

### 4.3 Manipulation via phase

In section 2.3 we discussed simple phase profiles which impart a velocity to soliton-trains. Setting the whole train in motion around a thin torus creates a stable, structured super-current in both repulsive and attractive cases. Measuring the velocity of the train is a way to verify the superfluidity of the gaseous BEC. One method would be combination of the new toroidal trap, created to Bose condense Cs , with optical phase-imprinting . Note that the superfluid properties of the attractive BEC are presently unknown.

One may study soliton interactions by treating members of the train separately. In figure 5 we show that planar solitons in quasi-1D interact in a manner similiar to solitons in 1D. In the repulsive case we have used an 0.3$`\pi `$ phase jump to send two density notches towards one another. As may be seen in the figure, they undergo elastic repulsive interaction twice over 50 natural time units.

In the attractive case we have used a real symmetric initial state with a phase ramp of 2$`\pi `$ on each bright soliton. Their first and third interactions show a characteristic smearing of the density which indicates a phase difference between zero and $`\pi `$ . In the second interaction they add coherently as they pass through each other, indicating a phase difference of zero. In their fourth interaction they repel. Therefore it is quite apparent that their relative phase has drifted by approximately $`\pi `$, due to the $`1\%`$ stochastic noise added to the initial state.

For computational reasons we chose phase profiles which caused interaction on a short time scale. But similar interactions could be observed on longer time scales by the appropriate choice of phase jump or ramp. Therefore solitons and stationary states in quasi-1D are connected in the same way as in 1D. Soliton interactions can be studied experimentally in the BEC.

By use of phase ramps one may also study how the BEC responds to shock. The repulsive case has the property of *phase rigidity*, related to symmetry-breaking of the condensate phase. Rigidity is common to many symmetry-breaking materials in condensed matter physics, from liquid crystal to chalk . The BEC exhibits this property by shattering into well-defined phase domains, the size of which are on the order of the healing length. In quasi-1D the borders between domains are made of density-notch solitons.

To illustrate this property we use longitudinal box rather than periodic boundary conditions, with volume 25x1x1. A large phase ramp of $`\varphi (x)=1.6\pi x`$ is applied to a two-density-notch stationary state. In figure 6 we show the time evolution of such a configuration. The first time slice shows the initial state. The second shows a shock wave, apparent in the density. The final time slice shows a highly excited state with well-defined phase domains separated by density notches. Because the kinetic energy of this state is very high the nonlinear term in the NLSE does not contribute signifigantly. Thus a linear beat phenomenon is seen in the density. We have used short time scales for convenience, but a smaller phase ramp would cause the same phenomenon to occur on a much longer, experimentally observable time scale.

## 5 Beyond quasi-one-dimension

In this section the transverse dimensions are extended beyond the quasi-1D regime in both two and three dimensions. Band and planar solitons become unstable, and shocking the condensate produces chaos.

### 5.1 Vortex creation

In figures 7 and 8 we show the decay of band and planar solitons into vortices. In both cases we begin with 2-soliton-train stationary states made in the same way as described in section 4. We use dimensions 25x12 in 2D and 25x12x12 in 3D with periodic boundary conditions in the longitudinal direction. White noise causes the solitons to decay into a pair of vortices of opposite circulation, conserving vorticity. The vortices are singly quantized, as may be seen in the figures by traversing a node in the phase. The colors pass once around the phase color circle, which is defined to be 2$`\pi `$ in extent. In the plots we chose times that would highlight the decay. However, the vortices are stable, spatially stationary, and do not interact over the full length of the study.

White noise in the BEC is principally due to uncondensed atoms and the number of such atoms is a function of temperature. Thus noise is nominally an experimental parameter. The time from creation of the initial stationary solitons to their decay into vortices is a function of the noise. With the correct choice of temperature and volume, stationary band or planar solitons, made via phase engineering , can be used to create vortices in a controlled manner.

### 5.2 Vortex interaction

Band and planar solitons respond to simple phase profiles in the same way as solitons in quasi-1D. The following method permits direct study of vortex interactions.

As explained above, the time from creation of the initial state to decay into vortices is a function of the temperature of the BEC. As in quasi-1D, they move toward each other with a velocity which depends on the phase jump. Thus two solitons can be made to decay into vortices shortly before collision. An example of this process in 2D is shown in figure 9. We use a phase jump of $`\pi `$/15, a noise level of $`0.1\%`$, and a container identical to the one in figure 7. In the first time slice a snake instability and a vortex pair have developed. By $`t=36`$ the vortices are strongly interacting. They have become a pair of vortex dipoles. By $`t=48`$ they have completed their interaction and are moving away from each other around the torus. The time slices in the figure are evenly spaced. Thus it can be seen that, unlike in quasi-1D, the collision is not elastic. The relative vortex velocity in the first three time slices, before the collision, is not the same as in the last two. In fact the two vortex dipoles no longer have the same speed.

Although we have not pictured it here, the evolution of the above collision is very intriguing. Despite the complicated density profile, it is apparent from the phase that the remnants of the two band solitons oscillate between single vortices, vortex dipoles, and band solitons. The relationship between the spectra of these various collective excitations is an unsolved problem.

### 5.3 Shock and phase rigidity

Explorations of 2D and 3D out of the quasi-1D regime have thus far have been with periodic boundary conditions in the longitudinal direction. If instead *box* boundary conditions are used, together with a phase ramp of the kind desribed in the previous section, the BEC becomes chaotic. As in quasi-1D, it exhibits the property of phase rigidity. In response to shock it breaks up into well-defined phase domains.

To illustrate this effect a volume of 25x12 with box boundary conditions in both dimensions was used. A large phase ramp of $`\varphi (x,y)=1.6\pi x+0.2\pi y`$ pushed a two-soliton, approximate stationary state into the walls at a skewed angle. The time evolution of the resulting state is shown in figure 10. The first time slice shows the initial state. The second time slice shows a shock wave. By the final two time slices the condensate appears to have fractured into a chaotic assortment of phase domains. The phase appears to be completely decoherent at scales larger than the healing length.

## 6 Conclusion

We have analytically described and numerically illustrated the experimental viability of soliton-train stationary states in quasi-one-dimension. Soliton-trains can be manipulated with simple phase profiles. In this way one can create a structured super-current or study soliton interactions.

In the repulsive BEC quasi-1D soliton-train solutions are stable. In the attractive BEC those in which the members of the train are separated by nodes are stable. Other solution-types are cyclically stable. The attractive BEC is stable against collapse in the quasi-1D regime.

Beyond the quasi-1D regime we have shown that band and planar solitons decay into vortices when perturbed by white noise. This is not only a new method of vortex creation but also allows for the direct study of vortex interactions. We have also shown that the phase of the BEC is rigid; when shocked it shatters chaotically.

We benefited greatly from extensive discussions with Charles Clark and Nathan Kutz. This work was supported by NSF Chemistry and Physics and the NSF Reseach Experience for Undergraduates program.

## References


format: text

text:
# 1 Introduction

## 1 Introduction

The study of phase separation and coarsening processes has been a subject of broad and growing interest in recent years . Following a quench into an ordered phase from a disordered one the typical linear size of domains of the ordered phase, $`\mathrm{}(t)`$, grows in time. Usually, at late times the system enters a scaling regime characterized by a single length scale $`\mathrm{}(t)`$ that grows as $`\mathrm{}(t)t^n`$. The value of the exponent $`n`$ depends on the symmetry of the order parameter and on the conservation laws of the system. In two dimensions and above it has been argued that for a scalar order parameter $`n`$ is either $`1/3`$  or $`n=1/2`$  depending on whether or not the dynamics conserves the order parameter, respectively.

In one dimension and for short-range interactions the situation is different due to the absence of long-range order at any temperature $`T>0`$. An exact solution of the non-conserving Glauber dynamics at zero temperature gives the expected scaling $`\mathrm{}(t)t^{1/2}`$ . On the other hand, at zero temperature under Kawasaki dynamics where the order parameter is conserved, the system gets trapped in metastable states with isolated domain walls and coarsening is arrested. Studies of such conserving one-dimensional models in a small $`T`$ limit have shown that the growth law $`\mathrm{}(t)t^{1/3}`$ still persists . However, for both conserving and non-conserving dynamics a solution of the deterministic continuum Ginzburg-Landau equations gives rise to a different answer. In these cases the domain walls interact with each other at large distance through exponentially small forces which in the absence of diffusion leads to the growth law $`\mathrm{}(t)\mathrm{log}(t)`$.

Another situation where a logarithmic growth law arises is when energy barriers proportional to the domain size have to be surmounted during the coarsening. A specific example of this is a three dimensional Ising model with next nearest neighbor interactions . However, in this model the slow coarsening is lost in two dimensions.

Coarsening processes in driven systems are less well understood. Here the dynamics does not obey detailed balance. This gives rise to many phenomena which do not occur in thermal equilibrium . For example, several one-dimensional driven systems have been shown to exhibit long-range order and spontaneous symmetry breaking even when the dynamics is local . A study of the coarsening process in a one-dimensional Ising model with conserved order parameter has shown that for a small driving field the average domain size grows as $`t^{1/2}`$ in contrast to the $`t^{1/3}`$ behavior of the non-driven case . However, in the coarsening of the two-dimensional conserved driven Ising model there is some evidence to suggest that after defining a properly rescaled isotropic domain size the growth law $`\mathrm{}(t)t^{1/3}`$ is retained .

Recently, several one-dimensional driven systems have been shown to exhibit coarsening and phase separation of a novel kind . The phase separation does not rely on microscopic rates tending to zero (as is the case for equilibrium systems in the limit of zero temperature). Instead the phase separation is achieved through a drive which stabilizes certain domain walls and allows ordered domains to be stable if the number of species of domains is greater than two. The mechanism (to be reviewed below) leads to a slow coarsening process in which $`\mathrm{}(t)\mathrm{log}(t)`$ .

In this paper we consider the coarsening dynamics of an extension of one of these models to two and higher dimensions. We show that the slow, logarithmic, coarsening is retained. In the ordered phase the various states are separated by domain walls, which in one dimension are point-like objects. On the other hand in higher dimensions these interfaces are fluctuating extended surfaces. The evolution of the coarsening domains may depend on the smoothness of these surfaces. It would thus be of interest to explore this issue in some detail. An analysis of the model introduced in this paper shows that the interfaces are smooth and have profiles described by a Fermi function with a linear potential (see Eq. 6). This is expected to be a generic feature for the class of models under consideration.

The one dimensional version of the model we study is defined on a ring of length $`L`$ where each site is occupied by one of the three types of particles $`A`$, $`B`$ or $`C`$ . The model evolves under a random sequential update procedure which is defined as follows: at each time step two neighboring sites are chosen randomly, and the particles at these sites are exchanged according to the rates

(1)
The model conserves the number of particles $`N_A,N_B`$ and $`N_C`$ of the three species. When $`q=1`$ the particles undergo symmetric diffusion and the system is disordered. For completeness we now review the explanation of why the system coarsens logarithmically in time for any $`q1`$ . This argument will provide a starting point for the analysis in higher dimensions. To be explicit we consider the case $`q<1`$. Due to the bias, an $`A`$ particle prefers to move to the left inside a $`B`$ domain and to the right inside a $`C`$ domain. Similarly the motion of the $`B`$ and $`C`$ particles in foreign domains is biased. Consider the evolution of the system starting from a random initial condition. The configuration is composed of a random sequence of $`A,B`$ and $`C`$ particles. Due to the bias a local configuration in which an $`A`$ domain is located to the right of a $`B`$ domain is unstable, and the two domains will exchange places on a relatively short time scale linear in the domain size. Similarly, $`AC`$ and $`CB`$ domain walls are unstable. In contrast domain walls of the type $`AB,BC`$ and $`CA`$ are long lived. Thus, after a short time the system will rearrange into a state of the type $`\mathrm{}AABBBCCAAABBCCC\mathrm{}`$ in which only stable domain walls are present. The evolution of this state will proceed by slow diffusion against a bias in which, for example, an $`A`$ particle crosses an adjacent $`C`$ domain. The time scale for such a process to occur is $`q^{\mathrm{}}`$ where $`\mathrm{}`$ is the typical domain size in the system. This suggests that the average domain size in the system grows as $`\mathrm{ln}t/|\mathrm{ln}q|`$. Eventually the system will phase separate into three domains of the three species of the form $`A\mathrm{}AB\mathrm{}BC\mathrm{}C`$ with small density fluctuations around the domain walls, leaving the bulk of the domains pure. In another model in the class this has been referred to as strong phase separation . In the thermodynamic limit these domains do not remix, implying breaking of translational invariance. For the specific case where the number of particles of each species is equal it was shown that the local dynamics obeys detailed balance with respect to a long-range Hamiltonian. Indeed, using this Hamiltonian it was rigorously proved that the model is completely phase separated into pure domains in the steady state .

The model we study is the simplest generalization of the above one-dimensional model to higher dimensionality. To be explicit, we consider the model in $`d=2`$ dimensions but the results presented in this paper are also valid in higher dimensions. We introduce a second lattice direction in which the particles perform unbiased diffusion. That is we study an $`L_x\times L_y`$ lattice, with periodic boundary conditions, in which at each time step two neighboring sites are chosen at random. If the two sites lie along the $`x`$ axis the particles perform a biased diffusion defined by Eq. 1 while if they lie in the $`y`$ axis the particles perform an unbiased diffusion. Thus, along the $`y`$ axis the dynamics is defined through the rates

(2)
Here for simplicity the transition rates are taken to be $`1`$, although choosing a different rate for hops in the $`y`$ direction would not affect the results obtained in this paper. We expect other generalizations of the model, where the motion in the $`y`$ direction is biased in a manner similar to that in the $`x`$ direction, to have the same generic behavior. We return to this point in the discussion.

The model conserves the total number of particles of each species. As will be demonstrated below, in two dimensions the rates do not satisfy detailed balance for $`q1`$ irrespective of the number of particles of each species. Thus the system is generically far from thermal equilibrium.

In the following we show, by numerical and analytical methods, that the slow logarithmic coarsening persists in the two dimensional model and also, we argue, in higher dimensions. In two dimensions the system forms on short time scales stripes of $`A,B`$ and $`C`$ particles which are aligned along the $`y`$ direction and are ordered along the $`x`$ direction in the form $`\mathrm{}AABBBCCAAABBCCC\mathrm{}`$. Thus the interfaces between three phases are lines whose fluctuations must be considered. If the interfaces are smooth the argument suggesting a logarithmic coarsening in the one-dimensional model should also apply in the two-dimensional model. That is, the average tunneling time of a domain would be exponentially small in the average domain size. However, if the interfaces are rough such that they can come close to each other it might be that different behavior would prevail. Thus, an important question is: how rough is the interface? What we find is that the interfaces are indeed smooth and the typical size of a stripe grows as $`\mathrm{}(t)\mathrm{log}(t)`$.

The paper is organized as follows: it is first shown in Section 2 that detailed balance cannot be satisfied in two dimensions (in contrast to the one-dimensional case). In Section 3 the coarsening of the model is studied numerically and a simple interface model which captures the essential physics of an interface between two phases is formulated. The model is solved and shown to have a smooth interface. In Section 4 we discuss the generality of the results both to other two dimensional models and to higher dimensions.

## 2 Lack of detailed balance in the model

To argue that the model presented above does not satisfy detailed balance (except for $`q=1`$) we first review the situation in the one-dimensional version of the model. In this case it was shown  that when the number of particles of each species is equal, that is when $`N_A=N_B=N_C`$, detailed balance is satisfied by the rates given in Eq. 1 for arbitrary $`q`$. However, when the number of particles of two species are not equal, say $`N_AN_B`$, detailed balance is not satisfied. To see this consider an arbitrary set of configurations $`1,2\mathrm{},k`$. Let $`W(ij)`$ be the transition rate from configuration $`i`$ to $`j`$. A necessary and sufficient condition for the existence of detailed balance  is that for any given set of $`k`$ states the following equality is satisfied:

$$W(12)W(23)\mathrm{}W(k1)=W(1k)W(kk1)\mathrm{}W(21).$$
(3)
That is for any closed loop in configuration space the product of rates going along one direction should be equal to the product of rates going in the opposite direction.

We now apply this criterion to show that in $`d=1`$ detailed balance is not satisfied when the three densities are unequal. Consider for simplicity the fully phase separated state $`A\mathrm{}AB\mathrm{}BC\mathrm{}C`$. Take the rightmost $`A`$ particle and move it to the right until it has traversed both the $`B`$ and the $`C`$ domains. The resulting configuration is a fully phase separated state translated by one lattice unit with respect to the starting configuration. Repeating this process $`N`$ times one returns to the starting configuration. The product of the microscopic rates involved in this process is $`q^{NN_B}`$. Carrying out a similar process but in the reverse direction leads again to the starting configuration. For this path the product of rates is $`q^{NN_C}`$. Thus according to the criterion (3) detailed balance is not satisfied when $`N_BN_C`$, or when any two of the densities are not equal. It is not difficult to show that when $`N_A=N_B=N_C`$ detailed balance is satisfied .

Next, consider the two-dimensional model. In this case whatever the number of particles of each species one can always choose as a starting configuration a state where on one of the rows the numbers of particles of each species are not equal. Repeating the argument above, given for the one dimensional model, by considering particle exchanges along this row we immediately see that due to the unequal number of particles detailed balance does not hold.

Therefore, in contrast to the one-dimensional model, whatever that number of particles of each species detailed balance is never satisfied in two dimensions. It is straightforward to see that this argument also implies lack of detailed balance in higher dimensional generalizations of the model.

## 3 Coarsening in two dimensions

To study the dynamics of the model in two dimensions we first show numerically that on a short time scale, the system indeed evolves into a striped state composed of sequences of $`A,B`$ and $`C`$ domains aligned in the $`x`$ direction. We then study the slow coarsening process of these stripes and argue that the typical length scale increases logarithmically with time.

### 3.1 Monte Carlo simulation

In order to study the short time behavior of the model and demonstrate the flow into a striped state, Monte Carlo simulations were performed for various values of $`L_x`$ and $`L_y`$. Starting from a random initial condition we studied the case with equal densities of particles of each species. A typical evolution of a system is presented in Fig. 1 for a lattice with $`L_x=300`$ and $`L_y=80`$ with $`q=0.15`$. One can clearly see that on short time scales a striped structure evolves in which stripes along the $`y`$ direction develop. As expected, the stripes are ordered in sequence $`A,B,C`$ along the $`x`$ direction. On short time scales, when the distance between neighboring domain walls is comparable to their width, the stripes are strongly fluctuating. At this stage the evolution of the structure through the formation of topological defects is relatively fast. As time progresses the asymptotic regime where the stripes are smooth is reached. In this regime the dynamics is dominated by the slow coarsening mechanism described in this paper. Due to the slow time scales involved in the simulations of the coarsening process our data does not allow us to demonstrate quantitatively that the coarsening is logarithmic in time. However in the next section arguments will be presented which support this growth law.

To analyze the data we have calculated the Fourier transform of the density profile along the $`x`$ direction for each value of $`y`$:

$$\mathrm{a}_y(k)=\underset{x=0}{\overset{L_x1}{}}\stackrel{~}{\mathrm{a}}_y(x)\mathrm{exp}\left(\frac{2\pi ikx}{L_x}\right).$$
(4)
Here $`\stackrel{~}{\mathrm{a}}_y(x)`$ is equal to $`1,0`$ or $`1`$ if the site $`(x,y)`$ is occupied by an $`A,B`$ or $`C`$ particle respectively. It thus corresponds to the density difference of the $`A`$ and $`C`$ species. We then calculated $`|\mathrm{a}(k)|`$, the average of $`|\mathrm{a}_y(k)|`$ for all values of $`y`$ and over one hundred simulations. The results are also presented in Fig. 1. One can clearly see that stripes of a characteristic size form. This is accompanied by the emergence of a typical Fourier mode. As time progresses the $`k`$ value of the typical Fourier mode decreases as expected in a coarsening system.

It is of interest to examine the considerations presented above for the existence of phase separation and slow coarsening in one dimension and to check whether the argument can be extended to two and higher dimensions. In $`d`$ dimensions the $`A,B,C`$ domains are separated by fluctuating $`d1`$ dimensional surfaces. For the argument for slow coarsening to go through unmodified, we require that these surfaces are smooth, that is of finite extent in the growth direction. Otherwise, if they were rough, the fluctuating distance between neighboring surfaces might significantly affect the coarsening process. Therefore, we have studied numerically the smoothness of the interfaces in this model. To do this we consider the density of, say, $`A`$ particles near the $`AB`$ domain wall. In Fig. 2 such profiles obtained from Monte Carlo simulations of the three species model are presented. The simulations were performed for a given $`q`$ at fixed value of $`L_x`$ and varying $`L_y`$. They were obtained by extracting the density profile of a given species of particles near a domain wall and averaging over $`40`$ sets of data. One can see that the interface is smooth since the density profile reaches a steady-state form of finite extent in the $`x`$ direction, as $`L_y`$ is increased.

We have thus established numerically that the model evolves towards a striped state and that the interfaces between the phases are smooth. Next, an interface model which captures the essential statistical properties of an interface is studied and shown to predict a smooth interface whose profile agrees very well with the numerics.

### 3.2 Interface model

To construct a simple model that well approximates the behavior of the actual interfaces we begin by assuming that the system is already in a striped configuration in which interfaces are well separated. The density of, say, the $`C`$ particles near an $`AB`$ interface is then negligibly small. This assumption will be shown to be consistent with a smooth interface. Also any current of particles is exponentially small in the domain size and can be neglected, implying that the average velocity of the interface vanishes. We may thus study the interface properties by considering only $`A`$ and $`B`$ particles placed on a lattice with closed (reflecting) boundary conditions in the direction of the bias. As we will show such a model reduces to a model of particles with excluded volume interaction placed in a gravitational-like field in the $`x`$ direction. For such systems the interface is smooth for any $`q1`$. This allows us to derive an exact form for the density profile near an interface which is valid for all dimensions. As will be argued later this result applies to a wide variety of models considered previously and is not specific to the model studied here. For example a similar profile is obtained in  for a one-dimensional model.

We define the interface model on an $`l_x\times l_y`$ lattice with periodic boundary conditions in the $`y`$ direction and closed boundary conditions in the $`x`$ direction. That is, the model is defined on a lattice with closed boundary condition in the direction of the bias. We choose the densities such that half the lattice sites of the lattice are occupied by $`A`$ particles and the others are occupied by $`B`$ particle. Therefore the center (average $`x`$ position) of the $`AB`$ interface is at $`l_x/2`$. The dynamics of the model is identical to that of the three species model (see Eqs. 1 and 2). However, all dynamics in which $`C`$ particles would participate are irrelevant as there are no such particles present.

In this model an $`A`$ can be thought of as a particle while a $`B`$ as a vacancy. Within this picture we have particles, placed on a lattice, which are biased to move in a preferred direction along $`x`$. Due to the closed boundary conditions, which imply a zero current, it is easy to show that the model satisfies detailed balance with respect to the steady state weight (unnormalized probability)

$$W(𝒞)=q^{_{x,y}xA_{x,y}},$$
(5)
where $`A_{x,y}=1(0)`$ if, the site ($`x,y`$) in configuration $`𝒞`$ is occupied (empty), and the sum runs over all lattice sites. The steady state weight describes particles placed in a linear potential with an excluded volume interaction. Note that the above expression is valid whatever the rates for exchange $`A`$ and $`B`$ in the $`y`$ direction.

The solution of the density profile is obtained in straightforward analogy to a Fermi gas. Since each lattice site ($`x,y`$) can be occupied by only a single particle it can be considered as a state of a Fermion with energy $`ϵ_{x,y}=x`$. The mean number of particles at site ($`x,y`$), $`n_{x,y}`$, is given by the Fermi distribution with a temperature $`T=1/|\mathrm{ln}q|`$:

$$n_{x,y}=\frac{1}{q^{(xl_x/2)}+1}.$$
(6)
As expected the distribution is independent of $`y`$. Here $`l_x/2`$ plays the role of a chemical potential marking the center of the profile. It is easy to check the condition that the average density in indeed $`0.5`$, namely $`_{x,y}n_{x,y}=l_xl_y/2`$. To see that this is satisfied for any value of $`q`$ one uses the relation

$$n_{l_x/2+\mathrm{\Delta },y}+n_{l_x/2\mathrm{\Delta },y}=1$$
(7)
where $`\mathrm{\Delta }=xl_x/2`$. The density profile of the interface is given by $`P(x)=n_{x,y}`$. It has a finite width, $`1/|\mathrm{ln}q|`$, as expected for a smooth interface.

A comparison of $`P(x)`$ of the interface model and the profile of the Monte Carlo simulations is given in Fig. 3. One can see that the interface model agrees very well with the simulations.

We have thus established that the interfaces in the model are smooth. This in turn implies that the argument which was presented for the one-dimensional version of the model is valid and predicts a logarithmic growth law for the domain size.

The interface studied above includes configurations with bubbles and overhangs. In the study of interfaces it is often assumed that bubbles and overhangs may be neglected without affecting the macroscopic properties, so that the interface may be described by a single valued function. In the present context we can study the validity of such an approximation since the profile both with (see Eq. 6) and without (see below) bubbles and overhangs can be computed. Thus we consider only those configurations in which the particles fill the left part of the system, up to a ‘height’ $`h_y`$ which is the number of particles in row $`y`$.

By performing the sum in Eq. 5 for each value of $`y`$ one obtains the following weight function:

$$W(𝒞)=q^{\frac{1}{2}_yh_y^2},$$
(8)
where a term of the form $`_yh_y`$ has been neglected, as it is a constant equal to the total number of particles in the system.

The properties of the interface can easily be obtained by working in the grand canonical ensemble. We thus introduce a chemical potential to respect the constraint that the total height of the system $`_yh_y=l_xl_y/2`$. The weight of a configuration in this ensemble is thus

$$W_{gc}(𝒞)=q^{\frac{1}{2}_y(h_yl_x/2)^2}.$$
(9)
The profile of the density along the $`x`$ direction is then given by

$$P(x)=\sqrt{\frac{2|\mathrm{ln}q|}{\pi }}_x^{\mathrm{}}𝑑hq^{\frac{1}{2}(hl_x/2)^2},$$
(10)
where the continuum limit is assumed and the upper limit of the integral was taken to be infinity. This gives

$$P(x)=\frac{1}{2}\mathrm{erfc}\left((xl_x/2)\sqrt{|\mathrm{ln}q|/2}\right),$$
(11)
where $`\mathrm{erfc}`$ is the complementary error function. As expected the model with no overhangs also predicts a smooth interface.

In Fig. 3 the resulting profile is compared with the one obtained by the numerical simulations. It is evident that although the two profiles agree in general, the fit is not as good as that with the Fermi function profile. This is to be expected since, for example, in the regions where the density of particles is low, the no overhangs assumption becomes a less accurate approximation. In these regions excluded volume does not play an important role. Due to the linear potential the profile decays as $`P(x)q^{(xl_x/2)}\mathrm{exp}(\alpha x)`$ with $`\alpha `$ a constant. However neglecting overhangs gives rise to a decay of the form $`P(x)\mathrm{exp}(\beta x^2)`$ with $`\beta `$ a constant.

## 4 Discussion

In this work we have studied a generalization to two dimensions of the $`ABC`$ model for phase separation . We have shown that the slow logarithmic coarsening of the one-dimensional case persists. Further, we have shown that the interface between the slowly coarsening domains is smooth and can be described by a Fermi function with a linear potential (see Eq. 6).

We next turn to the question of the generality of our results. First note that in the generalizations of the $`ABC`$ model to more than three species , and in the models introduced by Lahiri et al  and by Arndt et al  the mechanism which leads to coarsening in the system is the same as in the $`ABC`$ model. Namely, the presence of stable domain boundaries is due to a bias which drives particles of one species towards their own domain with a constant velocity. In the model of Lahiri et al which comprises two coupled rings each containing two species one has four different types of domain. The dynamics is then governed by four stable interfaces separating these domains, in analogy with the model considered in the present work. In the model of Arndt et al there are three types of domain separated by three interfaces. (One of the stable domain walls corresponds to $`q>0`$ in Eq. 6 and the other two correspond to $`q=0`$, i.e. a zero temperature Fermi function.) Therefore, we expect that the analysis of the present paper is applicable to extensions of these models to higher dimensions. This would imply that coarsening processes in these models are logarithmically slow.

The generalization of the model considered so far is highly anisotropic in the sense that the dynamics along the $`x`$ and $`y`$ direction are qualitatively different: while the diffusion is biased in one direction it is symmetric on the other. It is of interest to examine the case where the diffusion along both axes (or all axes in higher dimensions) has some bias. Let us consider such a generalization of the $`ABC`$ model to two dimensions where the exchanges along the $`y`$ direction have rates

(12)
while the exchanges along the $`x`$ axis are as in Eq. 1. On some coarse grained scale (where lattice effects may be ignored) the dynamics is effectively uniaxial. The preferred direction in the $`x`$$`y`$ plane for exchange of particles is at an angle $`\theta =\mathrm{arctan}\left(\frac{1q}{1r}\right)`$ with respect to the $`y`$ direction. The effective exchanges in the normal direction are symmetric. The effect of this would be for stripes to form perpendicular to the preferred direction. At large scales the model reduces to the model which we have studied in this paper. Similarly, generalizations of the $`ABC`$ model in higher dimensions are effectively uniaxial and domains form with interfaces which are $`d1`$ dimensional hyperplanes orthogonal to the preferred direction.

To test this scenario we carried out Monte Carlo simulations of the model in two dimensions with $`r=q`$. In Fig. 4 the coarsening process is illustrated. One sees that after a transient time stripes oriented perpendicular to the $`(1,1)`$ direction are formed confirming the above picture.

Studies of phase separation and coarsening in driven systems such as that of the present paper and the works reviewed in the introduction, have suggested that the processes involved are rather diverse. They range from relatively fast power law to slow logarithmic coarsening. For further examples of coarsening in related systems see . It would be very interesting to construct a more general framework within which all these coarsening processes could be understood and classified.

Acknowledgments: We thank Rashmi Desai for interesting discussions. MRE thanks the Einstein Center for support and DM acknowledges an EPSRC Visiting Fellowship. Support of the Israeli Science Foundation is gratefully acknowledged.


text:
# ROSAT Observations of the Vela Pulsar

## 1. INTRODUCTION

The Vela Pulsar, located at the center of a prominent supernova remnant and with characteristic age ($`p/2\dot{p}`$) of $`1.1\times 10^4`$ years, is “young”. The distance to the pulsar is only 500 pc (Frail and Weisberg, 1990) and the column density of absorbing material ($`2\times 10^{20}`$ atoms cm-2 or less) is low. Thus, soft x-rays from the pulsar can be detected and have been interpreted as thermal radiation from the surface of the hot neutron star. (Harnden et al. 1985, Ögelman et al. 1993).

Like the younger Crab Pulsar, the Vela Pulsar is a source of high energy gamma-rays (Kanbach et al 1994). It is well known that the Crab Pulsar exhibits approximately the same pulse structure at $`\gamma `$-ray, X-ray, and optical wavelengths, two rather narrow pulses, separated by $``$ 0.4 in phase. The Vela Pulsar, however, emits a pulse waveform which varies with energy (Grenier et al. 1988). At 100 MeV, there are two narrow pulses with width and spacing almost identical to those from the Crab Pulsar. In visible light, there are two broad pulses separated by only $``$ 0.2 in phase. X-ray pulsations from the Vela Pulsar were recently discovered by Ögelman et al (1993). They found the pulsed fraction to be $``$ 10% (in contrast to $``$ 100 % for the Crab Pulsar) and the pulse structure, complex.

Like other isolated pulsars, the period of the Vela Pulsar, has been observed to increase steadily and regularly with time. There are, however, in the Vela Pulsar, discontinuities or “glitches” in the timing which occur on average every 2.5 years. To observe this phenomenon, the radio signal is monitored daily (McCulloch et al 1988). During a glitch, the frequency of radio pulsations, $`\nu `$, increases suddenly by $`1`$ part in $`10^6`$ and the frequency derivative, $`\dot{\nu }`$, increases by $``$ 1 part in 102 and then recovers steadily from this transient behavior. Duration of the recovery is 10 - 100 days.

Glitches offer a rare opportunity to learn something about the interior of the neutron star and much has been learned from radio-timing data alone (Link et al 1992). It is believed that, during the glitch, rotational energy is dissipated in the stellar interior. The time scale of energy flow depends on the size of the star, the nature of the glitch mechanism, and the structure of the outer layers of the star. The amount of thermal energy generated and conducted to the surface could be enough to cause an observable increase in temperature of the surface (van Riper et al 1991). We attempted to observe this with ROSAT.

## 2. ROSAT OBSERVATIONS

The requirement that the ROSAT solar panels point at the sun $`\pm 15^{}`$ forbid observations of the Vela Pulsar except during the intervals, 22 April to 26 June and 25 October to 26 December. During these intervals, we obtained data with the ROSAT HRI 6 times. There were also two previous ROSAT PSPC observations. Parameters of the observations are listed in Tables 1 and 2, and exposure times range from 20 to 60 ks.. Figure 1 compares the dates of ROSAT observations with times of recent glitches.

The time structure of the ROSAT observations was not ideal for pulsar work. Typically the target was observed for 1000 - 3000 seconds during each 95-minute orbit but with some long intervals devoted to other targets. A typical data stream for an observation contains $``$ 10 hours exposure spread over $``$ 15 days time.

Figure 1: Time history of Vela Pulsar glitches and ROSAT observations. Observations are below the horizontal line. P indicates PSPC, H indicates HRI. Width of the pip shows time spread of data during that observation. Glitches are above the line. Height of the glitch-line is proportional to $`\mathrm{\Delta }\nu /\nu `$

In the longest exposure (60 ks spread over 14 days) the pulsed signal was strong enough so the pulsar frequency could have been identified using x-ray data alone, even without the known radio-measured pulsar frequency. In the shortest exposure, it was necessary to know the approximate radio-measured pulsar frequency to find the x-ray period.

We searched for the most significant x-ray period by epoch folding over a range bracketing the period calculated from a radio period measured within a few days of the ROSAT observation. In every case, the period calculated from the radio ephemeris agreed within a few ns with the most significant x-ray period in that range. The derived X-ray pulse shape, however, is very dependent on the value used for the period. A shift of a few ns makes a difference. We therefore used the Hartebeesthoek and the U. Tasmania radio data to calculate the expected period to 0.1 ns and used this period to fold the x-ray data. Table 1 lists the input radio data and table 2, the pulse periods calculated (and observed) at the start of each ROSAT observation.

We tried to use the radio observations to calculate the absolute phase of x-ray data but got errors of $``$ 0.2 in phase. To better determine the fine structure of the x-ray pulse, we aligned the HRI observation in phase using the highest peak and summed the observations. 40 phase bins were used. If this procedure were followed using random data, we would, of course, get a single, apparently significant, peak in the summed light curve, but there should not be other significant structure. Since the summed Vela Pulsar light curve shows a second sharp peak, we believe this procedure is valid and the observed structure to be real.

Figure 2 shows the light curve generated by adding the 3 longest pre-Oct-96-glitch observations (May 94, Oct 94, June 96). The pulse structure was clear in each of these and the phase alignment was done to an accuracy of .025 in phase. The same shape is apparent in all the longer individual HRI observations and, is similar to that reported by Ögelman et al. (1993) for the 1991 PSPC observations. There are 2, sometimes 3, narrow peaks superimposed on a broader structure. The Oct 95 and Dec 95 light curves were not as good statistically and the structure not as clear. The appearance of shorter observations is sometimes almost a sawtooth with one prominent peak and valley.

Figure 2: Pre-glitch Vela Pulsar light curve. HRI count rate vs. phase over two complete cycles. The 3 longest pre-glitch HRI observations have been phase-matched using the strongest peak, summed, and fit with an empirical curve. Errors shown are one sigma. Total observation time was 123 ks.

Figure 3: Post-glitch-12 light curve. October 96 data are compared with the pre-glitch empirical fit of Figure 2. Total observation time was 43 ks. This observation started 15 days after glitch 12 and is our best chance, in these data, to observe any post-glitch change in intensity or pulsed waveform. The post-glitch waveform is, within observational limits, the same as the average pre-glitch waveform.

Figure 3 shows the post-glitch October ’96 data alone. This observation started 15 days after glitch 12 and is our best chance, in these data, to observe any post-glitch change in intensity or pulsed waveform. The post-glitch waveform is, within observational limits, the same as the average pre-glitch waveform. The October ’94 observation started 59 days after glitch 11 (which was small) and we have included it in the sum of pre-glitch observations.

An independent technique was used to confirm that pulsed emission was present. The Princeton radio ephemeris was used to calculate the pulsar period at the midpoint of each HRI data set. Data were folded at this period and light curves similar to those illustrated were obtained.

## 3. INTERPRETATION OF PULSED SIGNAL

We assume that the pulse has two components: a broad pulse and two narrow peaks. There is also a large non-varying (probably thermal) component.

We used a $`10^{\prime \prime }`$ radius circular region to extract the pulsar data. The pulsar appeared slightly elongated in all observations because the “wobble” in the ROSAT pointing was not completely removed. This circle was large enough to encompass 85% of the events in the ROSAT point spread function. Table 3 lists the observed components of the HRI signal. Uncertainty in the total rate is due to counting statistics. An order of magnitude larger uncertainty in the signal from the pulsar is due to possible error in subtraction of a contribution from the surrounding diffuse nebula. (Harnden et al. 1985)

To compare pre and post-glitch data we made an empirical fit to the data shown in Figure 2. This fit is indicated by the solid line. Because it appears more square than a pure sine wave, the broad pulse form assumed is $`\mathrm{cos}2\pi \varphi (1/3)\mathrm{cos}6\pi \varphi `$. The two narrow peaks are assumed to be Gaussian with identical widths. The fit to Figure 2 is $`I=0.573+0.029[\mathrm{cos}2\pi \varphi 1/3\mathrm{cos}6\pi \varphi ]+0.071e^{\frac{(\varphi .214)^2}{2(.038)^2}}+0.056e^{\frac{(\varphi .535)^2}{2(.038)^2}}`$. $`\varphi `$ is phase (0-1) and the units of I are ROSAT HRI counts s-1. The two narrow peaks have FWHM of 0.09 in phase and are separated by 0.32 in phase. As summarized in table 3, after subtraction of a nebular component, 12% of the pulsar emission is pulsed with $`8`$% in the broad pulse and $`4`$% in the two narrow peaks.

The 1991 PSPC data were re-analyzed using the same techniques as for the HRI data. A 15′′ radius region was used to extract the pulsar signal and data were again folded at the radio-determined pulsar frequency. The April 1991 PSPC light curve is compatible with the fit to the HRI data. Epoch folding of the December 1991 data produced no clean indication of the pulsar which probably indicates a timing problem in our December 1991 data.

## 4. THE ABSOLUTE RATE

The countrate of the pulsar (and immediate surroundings) was monitored by extracting events from a circle 15′′ in radius, centered on the pulsar. To minimize the effect of variations in detector efficiency and background, the rate from the surrounding annulus with inner and outer radii of 15′′ and 80′′ was also monitored and used for normalization. This annulus contains a relatively bright region of diffuse emission surrounding the pulsar which should not vary on the time scale of this observation. The central region should contain 87% of the pulsar emission (from the ROSAT point response function in David et al 1998) and is large enough to allow for the observed attitude smearing of a few arcseconds.

Figure 4: Pulsar brightness vs time. Count rate from the pulsar and immediate surroundings for the 6 HRI observations. Top curve ($`+`$) is the rate taken from a small region around the pulsar. Middle curve (X) is from a surrounding annulus and includes diffuse emission only. Bottom curve ($`+`$) is the ratio of the two. Horizontal lines are averages weighted by exposure and glitches are indicated as bold vertical lines on the time axis. Error bars show $`1\sigma `$ counting statistics

These data are given in table 4 and plotted in Figure 4. Error bars show $`1\sigma `$ counting statistics only which are $`\pm 1`$ %. No significant variation in the pulsar flux is seen, either in the countrate or normalized countrate. We note that there is a 3% increase in rate after both glitches (or rather a 3% decrease in rate before the glitches) and take this as an upper limit to an observed effect. A 3.0% change in the rate shown in Figure 4 requires a 6.4% change in pulsar flux which, if attributed to black body emission, would correspond to a temperature increase of 1.6%. We also note that the accuracy inferred from Figure 4 is pushing the limit of the ROSAT HRI. Periodic calibrations with the supernova remnant N132D give apparent efficiency variations of up to $`\pm 5`$% for unknown reasons (David et al. 1998). Our normalized rates should avoid this as well as any efficiency change which might be associated with a detector gain (high voltage step) change on June 21, 1994 (David et al. 1998), between the first and second HRI observations.

## 5. A SIMPLE ANALYSIS

It is natural to assume that the narrow pulses are of non-thermal origin. The non-thermal gamma-ray pulses are both sharp and double and the phase separation of the two X-ray pulses is intermediate between the phase separation of the two optical and that of the two gamma-ray pulses (Ramanamurthy, 1994).

The broad pulse could also be non-thermal. The young pulsar, B0540-69 exhibits a single broad pulse at both X-ray and optical wavelengths (Deeter et al., 1999, Gouiffes et al., 1992, Manchester and Peterson, 1989). The young pulsar, B1509-58, also shows a single broad peak in the X-ray band with non-thermal spectrum (Kawai et al 1993). So it is interesting that the Vela pulse waveform shows both narrow Crab-like pulses and the broad pulse associated with two other young pulsars. A thermal origin, however, for either broad or narrow Vela-Pulsar pulses is not ruled out by these data.

Assuming the non-varying component to be thermal, we can derive the temperature and luminosity. (We note that the validity of a black-body approximation to the spectrum is untested and that the derived temperature may not be physical. Nevertheless, the result is useful.) Using a distance of 500 pc, a black body spectrum, a neutron star radiation radius of 10 km, and an interstellar column density of $`5\times 10^{19}`$ atoms/cm2, the PIMMS program was used to calculate the luminosity and surface temperature, $`T`$, from the ROSAT HRI count rate. Results are listed in Table 5. Because the luminosity varies as $`T^4`$, the 20% uncertainty in count rate of the non-varying component, if thermal, corresponds to a 5% uncertainty in temperature. A more realistic range for parameters can be derived by using extremes of the generally accepted ranges for distance (400-500 pc), neutron star radius (7-15 km), and interstellar column ($`N_H=0.22\times 10^{20}`$). Possible values for $`T`$ lie in the range $`511\times 10^5`$ K with corresponding bolometric luminosities of $`28\times 10^{32}`$ erg/s. The most likely fit to our data yields a neutron star temperature of $`8.5\times 10^5`$ K and bolometric luminosity of $`4\times 10^{32}`$ erg/s.

## 6. APPLICATION TO GLITCH PHYSICS

This null result can limit some values of parameters which describe neutron star structure. The parameters are largely unknown, are not independent, and can be linked through models. Theoretical models start with a model neutron star (a choice of mass and equation of state EOS). The remaining parameters include the initial (preglitch) temperature distribution, the amount of energy $`\mathrm{\Delta }`$E dissipated in the glitch and the location, or density $`\rho `$, where the energy dissipation takes place. This heat, released promptly during the glitch, will diffuse out to the neutron star surface and inward, into the neutron star core. To calculate the postglitch thermal signal on the neutron star surface, models of the heat capacity and thermal conductivity, as well as neutrino emissivities from the interior are required. At the age of the Vela pulsar, the interior of the neutron star is expected to be isothermal with a core temperature extending through the inner crust. The associated surface temperature is determined by conductivities in the outer crust (Gudmundsson, Pethick and Epstein 1981). Calculations of the expected postglitch thermal signal on the surface have been reported by Eichler and Cheng (1989), van Riper, Epstein and Miller (1991), by Chong and Cheng (1994) and by Hirano et al (1996) for various combinations of the parameters. The results are qualitatively similar. All calculations are quite sensitive to the initial value of the surface temperature. The fractional change in temperature is higher for lower initial temperatures. For a conservative evaluation of our null results, in terms of which model parameters can be ruled out, we take the initial surface temperature to be the highest value, covered in the model calculations, 106 K. This is also the appropriate choice of initial surface temperature for the Vela pulsar. The ROSAT measurements of Ögelman, Finley and Zimmermann (1993) yield surface temperatures of 1.5-1.6 $`\times `$ 106 K from black body fits to the point source and to the pulsed signal. The model-dependent analysis of the previous section yields a temperature of 0.85 $`\times `$ 106 K. The actual value of the surface temperature may be somewhat different but of the same order when neutron star atmosphere models are used.

Van Riper et al (1991) present their results in terms of tpeak, the time at which the glitch-induced temperature signal on the surface is maximum; the width $`\mathrm{\Delta }`$t of this signal and the maximum fractional change in the surface temperature. Contours of tpeak and $`\mathrm{\Delta }`$t for which the postglitch temperature enhancement signal would fall within the range of our observations, 15-22 days after the glitch, correspond roughly to neutron star radii R $`<`$ 13 Km. Thus for Ts,0 = 106 K our upper limits test the models only for stars with R $`<`$ 13 Km. The radius of a typical 1.4 M neutron star is 13.7 Km in the soft Baym, Pethick and Sutherland (1971) equation of state and 13.7 Km for the medium equation of state of Friedman and Pandharipande (1981). For a stiff equation of state like the Pandharipande, Pines and Smith (1976) model the radius of a 1.4 M star is 18.6 Km, so our results do not test the models if the neutron star EOS is actually stiff. If the energy release $`\mathrm{\Delta }`$E in the glitch was 1042 ergs, the thermal signal would violate our upper bound if the glitch dissipated energy at densities $`\rho 10^{13}`$ gm cm-3 or less so it is likely that $`\mathrm{\Delta }`$E $`<`$ 1042 ergs, conditional on the $`\rho `$ values. $`\mathrm{\Delta }`$E = 1043 ergs can be ruled out altogether, for no matter at what $`\rho `$ the glitch energy was dissipated, the signal on the surface would violate our upper bounds.

Chong and Cheng (1994) plot their results for 1.4 M neutron stars of different EOS, initial surface temperatures and $`\mathrm{\Delta }`$E , taking the density at which the glitch energy is dissipated to be $`\rho `$ = 1013 gm cm-3 for all cases. Taking their model for a soft BPS star with the core temperature Tc = 108 K, corresponding to an initial surface temperature of 1.27 $`\times `$ 106 K we find that $`\mathrm{\Delta }`$E = 1043 ergs is ruled out by our results. For the moderately stiff UT equation of state (Wiringa & Fiks 1988) with Tc = 108 K and Ts,0 = 106 K, the model prediction for the temperature enhancement at the times of our observations lie below the upper limits, so even $`\mathrm{\Delta }`$E = 1043 ergs is not ruled out. For the stiff PPS star model, with Tc = 108 K and Ts,0 = 8 $`\times `$ 105 K, $`\mathrm{\Delta }`$E = 1043 ergs is ruled out by our upper limits.

Hirano et al (1997) present similar results by marking in the AT $``$ $`\mathrm{\Delta }`$ Ts/Ts,0 vs tpeak plane the maximum AT points for different 1.4 M neutron star models and different densities $`\rho `$ of energy release. Results are given for various values of the initial surface temperature Ts,0 and energy release $`\mathrm{\Delta }`$E. To check if the predicted signals extend into our range of observation times, we utilized the temperature evolution for one-dimensional heat conduction, Eq.(4) of Hirano et al. which they quote as representative of the temperature evolution in their detailed numerical results. At Ts,0 = 106 K, $`\mathrm{\Delta }`$E = 1043 ergs is ruled out if the EOS is soft (the Baym, Pethick, Sutherland (1971) model). For a moderately stiff FP star, our results rule out $`\mathrm{\Delta }`$E = 1043 ergs if the energy was released at $`\rho `$ $`<`$ 3 $`\times `$ 1013 gm cm-3.

In addition to these model calculations, all of which assume that the glitch induced energy dissipation takes place in a spherically symmetric shell, a recent paper by Cheng, Li and Suen (1998) presents model calculations of the non-spherically symmetric case. These authors calculate the thermal signal on the neutron star surface for the case of glitch induced energy dissipation at densities 3 $`\times `$ 1012 gm cm-3 $`<`$ $`\rho `$ $`<`$ 3 $`\times `$ 1013 gm cm-3, and within a 2 deg $`\times `$ 2 deg solid angle. For $`\mathrm{\Delta }`$E = 1042 ergs a hot spot is found to emerge on the surface, keeping to the same solid angle range and peaking at 274 days after the glitch. During our observations, 15-22 days after the glitch, a 50 percent enhancement in temperature is predicted in the hot spot. The contribution to the total luminosity of this hot spot remains well below our upper limits. The absence of significant changes in the pulse shape of the Vela pulsar before and after a glitch does not lead to a test of the models. The hot spot is at the same angular position as the region of energy dissipation in the glitch. For many values of the angles from the rotation axis and from our line of sight, there would be no modulation in the observed signal. Furthermore, the energy dissipation might take place in an equatorial belt around the rotation axis, as is likely if the glitch is due to vortex line unpinning (a rotational instability in the neutron star crust superfluid).

Returning to the spherically symmetric models we note that the differences between the model predictions are not important for purposes of comparison with our current upper limits. (These differences are likely to be due to different surface layer compositions in the models in terms of the A and Z of the equilibrium nuclear species, effecting the thermal conductivities (Hirano et al 1997)). Within the uncertainties indicated by these differences, and the large number of unknown parameters, we conclude that our upper bounds on a thermal signal rule out glitches with energy release, $`\mathrm{\Delta }`$E, greater than 1043 ergs.

For superfluid unpinning models of the glitches the energy dissipated is given by:

$$\mathrm{\Delta }E=I_p\delta \mathrm{\Omega }\omega =I\mathrm{\Delta }\mathrm{\Omega }\omega .$$
(1)
In these models the glitch reflects a transfer of angular momentum, J= $`I_p\delta \mathrm{\Omega }=I\mathrm{\Delta }\mathrm{\Omega }`$, from a component of the neutron star, taken to be the pinned superfluid in the inner crust, to the rest of the star including the observed crust. Ip is the moment of inertia of the pinned superfluid and $`\delta \mathrm{\Omega }`$ the decrease in its rotation rate during the glitch. I is the moment of inertia of the observed crust and most of the neutron star, which is coupled rigidly to the crust. $`\mathrm{\Delta }\mathrm{\Omega }`$ is the observed jump in rotation frequency of the crust and $`\omega `$ is the difference in rotation rate between the superfluid and the normal crust. $`\omega `$ is related to the pinning forces. Our upper bound on $`\mathrm{\Delta }E`$ implies a very loose upper bound $`\omega `$ $`<`$ 100 rad s-1. Recent theoretical work (Pizzochero, Viverit & Broglia 1997) gives pinning energies which imply $`\omega `$ 1 rad s-1. The coupling between the pinned superfluid and the observed crust of the neutron star also entails a continuous rate of energy dissipation, realized though the thermal creep of vortices. In parallel to the glitch associated energy dissipation, this continuous rate is easily shown to be (Alpar et al 1984):

$$\dot{E}_{diss}=I_p|\dot{\mathrm{\Omega }}|\omega .$$
(2)
This should supply the thermal luminosity of an old pulsar. Observational upper bounds from the thermal luminosities of older pulsars (Alpar et al 1987, Yancopoulos, Hamilton & Helfand 1994) yield $`\omega `$ $`<`$ 1 rad s-1.

While crustquake models fail to explain the magnitude and rate of glitches of the Vela pulsar, crustquakes may trigger the vortex unpinning and angular momentum transfer that show up as glitches. The crust cracking that triggers a glitch (Ruderman 1991) entails an energy dissipation

$$\mathrm{\Delta }E=\mu \theta ^2$$
(3)
where $`\mu `$ denotes the elastic energy modulus and $`\theta `$ is the critical strain angle for the breaking of the lattice. For a volume of linear dimensions 105 cm $`\mu `$ = 1044 ergs. Our upper limit of $`\mathrm{\Delta }E`$= 1043 ergs then implies that $`\theta `$ $`<`$ 10-1 rad s-1. While terrestrial solids have $`\theta 10^4`$ or less, there have been speculations that the neutron star crust will have $`\theta `$ $`>`$ 10-2, a much larger strength resulting from its unscreened Coulomb interactions. Thus, observations that are sensitive to implied glitch energy dissipation rates $`\mathrm{\Delta }E`$ = 1041 ergs will be critical for the crust-breaking models of glitches as well as for the superfluid-unpinning model.

While the present upper limit does not give any stringent constraints on the glitch related energy release, a future detection could yield important clues into the neutron star EOS, since in all calculations, different EOS yield different timescales for the postglitch thermal signal. Post-glitch RXTE-PCA observations of the Vela pulsar, at one, four and nine days, and at about three months following the October 1996 glitch do not constrain models of glitch related energy dissipation (Gürkan et al. 1999). New instruments should have the capability to detect temperature rises from the Vela pulsar at the corresponding levels of sensitivity. The Chandra HRC detector is 3.4 times as sensitive as the ROSAT HRI and can thus collect data equivalent to that shown in Figure 2 in 36 ks. Two such observations, one day and one week after a glitch, should measure an increase in surface luminosity if the energy release is 1041 ergs or greater. This time frame does not cover all of parameter space but is probably the best place to search.

This work was supported by NASA through contract NAS8-39073 and grant NAG5-6853.

| Table 1 | | | |
| --- | --- | --- | --- |
|  | Rotation | Frequency |  |
| Date | Frequency | Derivative | Epoch |
| (UT) | (Hz) | $`10^{12}`$ (Hz s-1) | (MJD) |
| 14-26 May 94 | 11.197 366 280 36 | -15.5881 | 49492.274659780 |
| 25-28 Oct 94 | 11.197 160 398 75 | -15.6554 | 49653.5062011504 |
| 27 Oct -2 Nov 95 | 11.196 668 608 990 | -15.5921 | 50017.020877910 |
| 15-16 Dec 95 | 11.196 603 154 274 | -15.5858 | 50066.015970352 |
| 1-3 Jun 96 | 11.196 376 285 216 | -15.5777 | 50236.139720991 |
| 28 Oct - 4 Nov 96 | 11.196 198 755 884 | -15.5693 | 50385.205072337 |

| Table 2 | | | | | |
| --- | --- | --- | --- | --- | --- |
|  |  | Exposure | X-ray | Period |  |
| Date | Instrument | Live time | Period | Derivative | Epoch |
| (UT) |  | (s) | (s) | $`10^{13}`$(s s-1) | (MJD) |
| 14-26 May 94 | HRI | 62090 | 0.089 306 649 49 | 1.2433 | 49486.169164815 |
| 25-28 Oct 94 | HRI | 31545 | 0.089 308 327 69 | 1.2486 | 49650.775217647 |
| 27 Oct - 2 Nov 95 | HRI | 19409 | 0.089 312 279 83 | 1.2437 | 50017.896366423 |
| 15-16 Dec 95 | HRI | 20380 | 0.089 312 801 95 | 1.2432 | 50066.494671817 |
| 1-3 Jun 96 | HRI | 28834 | 0.089 314 611 67 | 1.2427 | 50235.015672685 |
| 28 Oct - 4 Nov 96 | HRI | 43194 | 0.089 316 027 86 | 1.2520 | 50384.625402581 |

| Table 3 | |
| --- | --- |
| Total rate within $`10^{\prime \prime }`$ radius circle (s-1) | 0.585 $`\pm `$ .005 |
| non X-ray and diffuse background | 0.0007 $`\pm `$ .0001 |
| contribution from surrounding nebula | 0.258 $`\pm `$ .066 |
| total signal from Vela PSR | 0.326 $`\pm `$ .066 |
| 2 narrow pulses | 0.0122 $`\pm `$ .003 |
| broad-pulse | 0.0277 $`\pm `$ .007 |
| steady (thermal?) | 0.286 $`\pm `$ .066 |

| Table 4 | | |
| --- | --- | --- |
| Date | Average Rate | Normalized Rate |
| May 94 | 0.707$`\pm `$.004 | 1.043$`\pm `$.008 |
| Oct 94 | 0.723$`\pm `$.005 | 1.073$`\pm `$.010 |
| Oct 95 | 0.709$`\pm `$.006 | 1.078$`\pm `$.011 |
| Dec 95 | 0.722$`\pm `$.006 | 1.093$`\pm `$.011 |
| Jun 96 | 0.699$`\pm `$.005 | 1.040$`\pm `$.010 |
| Oct 96 | 0.718$`\pm `$.004 | 1.067$`\pm `$.008 |

| Table 5 | | | | | |
| --- | --- | --- | --- | --- | --- |
|  | ISM | Neutron Star |  |  |  |
| Distance | Column | Radius | kT | T | Lbol |
| (pc) | ($`10^{20}\text{cm}^2`$) | (km) | (keV) | ($`10^5`$K) | ($`10^{32}\text{erg/s}`$) |
| 500 | 0.5 | 7 | .082 | 9.5 | 3.7 |
| 500 | 0.5 | 10 | .070 | 8.5 | 3.8 |
| 500 | 0.5 | 15 | .055 | 6.4 | 4.1 |
| 500 | 0.2 | 10 | .066 | 7.7 | 3.0 |
| 500 | 2 | 10 | .082 | 9.5 | 7.0 |
| extremes |  |  |  |  |  |
| 500 | 2 | 7 | .095 | 11 | 8.2 |
| 400 | 0.2 | 15 | .047 | 5.5 | 1.9 |


text:
# THE SCALED UNIVERSE II

## 1 Introduction

In a previous communication it was shown that the mysterious quantum prescription of microphysics has analogues at the much larger scale of stars, galaxies and superclusters. The common feature in all these cases is a Brownian type fractality. We now examine this circumstance in greater detail to obtain a rationale for this ”Scaled Quantum Mechanical” behaviour.
At the outset it must be observed that in a sense, the universe is not a continuum. Thus within an atom, the nucleus occupies a very tiny fraction of the volume, roughly $`10^{15}`$. Then there is a wide gap till we reach the orbiting electrons. Similarly there are intermolecular distances, interplanetary, interstellar, intergalactic…. distances which provide relatively huge gaps. This is not in the spirit of a uniform continuum.
We will now argue that this is because, for example the nucleons are bound together, so also the electrons and the nucleons are bound together, the atoms in the molecules are bound together… and so on with subsequent gaps, which leads to some interesting scale dependent consequences, all this in the context of a Brownian underpinning.

## 2 A Rationale for Scaled Effects

The starting point in was the fact that we have the following Random Walk type Relations:

$$Rl_1\sqrt{N_1}$$
(1)
$$Rl_2\sqrt{N_2}$$
(2)
$$l_2l_3\sqrt{N_3}$$
(3)
$$Rl\sqrt{N}$$
(4)
where $`N_110^6`$ is the number of superclusters in the universe, $`l_110^{25}cms`$ is a typical supercluster size $`N_210^{11}`$ is the number of galaxies in the universe and $`l_210^{23}cms`$ is the typical size of a galaxy, $`l_31`$ light year is a typical distance between stars and $`N_310^{11}`$ is the number of stars in a galaxy, $`R`$ being the radius of the universe $`10^{28}cms,N10^{80}`$ is the number of elementary particles, typically pions in the universe and $`l`$ is the pion Compton wavelength.
On this basis it was argued that we could introduce a ”Scaled” Planck Constant given by

$$h_110^{93}$$
(5)
for super clusters;

$$h_210^{74}$$
(6)
for galaxies and

$$h_310^{54}$$
(7)
for stars.
The relation (4) was observed nearly a century ago by Weyl and Eddington. It was shown (Cf.ref.) that far from being empirical this relation can be deduced on the basis of the fluctuational creation of particles from a background Zero Point Field or Quantum Vacuum, in a scheme which leads to a cosmology consistent also with Dirac’s large number coincidences. From this point of view the Random Walk character of equation (4) is not accidental, and this reasoning could be extended to equations (1),(2) and (3), in the light of equations (5), (6) and (7).
It may be observed that in all these cases we have a length, the Compton wavelength or its analogue which defines regions of matter separating relatively empty spaces.
Infact it was also argued in that these scaled ”Compton wavelengths” and scaled ”Planck Constants” arise due to the well known equation of gravitational orbits,

$$\frac{GM}{L}v^2$$
(8)
On the other hand equation (8) can be viewed as resulting from the Virial Theorem, where the velocity is replaced by the velocity dispersion.
This velocity $`v`$ would be different at different scales. For example for a black hole it would be the velocity of light, which would then give the Schwarzchild radius from (8). For galaxies it is $`10^7cms`$ per second. It is this circumstance that produces the above scales leading to fractality.
We could go one step further, because we expect that the same effect would apply to solar type systems: The planets and other objects are bound quite close to the sun compared to the interstellar distances. Infact we can verify that this is so for Kuiper Belt objects which have been studied in the recent past. In this case a typical size is $`5km`$, the distances are $`10^{15}cm`$, masses are $`10^{19}gms`$ their number is $`10^{10}`$ while an application of equation (8) shows that the velocities are $`10^5cms`$ per second which is indeed so. It can now be shown quite easily that this defines a scaled Planck Constant

$$h_410^{34}$$

In other words equation (8) characterizes ”Black Holes” with different maximal velocities.
Incidentally from (8) we could easily deduce that the angular momentum $`J`$ is given by

$$JM^2$$
(9)
It is quite remarkable that the equation (9) also applies to elementary particles and Regge trajectories.
Indeed this is no coincidence as can be seen as follows: It was shown that a pion can be considered to be a bound state of an electron and a positron within the context of Quantum Mechanical Black Holes. In this case we have instead of equation (8)

$$c^2\frac{e^2}{m_er}$$
(10)
where $`m_e`$ is the electron mass.
Whence we have

$$hm_\pi cr\frac{e^2\times 10^2}{c}10^{27}$$
(11)
It is remarkable that (11) gives us the value of Planck’s constant, given the pion and electron masses, and at the same time we can deduce an equation like (9) at the micro level also. Equation (11) is on the same footing as equations (5), (6) and (7). In this case it is the electromagnetic interaction which replaces the gravitational interaction in the equation (8).
Thus the rationale for the fractal structure seen above is in bound systems separated by relatively large and relatively empty spaces.