text:
# Transmittivity of a Bose-Einstein condensate on a lattice: interference from period doubling and the effect of disorder

## 1 Introduction

A quasi-onedimensional (1D) optical lattice created by the interference of two counterpropagating laser beams and superposed on a highly elongated magnetic trap provides an almost ideal periodic potential for an atomic gas. A mixture of a Bose-Einstein condensate and a degenerate fermion gas is an interesting system to study in such a confinement. The interactions can result in the spatial phase separation of the two components  or in the collapse of the fermionic component  and generally induce strong changes in the atomic dynamics , depending on the values of the boson-boson and boson-fermion coupling constants.

A number of experiments have examined properties of ultracold bosonic atoms and condensates in optical lattices: Bloch oscillations , Landau-Zener tunnelling , squeezed states , Josephson-like oscillations , a superfluid-to-Mott insulator transition , and 1D band structure  have all been observed. Theoretically, these effects have been studied by a variety of approaches such as the Gross-Pitaevski equation , the discretised non-linear Schrödinger equation , the Bogoliubov-de Gennes equations , the direct diagonalization of the Bose-Hubbard Hamiltonian , and quantum MonteCarlo simulations .

In this work we study the transport properties of the bosons in a boson-fermion mixture at zero temperature under the action of a constant driving force inside a 1D lattice. A constant acceleration can be imparted to the boson condensate by using a frequency-chirped optical lattice . For a suitable choice of the system parameters the tunnelling of the fermionic atoms across the lattice can be neglected and their presence simply modifies the well depths seen by the bosons. The interesting cases are those in which the fermionic component introduces additional periodicities or generates disorder in the lattice potential. All these modifications in the lattice potential could be more easily realized in the laboratory by using extra laser beams and/or speckle . However, we choose to focus on more sofisticated experimental situations where the additional periodicity or the disorder is induced by the fermions, with the aim to study the magnitude of these effects in such a system and to address our numerical results to experimental groups who investigate quantum transport in boson-fermion mixtures.

We describe the system by means of a 1D Bose-Hubbard Hamiltonian, using a renormalization of the scattering lengths to embody the transverse confinement . The Hamiltonian is not explicitly solved, but the whole lattice is reduced to a single dimer by a renormalization/decimation procedure . Incoming and outgoing leads are connected to the effective dimer and the steady-state transport coefficient of the bosons is inferred from the evaluation of the scattered wave function of the leads due to the presence of the dimer.

The paper is organized as follows. In Sec. 2 we present the physical system of specific interest and the model that we adopt to evaluate on-site and hopping energies of the 1D Bose-Hubbard Hamiltonian. In Sec. 3 we introduce a scattering matrix formalism for the case of out-of-equilibrium leads, with the aim of calculating the bosonic transmittivity through the system. The particle current is obtained from an extrapolation using a Landauer-like formula. The method is applied for two choices of the number of fermions, in cases when they are either regularly or randomly distributed inside the lattice. The numerical results and their physical interpretation are discussed in Sec. 4, and the main conclusions are given in Sec. 5. Two Appendices give technical details on the evaluation of some parameters entering the effective 1D Hamiltonian and on the recursive algorithm for the calculation of the scattering matrix, while a third Appendix discusses how period doubling for a pure Bose condensate can be achieved in a four-beam experimental set-up.

## 2 The model

We focus on a mixture of 87Rb and 40K atoms, which has the peculiarity of having a large and negative boson-fermion scattering length. Let us consider the case in which the number $`N_f`$ of (spin-polarized) fermions is much less than the number $`N_b`$ of bosons and at most equal to the number $`n_s`$ of available lattice sites. With a suitable choice of the laser frequency $`\omega _L`$, the lattice potential can become considerably deeper for the 40K atoms than for the 87Rb atoms, so that the fermions can be taken as localized inside the lattice potential wells.

The lattice potential for the two atomic species can be written as

$$U_{\mathrm{K},\mathrm{Rb}}(z)=U_{\mathrm{K},\mathrm{Rb}}^0\mathrm{sin}^2(kz)\frac{\mathrm{}\mathrm{\Omega }_{\mathrm{K},\mathrm{Rb}}^2}{4\delta _{\mathrm{K},\mathrm{Rb}}}\mathrm{sin}^2(kz)$$
(1)
where $`\delta _{\mathrm{K},\mathrm{Rb}}=\omega _L\omega _{\mathrm{K},\mathrm{Rb}}`$ is the detuning of the laser from the two atomic frequencies $`\omega _\mathrm{K}`$ and $`\omega _{\mathrm{Rb}}`$, $`k=2\pi /\lambda `$ is the wave number of the laser light and $`\mathrm{\Omega }_{\mathrm{K},\mathrm{Rb}}=d_{\mathrm{K},\mathrm{Rb}}E_0/\mathrm{}`$ is the Rabi frequency, which is determined by the laser electric field $`E_0`$ and by the atomic dipole $`d_{\mathrm{K},\mathrm{Rb}}`$. The period of the optical lattice is $`\lambda /2`$. The atomic dipole depends on the natural width $`\mathrm{\Gamma }_{\mathrm{K},\mathrm{Rb}}`$ of the transition  so that the ratio $`U_\mathrm{K}^0/U_{\mathrm{Rb}}^0`$ is given by

$$\frac{U_\mathrm{K}^0}{U_{\mathrm{Rb}}^0}=\frac{\mathrm{\Gamma }_\mathrm{K}}{\mathrm{\Gamma }_{\mathrm{Rb}}}\frac{\delta _{\mathrm{Rb}}}{\delta _\mathrm{K}}.$$
(2)
For the D2 lines of 40K and 87Rb, which are used in current experiments , the natural widths are almost equal ($`\mathrm{\Gamma }_\mathrm{K}\mathrm{\Gamma }_{\mathrm{Rb}}6`$MHz), while the atomic frequencies have values $`\omega _\mathrm{K}=390.83`$ THz and $`\omega _{\mathrm{Rb}}=384.227`$ THz. Choosing for example a laser wavelength $`\lambda =763`$ nm (or $`\omega _L=392.88`$ THz) the ratio between the depths of the two lattice potentials is $`U_\mathrm{K}^0/U_{\mathrm{Rb}}^0=4.26`$, so that the transport of the fermions can be neglected in a regime where tunnelling motions of the bosons are still allowed.

However, the presence of the localized fermions affects the transmittivity of the bosons by modifying the effective potential that the latter see. In the case of a negative mutual scattering length a well containing a fermion becomes deeper for the bosons inside it, so that these bosons will have a lower probability of hopping to the next well.

To study the effect of the fermions on the bosonic transport as indicated above, we build a 1D tight-binding Hamiltonian for the bosons and use a Green’s function approach to evaluate their transport coefficient through the lattice. The Bose-Hubbard Hamiltonian for the bosons is

$$H_I=\underset{i=1}{\overset{n_s}{}}\left[E_i|ii|+\gamma _i(|ii+1|+|i+1i|)\right].$$
(3)
Here, the parameters $`E_i`$ and $`\gamma _i`$ depend on the number of bosons in the lattice well labelled by the index $`i`$ and represent site energies and hopping energies, respectively.

We proceed to a 1D reduction of the Hamiltonian by introducing the transverse widths $`\sigma _{\mathrm{Rb},\mathrm{K}}`$ of the bosonic and fermionic wave functions in a cigar-shaped harmonic trap, together with a 1D condensate wave function $`\varphi _i(z)`$ in the $`ith`$ cell and the 1D fermion density $`n_f(z)`$. In a tight-binding scheme $`\varphi _i(z)`$ is a Wannier function for the bosons in the potential $`U_{\mathrm{Rb}}(z)`$ given in Eq.(1) and, according to the early work of Slater , can be written as a Gaussian function,

$$\varphi _i(z)=\varphi _i(0)\mathrm{exp}[(zz_i)^2/(2\sigma _{z\mathrm{Rb}}^2)],$$
(4)
where $`|\varphi _i(0)|^2`$ is the number of bosons in the lattice well $`i`$. The fermion density is similarly written as a superposition of Gaussian functions localized in a set of lattice wells labelled by an index $`i^{}`$,

$$n_f(z)\underset{i^{}}{}\mathrm{exp}[(zz_i^{})^2/\sigma _{z\mathrm{K}}^2].$$
(5)
The determination of the widths $`\sigma _{\mathrm{Rb},\mathrm{K}}`$ and $`\sigma _{z\mathrm{Rb},\mathrm{K}}`$ is carried out variationally, as described in A.

We can now evaluate the parameters entering the effective Hamiltonian. The site energies are given by

$$E_i=𝑑z\varphi _i^2(z)\left[\frac{\mathrm{}^2^2}{2m_{\mathrm{Rb}}}+U_{\mathrm{Rb}}(z)+\frac{1}{2}g_{bb}|\varphi _i(z)|^2+g_{bf}n_f(z)m_{\mathrm{Rb}}az+C_{\mathrm{Rb}}\right]$$
(6)
where $`m_{\mathrm{Rb}}`$ is the 87Rb mass, $`a=F/m_{\mathrm{Rb}}`$ is the acceleration due to a constant external force $`F`$ acting on the bosons, $`C_{\mathrm{Rb}}=\mathrm{}^2/(2m_{\mathrm{Rb}}\sigma _{\mathrm{Rb}}^2)+\frac{1}{2}m_{\mathrm{Rb}}\omega _{\mathrm{Rb}}^2\sigma _{\mathrm{Rb}}^2`$, and $`\omega _{\mathrm{Rb}}`$ is the radial frequency of the harmonic trapping potential. The parameters $`g_{bb}`$ and $`g_{bf}`$ are the strengths of the 1D boson-boson and boson-fermion interactions, which are given by

$$\{\begin{array}{c}g_{bb}=\frac{4\pi \mathrm{}^2}{m_{\mathrm{Rb}}}\frac{a_{bb}}{2\pi \sigma _{\mathrm{Rb}}^2}\hfill \\ g_{bf}=\frac{2\pi \mathrm{}^2}{m_r}\frac{a_{bf}}{\pi \left(\sigma _{\mathrm{Rb}}^2+\sigma _\mathrm{K}^2\right)}\hfill \end{array}$$
(7)
with $`a_{bb}`$, $`a_{bf}`$ the boson-boson and boson-fermion scattering lengths  and $`m_r`$ the boson-fermion reduced mass. Consistently with the model, the hopping energies $`\gamma _i`$ are given by

$$\gamma _i=𝑑z\varphi _i(z)\left[\frac{\mathrm{}^2^2}{2m_{\mathrm{Rb}}}+U_{\mathrm{Rb}}(z)+\frac{1}{2}g_{bb}|\varphi _i(z)|^2+g_{bf}n_f(z)+C_{\mathrm{Rb}}\right]\varphi _{i+1}(z).$$
(8)
This completes the determination of the effective 1D Hamiltonian for the bosons, that we shall use to evaluate the transmittivity for various distributions of the fermions in the lattice wells.

## 3 The transmittivity

In Fig. 1 we show a schematic representation of the Hamiltonian (3) in the absence of an external force. For such a system the boson site energies can be of two types, $`E_1`$ and $`E_2`$ depending on the absence or presence of a fermion in the well, and the boson hopping energies can take three values: $`\gamma _1`$ if the hopping is between a site with a fermion and a site without, $`\gamma _2`$ in the case of two neighboring wells which both contain a fermion, and $`\gamma _3`$ in the case of two neighboring wells without fermions. We shall consider the following situations: (i) the fermionic band is empty (EF), namely no fermions are present in the lattice; (ii) the fermionic band is filled (FF), namely there is one fermion per well; (iii) the fermionic band is half-filled ($`N_f=n_s/2`$) with every other well occupied in an orderly way (OHF); and (iv) the fermionic band is half-filled but the fermions are distributed on the lattice in a disordered way (DHF).

In the tight-binding regime the system of bosons behaves in the first two cases as a one-band system, which is described by a single site energy ($`E_1`$ in the first case and $`E_2`$ in the second) and by a single hopping energy ($`\gamma _3`$ and $`\gamma _2`$, respectively). The third case describes a two-band boson system with two alternating site energies ($`E_1`$ and $`E_2`$) and a single hopping energy ($`\gamma _1`$). There is in this case a folding of the dispersion relation in quasi-momentum space, with the opening of an energy minigap $`\mathrm{\Delta }E=|E_2E_1|`$ between two sub-bands of equal energy width. In the fourth case the bosons again form a two-band system of the type illustrated in Fig. 1 and the presence of disorder will cause a broadening of the two bands. We show in B how the third case could be realized in the laboratory by using a four-beam set-up instead of an ordered array of fermionic atoms at half occupancy. The forth case could be realized by adding some speckle to the lattice .

The effect of applying a constant force to the condensate can be described as a tilting of the bands in coordinate space, due to the induced position dependence of the site energies. This is shown in Fig. 2: the situation is analogous to that obtained in an electron gas on a lattice by applying to it a static, uniform electric field . The transmittivity of such a system can be calculated, following the proposal of Farchioni et al. , through an evaluation of the scattering matrix which explicitly takes into account the presence of the bias . For this purpose, let us connect the system to incoming and outgoing leads which mimic its coupling to the continuum by injecting and extracting a steady-state particle current, respectively. The bosons are initially in their lowest momentum state ($`p=0`$) and can leave the lattice only at $`p=h/\lambda `$ (in the one-band case) or at $`p=0`$ (in a small-minigap two-band system) if the coupling to the leads is effected as is shown in Fig. 2. Elastic Bragg scattering can also occur at both these points in quasi-momentum space as well as at the top of the lower sub-band in the split-band configuration. As we shall see, in the latter case this produces resonances in the transmittivity.

Let us define $`N_{\mathrm{out}}`$ (with $`N_{\mathrm{out}}n_s`$) as the number of hops after which the condensate reaches the highest energy point in the dispersion curve and can either leave the lattice or be Bragg reflected. The two leads have a band width $`4\gamma _L`$ and their centers at $`E_L^{\mathrm{in}}`$ and $`E_L^{\mathrm{out}}`$ are shifted by $`FN_{\mathrm{out}}\lambda /2`$ in order to optimize the coupling with the lattice. The outgoing wavefunction $`\varphi _{\mathrm{out}}^{\stackrel{~}{\kappa }}(z)`$ is an eigenstate of the outgoing lead and is connected to the incoming wavefunction $`\varphi _{\mathrm{in}}^\kappa (z)`$, which is an eigenstate of the incoming lead, by the relation

$$\varphi _{\mathrm{out}}^{\stackrel{~}{\kappa }}(z)=\mathrm{\Omega }^+(E)\varphi _{\mathrm{in}}^\kappa (z).$$
(9)
The wavevectors $`\kappa `$ and $`\stackrel{~}{\kappa }`$ of the two wavefunctions are uniquely identified by the relation

$$\kappa ,\stackrel{~}{\kappa }=\frac{2}{\lambda }\mathrm{arccos}\left(\frac{EE_L^{\mathrm{in},\mathrm{out}}}{2\gamma _L}\right).$$
(10)
The operator $`\mathrm{\Omega }^+(E)=𝕀+G^0(E)T(E)`$ in Eq. (9) involves the Green’s function $`G^0(E)=(EH_0)^1`$ of the two leads, as well as the scattering matrix $`T(E)=H_I[𝕀G^0(E)H_I]^1`$ where $`H_I`$ is the Hamiltonian acting on the first $`N_{\mathrm{out}}`$ sites of the lattice subject to the force $`F`$ (see C).

In the next step we renormalize the interaction Hamiltonian $`H_I`$ into an effective Hamiltonian $`\stackrel{~}{H}_I(E)`$ acting in the subspace $`\{|\mathrm{\hspace{0.17em}1},|N_{\mathrm{out}}\}`$ formed by the two edge sites that are directly connected to the leads. In this renormalization approach the scattering matrix is also projected on the subspace $`\{|\mathrm{\hspace{0.17em}1},|N_{\mathrm{out}}\}`$ and can be written as a $`2\times 2`$ matrix and evaluated with a recursive algorithm as described in C. In such a formalism the transmittivity coefficient $`𝒯(E)`$, which is defined as

$$𝒯(E)=\frac{lim_{n+\mathrm{}}n|\varphi _{\mathrm{out}}^{\stackrel{~}{\kappa }}\varphi _{\mathrm{out}}^{\stackrel{~}{\kappa }}|nv_{\mathrm{out}}}{lim_m\mathrm{}m|\varphi _{\mathrm{in}}^\kappa \varphi _{\mathrm{in}}^\kappa |mv_{\mathrm{in}}}$$
(11)
with $`v_{\mathrm{in}}`$ and $`v_{\mathrm{out}}`$ the velocities of the incoming and outgoing wavefunction, is related to the off-diagonal element $`T(E)_{1,N_{\mathrm{out}}}`$ of the scattering matrix by

$$𝒯(E)=4\frac{|T(E)_{1,N_{\mathrm{out}}}|^2}{\gamma _L^2}\mathrm{sin}\left(\kappa \frac{\lambda }{2}\right)\mathrm{sin}\left(\stackrel{~}{\kappa }\frac{\lambda }{2}\right).$$
(12)
It is easy to show that Eq. (12) is equivalent to the trace of Eq. (63) in the work of Paulsson .

The transmittivity determines the particle current passing through the lattice, and to calculate it we can exploit Landauer’s approach adapted to a bosonic gas with out-of-equilibrium leads. Again following Paulsson , the current $`j`$ can be written

$$j=\frac{N_b}{\pi \mathrm{}n_s}dE[f(E\mu _1)f(E\mu _2)]𝒯(E)$$
(13)
where $`f(E\mu )`$ is the Bose-Einstein distribution and $`\mu _{1,2}`$ are the chemical potentials of the reservoirs which act as the source and the sink of the bosons. Note that the normalization $`N_b/n_s`$ of the bosonic wavefunction plays the role of the electric charge in the Landauer formula . In the zero-temperature limit we have

$$f(E\mu _1)f(E\mu _2)\frac{f(E\mu )}{E}(\mu _2\mu _1)$$
(14)
and the density probability $`f(E\mu )/E`$ is a delta function $`\delta (E\mu )`$. Then Eq. (13) gives

$$j=\frac{N_b}{\pi \mathrm{}n_s}𝒯(\overline{\mu })\mathrm{\Delta }\mu $$
(15)
where $`\mathrm{\Delta }\mu =\mu _2\mu _1=N_{\mathrm{out}}F\lambda /2`$ is the difference in chemical potentials and $`\overline{\mu }=(\mu _1+\mu _2)/2`$ is their average value, corresponding to the dashed line in Fig. 2.

The transmittivity may also be expressed in terms of the number of bosons exiting from the lattice as a drop emitted under a constant driving force, as in the experiment of Anderson and Kasevich . The current is written in that context as

$$j=N_{\mathrm{drop}}\nu _B,$$
(16)
where $`\nu _B=(F\lambda )/(2h)`$ is the frequency of Bloch oscillations and $`N_{\mathrm{drop}}`$ is the number of bosons in the first drop. This can be calculated from Eqs.(15-16) as

$$\frac{N_{\mathrm{drop}}}{N_b}=\frac{N_{\mathrm{out}}}{n_s}𝒯(\overline{\mu }).$$
(17)
Equation (17) holds for a steady-state current and will not be applicable to the subsequent drops unless the condensate is continuously replenished.

## 4 Results

We have calculated the transmittivity of the model as a function of the acceleration $`a`$ for the four cases listed at the beginning of Sec. 3. The difference $`\mathrm{\Delta }\mu `$ of the chemical potentials of the two leads does not depend on $`a`$ and is fixed by the total band width $`4\gamma _L`$ of the leads, so that the current $`j`$ is directly proportional to the transmittivity (see Eq. (15)). We have chosen to focus on the typical experimental parameters used at LENS  for the 87Rb-40K mixture, setting $`N_b=10^5`$, $`n_s=200`$, $`U_{\mathrm{Rb}}^0=3.5Er_{\mathrm{Rb}}`$, and $`\omega _{\mathrm{Rb}}=2\pi \times 85.7`$ Hz. As previously remarked, we set $`\lambda =763`$ nm so that we can consider each 40K atom as localized in a lattice well. Moreover, for this choice of the laser wavelength, the detunings $`\delta _{\mathrm{K},\mathrm{Rb}}`$ are much larger than the natural widths $`\mathrm{\Gamma }_{\mathrm{K},\mathrm{Rb}}`$ and spontaneous emission can be neglected.

In Fig. 3 we show the transmittivity of a one-band system as a function of $`a`$ for the one-band cases in which either the fermions are absent (EF) or each well is occupied by one fermion (FF). The effect of the presence of the fermions should be that the boson-fermion attractions diminish the hopping probability of the bosons between two neighboring wells, so that the transmittivity should be lowered. However, for the case of only one fermion per well and for the usual values of the boson-fermion scattering length, this effect is not visible (see Fig. 3).

In Fig. 4 we show the transmittivity of the two-band system for $`N_f=100`$ fermions which are distributed in the lattice wells either regularly (OHF) or completely at random (DHF). The value of $`N_{\mathrm{out}}`$ in these calculations has been taken to run from 200 to 12: this corresponds to going from left to right in Fig. 4(a) and from right to left in Fig. 4(b). The transmittivity in Fig. 4 is on average strongly reduced relative to the cases shown in Fig. 3, but in the OHF case it shows a very prominent structure of peaks and troughs. The observed pattern arises from the interference occuring at $`p=h/2\lambda `$ (the left-hand edge of the Brillouin zone in the split-band configuration) between the boson wave packets that are Bragg scattered at $`p=h/2\lambda `$ in the lower sub-band and those that have tunnelled into the upper sub-band, travelled through it, and tunnelled back into the lower sub-band at $`p=h/2\lambda `$. The two times that govern the interference pattern are the Bloch oscillation period $`T_B=h/(F\lambda )`$ and the time $`\tau `$ taken by the bosons in coherently tunnelling twice across the minigap, which is proportional to $`(\mathrm{}N_b/\mathrm{\Delta }En_s)`$ . We expect constructive interference in the reflectivity at $`p=0`$ in the upper sub-band, i.e. destructive interference in the transmittivity, whenever the ratio $`T_B/\tau `$ is an integer number. This behaviour is illustrated in Fig. 4(b) by replotting the transmittivity data against $`T_B/\tau `$ with the choice $`\tau =(3\pi ^2/8)(\mathrm{}N_b/\mathrm{\Delta }En_s)`$. Of course this specific value of the proportionality constant in the tunnelling time, which is consistent with the Heisenberg principle, depends on the model Hamiltonian that we have assumed.

Returning to Fig. 4(a), the effect of disorder in the distribution of the fermions in the half-occupancy DHF case is to strongly diminish the contrast in the interference pattern. The disordered distribution that we have examined in Fig. 4 has been produced from a random number generator and the pattern is completely washed out in this case of maximum disorder. While the average occupancy of a site by a fermion is only one-half, the transmittivity of the bosons is greatly reduced relative to the FF case in Fig. 3. Evidently, the main effect at play is due to the fluctuations in the disordered distribution of the fermions, which tend to break the coherence of the Bose condensate.

## 5 Concluding remarks

In conclusion, we have studied the condensate current induced by a constant external force in a 1D optical lattice in the presence of various numbers of fermionic atoms, distributed over the lattice in orderly or random ways. We have focused on the 87Rb-40K mixture, with a choice of the laser generating the lattice potential such that the barriers for the fermions are much higher than those for the Bose condensate. In this regime fermionic transport can be neglected, but the presence of the fermions and their spatial distribution affect the bosonic transport via mean-field attractive interactions.

We have calculated the transmittivity and the condensate current without explicitly diagonalizing the Bose-Hubbard Hamiltonian, but by exploiting an out-of-equilibrium Green’s function approach adapted to the bosonic case. When there is one fermion per well, the current is essentially the same as in the absence of fermions and increases monotonically with the driving force. When every second well is occupied by a fermion, instead, the energy spectrum breaks into two sub-bands and the probability that the condensate may tunnel through the band minigap depends on the energy that it has acquired at the edge of the first sub-band and on the energy width of the minigap. Part of the condensate is back-scattered when it arrives at the edge of the lower sub-band and can interfere at the opposite edge of the Brillouin zone with the part that has tunnelled to the upper sub-band and travelled through it. This gives rise to an interference pattern showing marked maxima and minima in the transmittivity, whose location is an integer multiple of a parameter directly proportional to the energy minigap and inversely proportional to the driving force.

The interference pattern loses resolution when the distribution of the fermions in the lattice becomes disordered , up to being completely washed out in the case of full randomness. An ordered period doubling of the lattice could be realized, however, with a four-beam experimental set-up. Observation of resonances in the current as a function of the driving force would allow a measurement of the tunnelling time through the minigap.

This work was supported by INFM through the PRA2001-Photonmatter. ZA also acknowledges support from TUBITAK and from the Research Fund of the University of Istanbul under Project Number BYP-110/12122002.

## Appendix A Determination of radial and axial widths

We adopt the variational method proposed for Bose condensates by Salasnich et al. . Within a 1D nonpolynomial nonlinear Schrödinger equation approach, we find that the evaluation of the transverse widths $`\sigma _{\mathrm{Rb},\mathrm{K}}`$ requires the self-consistent solution of the Euler-Lagrange equations

$$\{\begin{array}{c}\frac{\mathrm{}^2}{2m_{\mathrm{Rb}}\sigma _{\mathrm{Rb}}^2}\frac{1}{2}m_{\mathrm{Rb}}\omega _{\mathrm{Rb}}^2\sigma _{\mathrm{Rb}}^2+\frac{1}{2}g_{bb}|\varphi _i(z)|^2+g_{bf}\frac{\sigma _{\mathrm{Rb}}^2}{\sigma _{\mathrm{Rb}}^2+\sigma _\mathrm{K}^2}n_f(z)=0\hfill \\ \frac{\mathrm{}^2}{2m_\mathrm{K}\sigma _\mathrm{K}^2}\frac{1}{2}m_\mathrm{K}\omega _\mathrm{K}^2\sigma _\mathrm{K}^2+g_{bf}\frac{\sigma _\mathrm{K}^2}{\sigma _{\mathrm{Rb}}^2+\sigma _\mathrm{K}^2}|\varphi _i(z)|^2=0.\hfill \end{array}$$
(18)
Assuming that $`\sigma _{\mathrm{Rb},\mathrm{K}}`$ are essentially constant along the axial position, we set in Eqs.(18) $`|\varphi _i(z)|^2N_b/(\sqrt{\pi }n_s\sigma _{z\mathrm{Rb}})`$ and $`n_f(z)1/(\sqrt{\pi }\sigma _{z\mathrm{K}})`$ or $`n_f(z)=0`$, depending on the presence or absence of a fermion in the corresponding well. The axial widths $`\sigma _{z\mathrm{Rb},\mathrm{K}}`$ entering the axial densities depend themselves on the radial widths $`\sigma _{\mathrm{Rb},\mathrm{K}}`$. Consequently we first solve Eqs.(18) for given values of $`\sigma _{z\mathrm{Rb},\mathrm{K}}`$, evaluated in the absence of interactions and then, having calculated the radial widths, we evaluate the correction to the axial widths.

Minimization of the bosonic and fermionic energies with respect to $`\sigma _{z\mathrm{Rb},\mathrm{K}}`$ yields the two coupled equations

$$\{\begin{array}{c}\frac{\mathrm{}^2}{2m_{\mathrm{Rb}}\sigma _{z\mathrm{Rb}}^2}\frac{1}{2}m_{\mathrm{Rb}}\omega _{z\mathrm{Rb}}^2\sigma _{z\mathrm{Rb}}^2+\frac{\mathrm{}^2}{2m_{\mathrm{Rb}}\sigma _{z\mathrm{Rb}}^2}\frac{g_{bb}N_b}{2\sqrt{2\pi }n_s\sigma _{z\mathrm{Rb}}}+\frac{g_{bf}}{\sqrt{\pi }}\frac{\sigma _{z\mathrm{Rb}}^2}{\left(\sigma _{z\mathrm{Rb}}^2+\sigma _{z\mathrm{K}}^2\right)^{3/2}}=0\hfill \\ \frac{\mathrm{}^2}{2m_\mathrm{K}\sigma _{z\mathrm{K}}^2}\frac{1}{}2m_\mathrm{K}\omega _{z\mathrm{K}}^2\sigma _{z\mathrm{K}}^2+\frac{g_{bf}N_b}{\sqrt{\pi }n_s}\frac{\sigma _{z\mathrm{K}}^2}{\left(\sigma _{z\mathrm{Rb}}^2+\sigma _{z\mathrm{K}}^2\right)^{3/2}}=0\hfill \end{array}$$
(19)
for the case of a well containing a fermion, while in the absence of a fermion we have

$$\frac{\mathrm{}^2}{2m_{\mathrm{Rb}}\sigma _{z\mathrm{Rb}}^2}\frac{1}{2}m_{\mathrm{Rb}}\omega _{z\mathrm{Rb}}^2\sigma _{z\mathrm{Rb}}^2+\frac{g_{bb}N_b}{2\sqrt{2\pi }n_s\sigma _{z\mathrm{Rb}}}=0$$
(20)
as already found by Chiofalo et al. . Here, the axial frequencies $`\omega _{z\mathrm{Rb},\mathrm{K}}`$ for a single well are expressed in terms of the potential depths $`U_{\mathrm{Rb},\mathrm{K}}^0`$ and of the recoil energies $`Er_{\mathrm{Rb},\mathrm{K}}=\mathrm{}^2k^2/2m_{\mathrm{Rb},\mathrm{K}}`$ as

$$\omega _{z\mathrm{Rb},\mathrm{K}}=\frac{1}{\mathrm{}}\sqrt{4U_{\mathrm{Rb},\mathrm{K}}^0Er_{\mathrm{Rb},\mathrm{K}}}.$$
(21)
The solution of Eqs. (18)-(20) provides the further input that is needed in Eqs. (6) and (8).

## Appendix B Period doubling with a four-beam set-up

Four laser beams can be used to create a 1D lattice with a double periodicity along the $`z`$ direction. Two beams are set in a counterpropagating configuration along the $`z`$ axis, while the other two beams are rotated by angles of 60 and 120 with respect to the $`z`$ axis. In the region of space where the four beams are superposed, the longitudinal electric field $`E`$ is proportional to

$$E\beta e^{i\alpha }(e^{ikz/2}+e^{i\varphi _1}e^{ikz/2})+(e^{ikz}+e^{i\varphi _2}e^{ikz}).$$
(22)
In Eq. (22) $`\beta `$ and $`\alpha `$ are the relative amplitude and phase of the rotated beams with respect to the counterpropagating pair, $`\varphi _1`$ is the relative phase of the rotated beams, and $`\varphi _2`$ the relative phase of the counterpropagating beams. All these parameters can be set by choosing suitable optical components.

Upon setting $`\varphi _1=\varphi _2=\pi `$ and $`\alpha =\pi /2`$, the resulting lattice potential will have the shape

$$U(z)\beta ^2\mathrm{sin}^2(kz/2)+\mathrm{sin}^2(kz)$$
(23)
and the double period enters through the well depths, that is the site energies, which will differ by a factor $`\beta ^2`$. Instead, by setting $`\varphi _1=\varphi _2=0`$ and $`\alpha =\pi /2`$ we have

$$U(z)\beta ^2\mathrm{cos}^2(kz/2)+\mathrm{cos}^2(kz)$$
(24)
and the potential hills (that is the hopping energies) have alternating values. In the case of a 87Rb condensate, a lattice potential of the type in Eq. (23) with $`\beta 0.1`$ is equivalent to the effective potential generated by a regular distribution of 40K atoms inside every other well of an optical lattice created with a single pair of laser beams.

## Appendix C The recursive algorithm

Figure 5 shows a schematic representation of the whole Hamiltonian for the system composed of a finite lattice and two leads, and of its decomposition into the sum of $`H_0`$ and $`H_I`$ for the evaluation of the scattering matrix. We show below how $`H_I`$ can be reduced to an effective Hamiltonian $`\stackrel{~}{H}_I`$ for a dimer in a subspace $`\{|\mathrm{\hspace{0.17em}1},|N_{\mathrm{out}}\}`$ and how this leads to an expressions for the scattering matrix and for the outgoing wavefunctions.

Let $`n`$ be the site index number. To calculate the transmitted wavefunction we have to evaluate

$$n|\varphi _{\mathrm{out}}^{\stackrel{~}{\kappa }}=\underset{m}{}n|𝕀+G^0(E)T(E)|mm|\varphi _{\mathrm{in}}^\kappa $$
(25)
for $`n>N_{\mathrm{out}}`$. Since the wavefunctions $`\varphi _{\mathrm{in}}^\kappa `$ and $`\varphi _{\mathrm{out}}^{\stackrel{~}{\kappa }}`$ are defined in two disconnected spaces, Eq. (25) becomes

$`n|\varphi _{\mathrm{out}}^{\stackrel{~}{\kappa }}`$ $`=`$ $`{\displaystyle \underset{m1,l}{}}n|G^0(E)|ll|T(E)|mm|\varphi _{\mathrm{in}}^\kappa `$ (26)
$`=`$ $`n|G^0(E)|N_{\mathrm{out}}N_{\mathrm{out}}|T(E)|\mathrm{\hspace{0.17em}1}1|\varphi _{\mathrm{in}}^\kappa `$
$`=`$ $`\sqrt{2}[G^0(E)]_{n,N_{\mathrm{out}}}[T(E)]_{N_{\mathrm{out}},1}\mathrm{sin}(\kappa \lambda /2)u_\kappa (\lambda /2),`$
where we have set $`1|\varphi _{\mathrm{in}}^\kappa =u_\kappa (\lambda /2)[\mathrm{exp}(i\kappa \lambda /2)\mathrm{exp}(i\kappa \lambda /2)]/i\sqrt{2}`$ with $`u_\kappa (\lambda /2)=\varphi _i(z)`$. It can be shown that Eq. (26) is equivalent to Eq. (34) in the work of Paulsson  by making use of the relation $`G^0(E)T(E)=G(E)H_I(E)`$.

The Green’s function element $`[G^0(E)]_{n,N_{\mathrm{out}}}`$ in Eq.(26) determines the coherence between site $`n`$ and site $`N_{\mathrm{out}}`$ on the chain for the outgoing lead. It can be written as

$$[G^0(E)]_{n,N_{\mathrm{out}}}=\frac{\gamma _L^{nN_{\mathrm{out}}}}{|\gamma _L^{n+1N_{\mathrm{out}}}|}e^{i\stackrel{~}{\kappa }(n+1N_{\mathrm{out}})\lambda /2}$$
(27)
in terms of quantities defined in the main text. The scattering matrix element $`[T(E)]_{N_{\mathrm{out}},1}`$ can be evaluated from the effective Hamiltonian $`\stackrel{~}{H}_I(E)`$ as schematically illustrated in Fig. 5:

$$[T(E)]_{N_{\mathrm{out}},1}=N_{\mathrm{out}}|\stackrel{~}{H}_I(E)[𝕀_2\stackrel{~}{G}^0(E)\stackrel{~}{H}_I(E)]^1|\mathrm{\hspace{0.17em}1}$$
(28)
where $`𝕀_2`$ is the $`2\times 2`$ identity matrix and $`\stackrel{~}{G}^0(E)`$ is the Green’s function of the leads projected onto the subspace $`\{|\mathrm{\hspace{0.17em}1},|N_{\mathrm{out}}\}`$. The only non-zero elements of the $`2\times 2`$ matrix $`\stackrel{~}{G}^0(E)`$ are the diagonal ones, which are given by $`𝒢(EE_L^{\mathrm{in}},\gamma _L)`$ and $`𝒢(EE_L^{\mathrm{out}},\gamma _L)`$, with

$$𝒢(EE_s,\gamma _L)=\frac{1}{2\gamma _L^2}[EE_s\sqrt{(EE_s)^24\gamma _L^2}].$$
(29)
Finally, the elements of the matrix $`\stackrel{~}{H}_I(E)`$ in Eq. (28) can be calculated by decimation of the sites (2, 3, $`\mathrm{}`$, $`N_{\mathrm{out}}1`$) and can be written in the form of continued fractions. That is,

$$[\stackrel{~}{H}_I(E)]_{1,1}=[H_I]_{1,1}+\frac{[H_I]_{1,2}^2}{E[H_I]_{2,2}{\displaystyle \frac{[H_I]_{2,3}^2}{E[H_I]_{3,3}\mathrm{}}}}$$
(30)
and

$$[\stackrel{~}{H}_I(E)]_{1,N_{\mathrm{out}}}=[H_I]_{1,2}\frac{1}{E[H_I]_{2,2}}[H_I]_{2,3}\frac{1}{E[H_I]_{3,3}{\displaystyle \frac{[H_I]_{3,2}^2}{E[H_I]_{2,2}}}}[H_I]_{3,4}\mathrm{},$$
(31)
with $`[\stackrel{~}{H}_I(E)]_{N_{\mathrm{out}},1}=[\stackrel{~}{H}_I(E)]_{1,N_{\mathrm{out}}}`$ and $`[\stackrel{~}{H}_I(E)]_{N_{\mathrm{out}},N_{\mathrm{out}}}`$ given similarly to $`[\stackrel{~}{H}_I(E)]_{1,1}`$ in Eq. (30).

## References


format: text

text:
# SIS epidemics with household structure: the self-consistent field method

## 1 Introduction

In this paper we will give a treatment of an SIS infection model in a population partioned into households. Our main emphasis will be the use of a method borrowed from statistical physics, the self-consistent field (SCF) technique.

The susceptible-infected-susceptible (SIS) model is one of the basic themes in mathematical epidemiology . In this idealization of the spread and persistance of an infection, $`N`$ individuals are initially either susceptible, $`S`$, or infected, $`I`$. Each one of the infected is coupled to a certain number of other agents. Each $`I`$ infects all of the $`S`$ with which it has contact with a rate $`B`$, and each $`I`$ recovers with rate $`\gamma `$, and immediately becomes susceptible again. Schematically, we can write:

$$S+I2IIS$$

Clearly, at long times the system might achieve an equilibrium where the number of infected fluctuates around a mean. We characterize this state by denoting the mean fraction of population which is infected by $`f`$, i.e.,

$$f=_t/N,$$
(1)
where $`<>_t`$ denotes a time average. We refer to this quantity as the endemic level.

The model is defined by the two parameters, $`B,\gamma `$, and by the contact structure of the population. Two cases of the latter have been studied in some detail. In the first, each individual is equally coupled to all the others in a large population. This random mixing model is rather easy to treat, and a large number of results are known . For example, in the deterministic limit, $`I`$ obeys:

$`dI/dt`$ $`=`$ $`bIS/N\gamma I`$
$`S`$ $`=`$ $`NI`$ (2)
where we have redefined $`B=b/N`$. It is easy to show that for $`b>b_c\gamma `$, and for any intial condition with $`I>0`$, the epidemic goes to the endemic level, $`f=1\gamma /b`$. For $`b<b_c`$ the infection dies out.

In the second well-studied case the individuals are arranged in some geometric way and couple only to their nearest neighbors. This is called the contact process . In the last section of this paper we give numerical results on a model which interpolates between the two cases.

Ball  introduced the SIS household model, which puts more structure into the random mixing case by having two levels of mixing. In this idealization the population, $`N`$, is partitioned into $`m`$ households with $`n`$ members. The infection rate within the household, $`w`$, is different from that with individuals outside, $`B`$. This followed work on SIR models with similar structure . Ball’s results are for the case $`m\mathrm{}`$ with fixed $`n`$, that is, a large population of small households, which is clearly of practical interest.

In this paper we extend Ball’s work from a different point of view inspired by techniques in statistical physics. These will be explained in the next section. Briefly, the idea is that for random mixing and a large population, we can consider that each household is acted on by an average force of infection from outside. This is quite similar to methods used in the theory of magnetism , inspired by the work of P. Weiss. Weiss introduced the ‘molecular field’, which is proportional to the mean orientation of the neighbors of a given magnetic ion. The orientation of the ion itself is then determined, and the unknown molecular field is gotten by demanding that the neighbors have the same orientation. (Note that this step assumes that there is no correlation between neighbors.) This step is known as self-consistency, hence the name, self-consistent field. As we will see, $`f`$ plays the role of the self-consistent field. The technique itself is called self-consistent field (or sometimes mean-field) theory.

Some of our results were already obtained in , e.g., the expression for $`b_c`$ for households of size $`n`$, Eq. (26), below. Other explicit results are new. We think that one of the main values of this paper is that this point of view allows us to look at the solutions in a new way, and gives considerable insight.

In the last section we continue the analogy to magnetic systems by doing numerical simulations for finite $`m`$, and for systems where the random-mixing assumption is not valid. We find, as in magnetism, that the SCF assumption is remarkably accurate in situations where we would expect it to fail.

## 2 Stochastic formulation and steady-state solutions

Consider an SIS epidemic in a population partitioned into $`m`$ households each consisting of $`n`$ individuals. We suppose that within a household the infection rate per contact is $`w`$ and between different households $`B`$. We write $`B=b/[(m1)n]`$. (In the notation of Ball , $`\lambda _w/n=w,\lambda _b=b`$ ). We suppose that $`b`$ remains finite in the limit $`m\mathrm{}`$. The recovery rate is $`\gamma `$, and a recovered individual can be reinfected at once.

### 2.1 Self-consistent field equations

We begin intuitively by considering the case $`n=2`$ and try to formulate the Chapman-Kolmogorov forward equations for the probability, $`P_k`$, to have $`k`$ infected in a typical household:

$`{\displaystyle \frac{dP_0}{dt}}`$ $`=`$ $`\gamma P_12bfP_0`$
$`{\displaystyle \frac{dP_1}{dt}}`$ $`=`$ $`2bfP_0^j(\gamma +w+bf)P_1+2\gamma P_2^j`$
$`{\displaystyle \frac{dP_2}{dt}}`$ $`=`$ $`(w+bf)P_1^j2\gamma P_2`$ (3)
$`f`$ $`=`$ $`(P_1+2P_2)/2`$ (4)
The terms in $`b`$ are the mean force of infection between households, and $`f`$ has the interpretation of the fraction of a typical household that is infected. The point is that the only communication between households is by infection, and for a large number of such contacts, only the average number of infected in the other households is important. In effect, we have assumed that the households are independent. Note that $`_kP_k=1`$, and $`0f1`$.

The strategy for solving these equations (in the steady state) is to find the $`P_k`$ with $`f`$ as a parameter, and solve for $`f`$ self-consistently, using Eq. (4), at the end of the calculation. This is a heuristically appealing method. For magnetic models, the approximation is exact in the case of infinite-range interactions (in our case, all households interact equally with all others). Even for short-range interactions (i.e., households interacting with a ‘neighborhood’) the approximation is often surprisingly useful.

In fact, a rigorous proof of the validity of Eqs. (3, 4) is available for this case as well. This is given by Ball  who uses the work of Kurtz  to prove the following: let $`X_k^{(m)}`$ be the number of households in the population with exactly $`k`$ infected. Define $`P_k=lim_m\mathrm{}X_k^{(m)}/m`$. Then $`P_k`$ satisfies Eqs. (3, 4) for any initial condition with finite density, i.e., with $`P_k(0),k>0`$ not all equal to $`0`$. In the last section we will investigate numerically how good this approximation is for finite $`m`$ and for finite range of interaction.

We note, for later use, that we can give another interpretation to Eq. (3). Suppose we consider a single household with 2 individuals subject to a *fixed* external force of infection. Then we can think of $`P_k`$ as the probability to have $`k`$ infected. From standard results on stochastic processes, there is a unique steady-state solution to these equations for any positive $`f,b,w`$. Further, if $`f=0`$ the only steady-state solution is $`P_0=1,P_1=P_2=0`$; that is, the infection becomes extinct with probability unity.

In the general case the time evolution of the vector of probabilities $`𝐏=(P_0,P_1,\mathrm{},P_n)^T`$ (T denotes a column vector) to have exactly $`k`$ infected in the household of $`n`$, is given by:

$$\frac{d𝐏}{dt}=𝐌(f)𝐏$$
(5)
where $`𝐌(f)`$ is a square tridiagonal matrix with nonzero elements

$`m_{k,k1}`$ $`=`$ $`\left(nk+1\right)\left[\left(k1\right)w+fb\right]`$ (6)
$`m_{k,k}`$ $`=`$ $`\left[k\gamma +(nk)kw+(nk)fb\right]`$
$`m_{k,k+1}`$ $`=`$ $`\left(k+1\right)\gamma `$
($`0kn`$). The equations again contain $`f`$, the mean level of infection, as a parameter.

$$f=\frac{1}{n}\underset{k=0}{\overset{n}{}}kP_k$$
(7)
For the reasons given above, these equations are valid in the limit $`m\mathrm{}`$.

We will be interested in the stationary solutions of Eqs. (5, 7). It is clear by inspection that there is always a trivial steady-state solution $`f=0,P_0=1,P_k=0,k>0`$ corresponding to extinction of the infection. We will find that above a certain threshold there is another solution, and we will show that when this occurs the trivial solution is unstable.

### 2.2 Steady-state solutions

In the case $`n=2`$ it is a matter of simple algebra to find the steady-state solution:

$`P_1`$ $`=`$ $`{\displaystyle \frac{2fb\gamma }{\gamma ^2+2fb\gamma +(fb)^2+fbw}}`$
$`P_2`$ $`=`$ $`{\displaystyle \frac{(fb)^2+fbw}{\gamma ^2+2fb\gamma +(fb)^2+fbw}}`$ (8)
Using Eq. (4) we get a cubic polynomial equation to solve for the endemic level, $`f`$. After removing the trivial solution we have:

$$(fb)^2+f(2b\gamma +wbb^2)+(\gamma ^2b\gamma wb)=0$$
(9)
The larger root of this equation is:

$$f=\frac{(bw2\gamma )+\sqrt{(b+w)^2+4w\gamma }}{2b}$$
(10)
This expression is acceptable, i.e., positive, provided:

$$b(\gamma +w)/\gamma ^2>1$$
(11)
Thus Eq. (11) serves as a threshold condition for the existence of a non-trivial endemic level.

For $`n=3`$ it is not difficult to find a cubic equation for $`f`$ analogous to Eq. (9) (in this expression we set $`\gamma =1`$, which sets the unit of time):

$`f^3b^3f^2(b^33wb^23b^2)`$ $``$ $`f(b(2b+3bw)2bw^23b3bw)`$
$`b(1+2w+2w^2)+1`$ $`=`$ $`g(f)=0`$ (12)
By inspection $`g(1)>0`$. If $`g(0)=b(1+2w+2w^2)+1<0`$ then there is at least one root in the interval $`[0,1]`$. (In fact, numerics shows that there are also two negative roots). Thus if:

$$b(1+2w+2w^2)>1$$
(13)
we have a non-trivial solution, and otherwise all the roots of Eq. (12) are all negative, and only the trivial solution is acceptable.

In the general case we look for the stationary solution of Eq. (5) by setting $`d𝐏/dt=0`$. We get a homogeneous system of equations, from which the endemic level has to be established self-consistently. Recall that $`_kP_k=1`$, so that we can reduce the number of equations by one by eliminating $`P_0`$. Denote the vector of length $`n`$, $`𝐩=(P_1,\mathrm{},P_n)^T`$. The resulting system of equations for the stationary state is inhomogeneous. It reads

$`𝐒(f)𝐩`$ $`=`$ $`nfb𝐮`$
$`𝐮`$ $`=`$ $`(1,0,0,\mathrm{})^T`$ (14)
The matrix $`𝐒`$ is given by:

$$\left(\begin{array}{cccccc}m_{1,1}nfb\hfill & m_{1,2}nfb\hfill & nfb\hfill & nfb\hfill & \mathrm{}\hfill & nfb\hfill \\ m_{21}\hfill & m_{22}\hfill & m_{23}\hfill & 0\hfill & \mathrm{}\hfill & 0\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ 0\hfill & \mathrm{}\hfill & m_{k,k1}\hfill & m_{k,k}\hfill & m_{k,k+1}\hfill & 0\mathrm{}\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ 0\hfill & 0\hfill & \mathrm{}\hfill & 0\hfill & m_{n,n1}\hfill & m_{nn}\hfill \end{array}\right)$$
(15)
The matrix $`𝐒(f)`$ is linear in $`f`$:

$$𝐒(f)=𝐖+fb𝐁$$
(16)
where both matrices, $`𝐖`$ and $`𝐁`$ are $`f`$-independent. The structure of $`𝐖`$ is:

$$\left(\begin{array}{cccc}\gamma (n1)w\hfill & 2\gamma \hfill & 0\hfill & \mathrm{}\hfill \\ (n1)w\hfill & 2\gamma 2(n2)w\hfill & 3\gamma \hfill & \mathrm{}\hfill \\ 0\hfill & 2(n2)w\hfill & 3\gamma 3(n3)w\hfill & 4\gamma \hfill \\ 0\hfill & 0\hfill & 3(n3)w\hfill & \mathrm{}\hfill \\ 0\hfill & \mathrm{}\hfill & (n1)w\hfill & n\gamma \hfill \end{array}\right)$$
(17)
$`𝐖`$ depends only on the parameters $`\gamma `$ and $`w`$, characterizing the infection within a household. Note that it is nondegenerate: $`det𝐖=(\gamma )^nn!`$. To see this, add to each row the sum of all the rows below. The result is a lower triangular matrix with diagonal elements in the $`n`$-th row equal to $`n\gamma `$. Moreover, $`𝐁`$ is a matrix with integer components:

$$𝐁=\left(\begin{array}{ccccc}(2n1)\hfill & n\hfill & n\hfill & n\hfill & \mathrm{}\hfill \\ (n1)\hfill & (n2)\hfill & 0\hfill & 0\hfill & \mathrm{}\hfill \\ 0\hfill & (n2)\hfill & (n3)\hfill & 0\hfill & \mathrm{}\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ 0\hfill & 0\hfill & \mathrm{}\hfill & 1\hfill & 0\hfill \end{array}\right).$$
(18)
We give the explicit forms of the matrices $`𝐖`$ and $`𝐁`$ for $`n=2`$ and $`n=3`$, setting $`\gamma =1`$: For $`n=2`$ one has

$$𝐖=\left(\begin{array}{cc}1w\hfill & 2\hfill \\ w\hfill & 2\hfill \end{array}\right)\text{ and }𝐁=\left(\begin{array}{cc}3\hfill & 2\hfill \\ 1\hfill & 0\hfill \end{array}\right),$$
(19)
For $`n=3`$ we get

$$𝐖=\left(\begin{array}{ccc}12w\hfill & 2\hfill & 0\hfill \\ 2w\hfill & 22w\hfill & 3\hfill \\ 0\hfill & 2w\hfill & 3\hfill \end{array}\right)\text{ and }𝐁=\left(\begin{array}{ccc}5\hfill & 3\hfill & 3\hfill \\ 2\hfill & 1\hfill & 0\hfill \\ 0\hfill & 1\hfill & 0\hfill \end{array}\right).$$
(20)
From the remarks above, it is clear that there is a non-zero solution, $`𝐩`$, to Eq. (14) for any $`f>0`$, so that $`𝐒^1`$ exists. Thus:

$$𝐩=nfb𝐒^1𝐮$$
(21)
From this equation a closed form for the endemic level $`f`$ follows. Multiplying both sides of Eq.(21) by the vector $`𝐧_\mathrm{𝟏}=(1,2,\mathrm{},n)`$ and noting that $`f=n^1𝐧_\mathrm{𝟏}𝐩`$ we get:

$$f=fb𝐧_\mathrm{𝟏}𝐒^1𝐮$$
(22)
which is a closed algebraic equation for $`f`$ which always has the trivial solution $`f=0`$. Factoring out this out, we get an equation for the nontrivial solution:

$$b𝐧_\mathrm{𝟏}𝐒^1𝐮=1.$$
(23)
As in Eqs. (11, 13) there is a condition on the parameters for the solution to be positive. For given $`w`$ we define a value $`b_c`$ such that for $`b>b_c`$ we have such an $`f`$. This can be found from Eq.(23) by noting that at $`b_c`$ the second solution crosses $`f=0`$. At that point, $`𝐒=𝐖`$, and we get:

$$b_c=1/\left[𝐧_\mathrm{𝟏}𝐖^1𝐮\right].$$
(24)
To get an explicit expression for $`b_c`$ consider the vector $`𝐯=𝐖^1𝐮`$, whose elements are $`𝐯_k=\left[𝐖^1\right]_{1,k}`$. As we will show,

$$𝐯_k=\frac{(n1)!}{k(nk)!}w^{k1},$$
(25)
so that the threshold is given by

$$b_c^1=𝐧_\mathrm{𝟏}𝐯=\underset{k=1}{\overset{n}{}}\frac{(n1)!}{(nk)!}w^{k1}.$$
(26)
To prove Eq.(25) we note that the solution to $`\mathrm{𝐖𝐯}=𝐮`$ is unique so that we need only insert Eq. (25) into $`\mathrm{𝐖𝐯}=𝐮`$ . The first element is

$$(1(n1)w)(n1)w=1,$$
(27)
while the $`k`$-th element ($`k>1`$) is:

$`(k1)(n+1k)w{\displaystyle \frac{(n1)!}{(n+1k)!(k1)}}w^{k2}`$
$`+(kk(nk)w){\displaystyle \frac{(n1)!}{(nk)!k}}w^{k1}+`$
$`+(k+1){\displaystyle \frac{(n1)!}{(nk1)!(k+1)}}w^k=`$
$`=`$ $`(n1)!\{[{\displaystyle \frac{1}{(nk)!}}{\displaystyle \frac{1}{(nk)!}}]w^{k1}+`$
$`+[{\displaystyle \frac{1}{(nk1)!}}{\displaystyle \frac{1}{(nk1)!}}]w^k\}`$
$`=`$ $`0.`$
Once we obtain $`𝐩`$ we can find other moments of the distribution by multiplying Eq.(21) by $`𝐧_𝐪=(1,2^q,3^q,\mathrm{},n^q)`$. We present numerical results on $`𝐧_\mathrm{𝟐}𝐩`$ below.

We now discuss the behavior of the solution in the vicinity of the transition point. Eq.(21) can be rewritten:

$$𝐩=nfb\left(𝐈+fb𝐖^1𝐁\right)^1𝐖^1𝐮,$$
(29)
Then we formally expand in a geometric series:

$$𝐩=nfb\left(𝐈fb𝐖^1𝐁+(fb)^2\left(𝐖^1𝐁\right)^2\mathrm{}\right)𝐖^1𝐮.$$
(30)
Multiplying both sides of Eq.(30) by $`n^1𝐧_\mathrm{𝟏}`$ we get the following equation for $`f`$:

$$f=fb\left(a_0fba_1+(fb)^2a_2+\mathrm{}+(1)^m(fb)^ma_m+\mathrm{}\right)$$
(31)
where

$$a_m=𝐧_\mathrm{𝟏}\left(𝐖^1𝐁\right)^m𝐖^1𝐮$$
(32)
from which $`f`$ can be obtained with any necessary accuracy. Note that according to Eq.(24) the coefficient

$$a_0=𝐧_\mathrm{𝟏}𝐖^1𝐮=1/b_c.$$
(33)
Eq.(31) is especially useful in the vicinity of transition, where $`f`$ is small. Both the point of transition and the endemic level near the threshold can be determined.

Eq.(31) always possesses a trivial solution $`f=0`$. An equation for a nontrivial one then reads:

$$G(f)=fb^2a_1f^2b^3a_2+\mathrm{}+(1)^{m+1}f^mb^{m+1}a_m+\mathrm{}=1+ba_0.$$
(34)
The formal solution of this equation involves the inverse of $`G`$:

$$f=G^1(1+ba_0)=G^1(1b/b_c).$$
(35)
where Eq (33) has been used. Using the relationship between the Taylor series of a function and its inverse we get, in the vicinity of $`b_c`$:

$$f=\beta _1\tau +\beta _2\tau ^2\mathrm{}$$
(36)
where $`\tau =(bb_c)/b_c`$ and the coefficients are:

$`\beta _1`$ $`=`$ $`{\displaystyle \frac{1}{b^2a_1(w)}},`$
$`\beta _2`$ $`=`$ $`{\displaystyle \frac{a_2(w)}{b^3a_1^3(w)}},`$ (37)
$`\beta _3`$ $`=`$ $`{\displaystyle \frac{1}{b^4a_1^5(w)}}\left(2a_2(w)a_1(w)a_3(w)\right),`$
$`\mathrm{}`$
In the vicinity of the transition one can take $`b=b_c(w)`$; the behavior of $`f`$ is dominated by the linear term:

$$f=\beta _1\tau =\frac{1}{b^2a_1(w)}\tau $$
(38)
Explicit calculations show that $`a_1<0`$.

### 2.3 Stability of the disease-free equilibrium

Restoring the time-dependence in Eq.(14) we obtain an equation for the time evolution of the reduced probability vector $`𝐩`$:

$$\frac{d}{dt}𝐩=(𝐖+fb𝐁)𝐩+nfb𝐮.$$
(39)
Since $`𝐩`$ depends on $`f`$, Eq.(39) is a nonlinear equation describing the evolution of $`𝐩`$. In the previous section we have shown that this equation has two equilibria provided $`b>b_c`$, where $`b_c`$ is defined in Eq. (26). We will now show that under the same conditions the trivial solution, $`f=0,𝐩=\mathrm{𝟎}`$ is unstable.

In the vicinity of the trivial equilibrium, we can omit the term $`fb𝐁`$ in brackets in Eq. (39). The linearized equation reads

$$\frac{d}{dt}𝐩=\mathrm{𝐖𝐩}+b(𝐧_\mathrm{𝟏}𝐩)𝐮\mathrm{𝐓𝐩}$$
(40)
where

$$𝐓=𝐖+b𝐍,$$
(41)
The matrix $`𝐍`$ has $`𝐧_\mathrm{𝟏}`$ as its first row and all other elements zero.

The determinant of $`𝐓`$ is linear in $`b`$. To see this, note that the explicit form of $`𝐓`$ is:

$$\left(\begin{array}{cccc}\gamma (n1)w+b\hfill & 2\gamma +2b\hfill & 3b\hfill & \mathrm{}\hfill \\ (n1)w\hfill & 2\gamma 2(n2)w\hfill & 3\gamma \hfill & \mathrm{}\hfill \\ 0\hfill & 2(n2)w\hfill & 3\gamma 3(n3)w\hfill & \mathrm{}\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ 0\hfill & \mathrm{}\hfill & (n1)w\hfill & n\gamma \hfill \end{array}\right)$$
(42)
If we expand the determinant about the first row, each term contains $`b`$ only once. It can thus be written $`det𝐓=det𝐖+cb`$, with $`c`$ a constant. However, the determinant of $`𝐓`$ is the product of its eigenvalues. Since for $`b=0`$ (separate households) the trivial solution is stable, all of the eigenvalues of $`𝐖`$ are negative; $`det𝐖`$ is nonzero, and its sign is positive for $`n`$ even and negative for $`n`$ odd. From the linearity in $`b`$ we see that $`det𝐓`$ changes sign at $`b=c^1det𝐖`$ which corresponds to the appearance of a positive eigenvalue.

To find the $`b`$ for which $`det𝐓=0`$ write Eq.(41) as:

$$\mathrm{𝐓𝐖}^1=𝐈+b\mathrm{𝐍𝐖}^1.$$
(43)
The determinant of the matrix in the l.h.s. of this expression is

$$det𝐓/det𝐖=\left[(\gamma )^n/n!\right]det𝐓$$

and vanishes when $`det𝐓=0`$. But, since $`𝐍`$ has non-zero elements only in the first row, the matrix on the r.h.s. of Eq. (43) is:

$$\left(\begin{array}{ccccc}1+bq_1\hfill & bq_2\hfill & bq_3\hfill & \mathrm{}\hfill & bq_n\hfill \\ 0\hfill & 1\hfill & 0\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ 0\hfill & 0\hfill & 1\hfill & 0\hfill & \mathrm{}\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ 0\hfill & 0\hfill & \mathrm{}\hfill & 0\hfill & 1\hfill \end{array}\right),$$
(44)
where $`q_i`$ is the $`i`$-th element of the vector $`𝐧_\mathrm{𝟏}𝐖^1`$. By inspection, the determinant of this matrix is $`1+q_11+b𝐧_\mathrm{𝟏}𝐖^1𝐮`$. It vanishes at $`b=1/(𝐧_\mathrm{𝟏}𝐖^1𝐮)`$, i.e. exactly at $`b=b_c`$, see Eq.(24). Therefore, for $`b>b_c`$ a positive eigenvalue appears, and the trivial solution loses its linear stability.

### 2.4 Numerical results

It is quite simple to work out the dependence of $`f`$ on the parameters $`w,b,n`$. We wrote a Matlab program which solved Eq. (23). In Figure (1) we show $`f`$ as a function of $`b`$ for several $`w`$’s for $`n=3`$. In this and the following figures we have taken $`\gamma =1`$. For other values of $`\gamma `$ one can use the obvious scaling relation $`(w,b,\gamma )(w/\gamma ,b/\gamma ,1)`$.

Recall that $`f=𝐧_\mathrm{𝟏}𝐩`$ can be interpretated as $`<i>/n`$ where $`i`$ is the number infected in a household, and $`<>`$ is the average over the distribution $`𝐏`$. If we calculate $`𝐧_\mathrm{𝟐}𝐩`$ we get a measure of the fluctuation of $`i`$ in a household. In Figure (2) we show the standard deviation, $`\sigma =[<i^2><i>^2]^{1/2}=[𝐧_\mathrm{𝟐}𝐩(nf)^2]^{1/2}`$ as a function of $`b`$ for several $`w`$’s, for $`n=3`$. In a similar way, in Figures (3, 4) we show $`f`$ and $`\sigma `$ for $`w=0.2`$ and $`n=1,2,\mathrm{},8`$.

## 3 Simulations: finite systems and finite range

The considerations in the previous sections apply in the limit $`m\mathrm{}`$. For any finite $`m`$ the only steady-state solution possible is the trivial one, $`f=0`$; with probability unity the infection will become extinct. However, the time to extinction may be very long, $`O(e^m)`$ if the system is above $`b_c`$. See, for example,  for the case $`n=1`$. For times short compared to the extinction time the system will be in a quasi-stationary state which resembles the infinite system.

We have done preliminary investigations of this and related behavior using computer simulations. We represent the population as follows: imagine a population of $`m`$ households each containing $`n`$ agents arranged on a ring. All of the members of the same household can infect one another with probability/unit time $`w`$. The communication with other households is variable range: we allow each household to infect all of the other households in a range $`L`$ on either side. See Figure (5). That is, the number of outside contacts is $`q=2Ln`$, and the infection probability is $`b/q`$. The equivalent of the random mixing case of the previous sections is $`2L=m1`$; in this case all agents are coupled to each other. The case $`n=1,L=1`$ is the well-known contact process .

Our computer algorithm is quite simple: we suppose that $`\gamma ,w,b`$ are all less than unity. We start by marking a certain fraction of agents as infected, which establishes $`I(0)`$ where $`I`$ is the total number infected in the population. For each computer time step we choose one agent, in household $`k`$, at random among the infected, and find all the other agents that are coupled to the one in question. These will be in households with index $`kL,kL+1,\mathrm{},k+L`$, where the indicies are to be interpreted with periodic boundary conditions. With probability $`w`$ or $`b/q`$ (depending whether the other agent is in $`k`$ or outside) the other agent is infected. Then another agent is chosen at random among the infected, and it is allowed to recover with probability $`\gamma `$. Each computer time step corresponds to an increment of ‘clock time’, $`dt=1/I`$. A typical result for $`I(t)`$ is shown in Figure(6). It is worth noting that this algorithm is quite fast – it is not difficult to treat as many as 3000 agents on a simple workstation.

In this section we will give a few results of the simulations. We have found that the SCF is an excellent representation of the behavior in large systems with long-range interactions, but that for smaller systems and shorter ranges the approximation becomes poorer, but still surprisingly good. A systematic investigation of these effects will be reserved for a future publication: here we give a few preliminary results.

### 3.1 Finite size effects

For large enough $`m`$ the simulation results agree very well with the SCF results. In Figure(6) the average of the time series agrees almost perfectly with the results from the previous sections: the numerical average of the simulation data differs from the SCF by less than the width of the line.

In addition we can try to understand the fluctuations of the time series by the following reasoning: if we imagine that the different households fluctuate more or less independently (this is the assumption under which the SCF method is valid), then the standard deviation of the total should be $`\sigma _{tot}=\sigma \sqrt{m}`$. As we see in the figure, this is a good approximation. However, we have observed that there is a large finite-size effect for this statistic. For $`m=301`$ the simulation result for $`\sigma _{tot}`$ is almost twice the prediction.

Another way to look at the standard deviation is to simply record the time series of the number infected in a given household, and directly find $`\sigma `$. We did this by averaging over 10 randomly chosen households in the system. We find very good agreement with the analytic treatment: for example, for $`w=0.2,b=1.,n=3,m=501`$ we find $`\sigma =0.82`$. The analytic result is $`0.84`$, well within our estimated error.

As we reduce $`m`$ we find a surprising result: within the numerical accuracy of our simulation we can find no change in $`f`$. Rather, at a certain small size, $`m=50`$ for the parameters in the last paragraph, the infection dies before we can collect adequate statistics.

### 3.2 Finite range effects

For the case of small $`L`$ the SCF ceases to be valid. In magnetic systems, locality of the interaction is well known to increase fluctuations and to change the threshold. We find related effects here, but the range needs to be very small indeed for the SCF result to be very inaccurate.

In Figure (7) we show simulation results for $`f(L)`$, the endemic level as a function of the range of the interaction between households, $`L`$. We show three different sets of parameters for $`n=3`$ and compared with the SCF ($`L\mathrm{}`$). The dependence is remarkable: for the parameter values ($`w=1.,b=1.`$), finite range has hardly any effect. The other two sets do show a decrease of $`f`$, and an apparent shift in the threshold. For $`w=0.2`$ $`b_c=0.67`$ in the SCF, well below $`b=1.`$ However, for $`L<5`$, the system is below threshold, and the epidemic dies quickly. A similar effect occurs for $`w=1.,b=0.3`$, which is below $`b_c`$ for $`L<10`$. (For large $`L`$, $`b_c=0.2`$ for these parameters.) Another way to put this is that there is a characteristic length $`L_c10`$ for these two sets of parameters.

These results are reminiscent of effects that occur in phase transitions in statistical physics. In fact, the contact process is fruitfully viewed as an example of a continuous non-equilibrium phase transition . Near the threshold (i.e. the critical point) of a continuous phase transition, the system develops a diverging correlation length. For total system sizes greater than this length finite size effects are small. Thus we might guess that the correlation lengths for two of the sets of parameters could be of order $`L_c`$. We have not investigated this matter in detail, but it would repay consideration.

## 4 Summary

In this paper we have given a number of explicit results for the SIS household model. Our point of view was to emphasize the SCF method as a guide to solution and insight. The numerical results give rise to a number of questions that are probably not accessible to rigorous proof, but for which numerical methods could be pursued. For example, we have speculated that there are critical fluctuations  in this system. That would mean that near threshold there would be large correlations in adjacent households for a finite range system, a result that could be of considerable interest.

There is also a practical aspect to this method. SCF theory is quite flexible, and allows considerable complication to be added to the model without substantially increasing the complexity of the solution. Suppose our household contains two types of members, for example, children who go to school and contract infections, and adults who do not. This would require a few changes: for a 3 person household with two adults and a child, the analog of the vector $`𝐏`$ would be of length 6 instead of 4, and the complete solution would be no worse than the numerical inversion of a $`5\times 5`$ matrix.

We would like to thank Frank Ball, Mark Newman, Carl Simon, and Jim Koopman for helpful discussions. GG was supported by the REU summer program of the Physics Department of the University of Michigan, and IMS would like to thank the Michigan Complex Systems program for support.


text:
# A Major SGR-like Outburst and Rotation Glitch in the No-Longer-So-Anomalous X-ray Pulsar 1E 2259+586

## 1 Introduction

Radiation from some isolated young neutron stars has been suggested to be powered by the decay of an ultra-high magnetic field (Thompson & Duncan, 1995, 1996). Only this “magnetar” model has been able to explain the properties of the rare class of soft gamma repeaters (SGRs), whose hallmark is repeated X-ray and gamma-ray bursts (e.g. Kouveliotou et al., 1998). A different rare class, Anomalous X-ray Pulsars (AXPs), are characterized by soft X-ray pulsations (Mereghetti & Stella, 1995), and have also been proposed as magnetars (Thompson & Duncan, 1995), from the similarity of their pulsations with those of SGRs in quiescence. However, competing models in which AXPs are accretion powered, via a fall-back disk rather than an unseen binary companion, have been proposed (Chatterjee et al., 2000). Recently, the detection of two weak X-ray bursts from the direction of AXP 1E 1048.1$``$5937 argued for AXPs being magnetars (Gavriil et al., 2002). That AXP had previously been identified as one that was most likely to burst on the basis of its unstable timing behavior which was reminiscent of that seen in SGRs, and because of its SGR-like spectrum (Kaspi et al., 2001).

1E 2259+586, a 7-s AXP in the supernova remnant CTB 109 (Fahlman & Gregory, 1981), in contrast to 1E 1048.1$``$5937, has shown remarkably stable timing behavior and pulsed X-ray fluxes in the past 5.6 yr (Kaspi et al., 1999; Gavriil & Kaspi, 2002). It also has the smallest inferred surface dipolar magnetic field of all AXPs (and SGRs for which it has been determined), and has an X-ray spectrum softer than those of the SGRs. Past observations have suggested, however, that the pulsar may experience epochs of activity, including flux, timing, and pulse profile variations (Iwasawa et al., 1992; Corbet et al., 1995; Baykal & Swank, 1996).

Here we report a major SGR-like outburst from the AXP 1E 2259+586, in which over 80 X-ray bursts were detected along with a variety of significant changes to the pulsed and persistent emission.

## 2 Observations and Results

The 1E 2259+586 outburst was detected in an observation that was made as part of a long-term Rossi X-ray Timing Explorer (RXTE) AXP monitoring program (Kaspi et al., 1999; Gavriil & Kaspi, 2002). Unexpectedly, bursts were seen during a 14.4 ks observation on June 18, 2002 (UT 15:39). The total on source exposure time was 10.7 ks. Data were taken with the Proportional Counter Array (PCA) in GoodXenonwithPropane mode, which records photon arrival times with 1-$`\mu `$s resolution, and bins photon energies into 256 channels. In subsequent analysis, photon arrival times at each epoch were adjusted to the solar system barycentre. The resulting time series were analysed in a variety of ways. Figure 1 shows the light curve binned with 125 ms time resolution, along with time series of several properties of the pulsed and persistent emission (see below). The decreasing burst rate and flux throughout our observation clearly indicates that we observed only the end of an event that commenced prior to the start of our observations. Lightcurves in the three operational PCUs look similar. Only the largest burst showed any excess in the PCA Standard 1 “Remaining Counts,” however, the flux correction due to deadtime is minimal ($``$10%).

To identify all bursts, the time series made using photons in the 2–20 keV range were searched for significant excursions from the mean count rate by comparing each time bin value with a windowed 7-s running mean. Bursts were identified assuming Poissonian statistics, and by combining probabilities from the separate PCUs. A detailed description of the burst search algorithm and subsequent analysis will be presented elsewhere. Here we summarize the main results. The bursts range in duration from 2 ms to 3 s and show a variety of structures. A few ($``$10) have clear fast-rise, exponential-tail morphology, however the rest appear approximately symmetric. Their spectra are all well modeled by simple power laws, with indices in the range $``$0.25 to $``$2.4. Burst fluences in the 2–10 keV band are in the range $`3\times 10^{11}`$ to $`5\times 10^9`$ erg cm-2, corresponding to energies $`3\times 10^{34}`$ to $`5\times 10^{36}`$ erg, assuming isotropic emission and a distance of 3 kpc to the source (Kothes et al., 2002). The sum total of all burst fluences is $`3.2\times 10^8`$ erg cm-2, corresponding to energy $`3.4\times 10^{37}`$ erg (2–10 keV). Peak burst fluxes in the 2–10 keV band are in the range $`1\times 10^9`$ to $`4\times 10^7`$ erg s-1 cm-2, corresponding to peak luminosities of $`1\times 10^{36}`$ to $`4\times 10^{38}`$ erg s-1.

Follow-up RXTE observations on June 20 revealed no further bursts, nor have any of the 15 observations, each of duration 7–8 ks, obtained every $``$10 days since. Neither target-of-opportunity observations obtained with the XMM-Newton satellite on June 21, nor XMM-Newton observations scheduled fortuitously 7 days prior to the burst, revealed any additional bursts (work in preparation). The RXTE All Sky Monitor observed the field on June 18 at UTs 03:50 and 14:43 for $``$90 s per observation but detected no enhanced flux, with 99% confidence upper limits of $`1\times 10^9`$ erg s-1 cm-2 (2–10 keV).

In spite of the large (1) field of view of the PCA, we are certain that the AXP is the origin of the bursts, as many properties of the pulsed emission were simultaneously observed to change dramatically (Fig. 1). The persistent flux evolution was determined as follows. A spectral analysis was done using the XSPEC software package v11.2.0 111http://xspec.gsfc.nasa.gov, in which the preburst PCA data for 1E 2259+586 were modeled using the best available background models made with the FTOOL pcabackest (http://heasarc.gsfc.nasa.gov/docs/xte/recipes/p2.html), the pulsar spectrum as determined using a XMM-Newton observation of the pulsar made 1 week before the burst observation (work in preparation), and an additional component to account for the remaining emission in the PCA field-of-view. The burst data were modeled by a blackbody plus power-law component, while holding the pulsar equivalent neutral hydrogen column density and the remaining emission model fixed. The resulting persistent fluxes are show in Figure 1 (second panel from the top). The pulsed flux evolution, also show in in Figure 1 (second panel from the top) was calculated by first folding $``$200 s long data segments with the spin ephemeris (Table 1), then summing the first six harmonics of the normalized Fourier powers of the resulting pulse profiles. The total 2–10 keV fluence over and above the quiescent flux is $`2\times 10^6`$ erg cm-2, two orders of magnitude above that in the bursts. As is also seen in Figure 1, the pulsar spectrum clearly hardened during the outburst, and relaxed back toward the quiescent spectral paramters during the course of the observation. The fitted blackbody radius remained approximately constant throughout. Interestingly, the same cannot be said of the ratio of power-law to blackbody flux; as is seen in the bottom panel of Figure 1, the latter continued evolving away from the quiescent state during our observation.

A significant change in the pulse morphology was observed at the burst epoch, as shown in Figure 2. During the outburst, the amplitudes of the peaks relative to the pre- and post-outburst profiles are clearly reversed. The relative phase displayed above is that successfully used in our timing analysis. The profile change is similar in different energy bands. This different pulse profile persisted for at least 2 days following the outburst, and gradually returned to its pre-outburst morphology after $``$6 days (Gavriil et al. in preparation).

The star underwent a sudden spin up or “glitch” at the outburst epoch. For details regarding how the timing analysis was done, see Gavriil & Kaspi (2002). For the burst and immediate post-outburst data, the pulse profile changes described above resulted in obvious phase jumps corresponding to the two peaks being swapped in the cross-correlation. The glitch epoch was determined by requiring zero phase jump between pre- and post-outburst ephemerides. The RXTE data obtained during and after the burst are well characterized (rms residual 1.5% of the period for the full data set) by $`\mathrm{\Delta }\nu /\nu =(4.10\pm 0.03)\times 10^6`$, similar to that observed in radio pulsar glitches (Lyne & Smith, 1990). The best-fit glitch epoch is consistent at the $`<1\sigma `$ level with having occurred during our observation. Additionally, the spin-down rate can be modeled as having approximately doubled abruptly. Precise spin parameters are given in Table 1. Residuals from the last $``$100 days of timing suggest that the spin-down rate may have relaxed back to near its pre-burst value by $``$60 days post-outburst, however additional observations are required to confirm this.

Target-of-opportunity near-infrared observations were made using the NIRI instrument at the 8-m Gemini-North telescope in Hawaii on June 21 at UT 14:44 using a $`K_s`$ filter (0.15 $`\mu `$m wide centred on 2.15 $`\mu `$m). The observation had total exposure 1530 s with 0.7′′ seeing and light cirrus. The data were reduced using the Gemini IRAF package and photometry performed using standard IRAF procedures. The proposed infrared counterpart (Hulleman et al., 2001) of 1E 2259+586 had magnitude $`20.36\pm 0.15`$, $`1.33\pm 0.22`$ mag (factor of $`3.40_{0.62}^{+0.77}`$) brighter 3 days after the outburst than was observed in a 2000 Keck telescope observation (Hulleman et al., 2001). A second Gemini/NIRI observation was obtained on June 28 at UT 14:51, with 900 s of exposure and 0.55′′ seeing. This time, the AXP counterpart had faded to magnitude $`21.14\pm 0.21`$, for a difference relative to the 2000 Keck observation of $`0.56\pm 0.29`$ mag (factor of $`1.67_{0.39}^{+0.52}`$ in brightness). Photometric measurements on 7 reference objects in the field agreed with those obtained at Keck to within 0.007 mag and 0.028 mag for the first and second nights, respectively.

Target-of-opportunity radio observations were also made using the Very Large Array in New Mexico on June 20, 2002. The 1E 2259+586 field was observed for 2420 s in B array at a central observing frequency of 1424.3 MHz. After standard calibration, imaging and cleaning using the MIRIAD software package, an rms noise level of 15 $`\mu `$Jy/beam ($`4.6^{\prime \prime }\times 3.9^{\prime \prime }`$ beam) was achieved. No emission was detected. We place a 3$`\sigma `$ upper limit of $`50`$ $`\mu `$Jy on the radio flux for this epoch.

## 3 Discussion

The X-ray phenomenology we have observed in this major AXP outburst is all reminiscent of that seen in SGR bursts. The short bursts (Fig. 1) are very similar to short SGR bursts. The long, thermally evolving tail is similar to that seen in a handful of SGR bursts (Lenters et al., 2003). A timing anomaly in SGR 1900+14 was seen at the time of the giant flare in 1998 (Woods et al., 2002), as was a pulse profile change and enhanced pulsed and persistent flux (Göğüş et al., 2001). Thus, this AXP has shown uniquely SGR-like bursting behavior. 1E 2259+586 showed the most stable timing behavior of all AXPs in the 5.6 yr prior to this event (Gavriil & Kaspi, 2002), while 1E 1048.1$``$5937, the only other AXP seen to burst (Gavriil et al., 2002), showed the least stable behavior, as well as the hardest AXP spectrum (Kaspi et al., 2001). Thus it seems any AXP can burst.

The properties of the outburst solve a number of previously outstanding AXP problems. A similar pulse profile change was claimed previously in data for 1E 2259+586 from the Ginga mission in 1989 (Iwasawa et al., 1992). The archival Ginga data show no evidence of bursts. The Ginga observation probably took place just after an outburst, consistent with the reported timing anomaly at the same epoch (Iwasawa et al., 1992). This suggests that such outbursts occur on decade time scales (Heyl & Hernquist, 1999). In addition, previously reported large X-ray flux variations in 1E 2259+586 and 1E 1048.1$``$5937 (Iwasawa et al., 1992; Baykal & Swank, 1996; Oosterbroek et al., 1998) that were called into question by the flux stability observed in the $``$5 yr prior to June 2002 (Gavriil & Kaspi, 2002; Tiengo et al., 2002) are now more understandable as enhanced emission due to bursting episodes. Pulsations from the AXP candidate AX J1845$``$0258 have been detected only once, in spite of multiple observations (Vasisht et al., 2000). This may have been following a similar outburst (though there is no evidence for bursts in the archival ASCA data).

The 1E 2259+586 outburst likely resulted from a sudden event in the stellar crust, such as a crustal fracture, which simultaneously affected both the superfluid interior and the magnetosphere. The large spin-up can be explained by the coupling of the faster-rotating superfluid inside the star with the crust, following the unpinning of angular momentum vortices from crustal nuclei. The fractional frequency increase is similar to that observed in many radio pulsars (Lyne & Smith, 1990) but is smaller than could have been detected in SGR timing data. However, the factor of two increase in the spin-down rate is unprecedented for radio pulsars, though possibly not for SGRs (Woods et al., 2002). According to glitch theory (Alpar et al., 1993), when the glitch occurs, a portion of the superfluid decouples, decreasing the effective moment of inertia of the star. For fixed external torque, an increase in spin-down rate results. For radio pulsars, the decoupled portion amounts to $``$1% of the stellar moment of inertia, corresponding to the observed $``$1% increases in spin-down rates (Alpar et al., 1993). For 1E 2259+586, however, most of the stellar moment of inertia would have had to decouple. This could imply a decoupling of core, as opposed to crustal, superfluid.

Alternatively, the external torque could have changed, due to a restructuring of the magnetosphere. Indeed the enhanced X-ray luminosity is too large to be explained as energy dissipated by vortex unpinning (Thompson & Duncan, 1996) or crustal elastic energy (Ruderman, 1991). A decaying magnetar-strength magnetic field can cause severe stress on the crust. A large-scale fracture could trigger vortex unpinning, and, simultaneously, shift magnetic field footpoints, resulting in a magnetospheric reconfiguration (Thompson & Duncan, 1995). The pulse profile variation is unlikely to be a result of the change in magnetospheric structure (Thompson et al., 2002), since the torque change is much longer-lived. Rather, the profile change probably occurred at the surface; the effective blackbody radius of $``$1 km as determined from the spectral fits supports a localized enhancement.

Notably, there is no evidence for an accompanying giant soft gamma-ray flare, as might be expected from a sudden restructuring of the surface magnetic field of a magnetar (Woods et al., 1999). From the Interplanetary Network spacecraft, an upper limit on the fluence of a soft gamma-ray flare from 1E 2259+586, near the time of the X-ray outburst, is $`5\times 10^7`$ erg cm-2 (25-150 keV) on time scales 0.25-0.5 s (K. Hurley 2002, personal communication). This corresponds to an energy $`5\times 10^{38}`$ erg, six orders of magnitude below that released in the giant SGR flares. This is consistent with the absence of any radio emission post-outburst from 1E 2259+586, as well as with the absence of bright soft gamma-ray flares from this source in the past. Thompson & Duncan (1996) showed that the absence of a large flare in the event of a glitch requires the neutron star crust to have deformed plastically. This demands a magnetic field roughly two orders of magnitude greater than that implied by the spin-down of 1E 2259+586. This could be explained by higher order multipole moments which are negligible at the light cylinder, where the spin-down torque arises.

The characteristic age ($`P/2\dot{P}100200`$ kyr) of 1E 2259+586 is much larger than the inferred age ($`10`$ kyr) of the supernova remnant CTB 109 in which it resides (Rho & Petre, 1997). It is tempting to explain this discrepancy as being due to the pulsar having episodes of transient accelerated spin-down such as we observe post-outburst. However at least in this instance, the increased spin-down rate could be roughly compensated by the sudden spin-up.

The near-infrared enhancement post-outburst is intriguing. Currently, the magnetar model does not address the origin of such emission. In conventional rotation-powered pulsars, infrared emission is thought to arise from a population of synchrotron radiating electron/positron pairs in the outer magnetosphere (see Lyne & Smith (1990) and references therein). An enhancement is therefore consistent with a change in the magnetospheric field structure suggested by the torque change. Future observations can test this by comparing infrared variations to the torque evolution.

We thank M. Lyutikov, R. Manchester, and C. Thompson for useful discussions. We thank F. Nagase and B. Paul for help with the archival Ginga data. This work was supported in part by NSERC, CIAR, and NASA/LTSA and made use of data obtained through the HEASARC Online Service, provided by the NASA/Goddard Space Flight Center. We thank J. Swank and the RXTE scheduling team. The Gemini Observatory is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the Particle Physics and Astronomy Research Council (United Kingdom), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), CNPq (Brazil), and CONICET (Argentina). We thank F. Jansen for scheduling Target-of-Opportunity observations with XMM-Newton and B. Clark for scheduling the Ad Hoc observation with the VLA. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.