text:
# In Hot Pursuit of the Hidden Companion of 𝜂 Carinae: An X-ray Determination of the Wind Parameters

## 1 Introduction

The superluminous star $`\eta `$ Carinae (HD 93308, HR 4210) continues to be extensively studied over a host of different wavelengths, yet remains intriguingly enigmatic. It is amongst the most unstable stars known. In the 1840s, and again in the 1890s, it underwent a series of giant outbursts (e.g. Viotti V1995 (1995)) which ejected large masses of material into the surrounding medium. HST images of the resulting bipolar nebula (e.g. Morse et al. MDB1998 (1998)), known as the Homunculus, show it to be amongst the most spectacular in our Galaxy. The central object is now largely obscured by dust, and the cause of the outbursts and the nature of the underlying star (at outburst and today) remain speculative. The source continues to show brightness fluctuations and emission-line variations. Further details of $`\eta `$ Car can be found in the review by Davidson & Humphreys (DH1997 (1997)).

In recent years evidence for binarity in this system has been accumulating. Damineli (D1996 (1996)) first noted a 5.5 yr period in the variability of the He I 10830 $`\mathrm{\AA }`$ line. Further photometric and radial velocity studies (Damineli et al. DCL1997 (1997), D2000 (2000)), X-ray observations (Tsuboi et al. TKSP1997 (1997); Corcoran et al. CFPSD2000 (2000) and references therein), and radio data (Duncan et al. D1995 (1995), DWRL1999 (1999)) have supported the 5.5 yr period and the binary hypothesis. However, the ground based radial velocity curve was not confirmed by higher resolution spectra with STIS, indicating that at least the time of periastron passage is not well defined by the UV and optical spectra (Corcoran et al. 2001a ). A comparison of the abundances from the central object(s) and the composition of the Homunculus nebula has also determined that there are at least two stars in this system (Lamers et al. LLPW1998 (1998)).

$`\eta `$ Car is often classified as a luminous blue variable (LBV). These are massive stars believed to be in a rapid and unstable evolutionary phase in which many solar masses of material are ejected into the interstellar medium over a relatively short period of time ($`10^4`$ yr). LBVs are regarded as a key phase in the evolution of massive stars, during which a transition into a Wolf-Rayet star occurs (e.g. Langer et al. L1994 (1994); Maeder & Meynet MM2001 (2001)). Due to their rarity and complex nature however, we unfortunately still have no definitive theory for mass-loss during the LBV stage. The majority of proposed mechanisms to drive LBV instabilities, the onset of higher mass-loss rates and underlying eruptions, are concerned with the importance of radiation pressure within the outer envelope of the LBV, and for example utilize pulsational instabilities (e.g. Guzik et al. GCD1999 (1999)), dynamical instabilities (e.g. Stothers & Chin SC1993 (1993)), or presuppose Eddington-like instabilities. The latter could arise from an enhancement in opacity as the star moves to lower temperatures (e.g. Lamers Lam1997 (1997)), or from the influence of rotation (e.g. Langer Lan1997 (1997); Zethson et al. Z1999 (1999)). Alternatively, the possibility that binarity plays a fundamental role in explaining observed LBV outburst properties has also been considered (Gallagher G1989 (1989)), though most LBVs are not known binaries. Clearly, determining the wind and stellar properties of LBV stars is paramount (see, for example, the discussions in Leitherer et al. L1994 (1994) and Nota et al. N1996 (1996)).

An important question is the degree to which binarity influences the properties of LBVs (i.e. do LBVs in binaries evolve differently than single LBVs?). So while the presence of a companion can be exploited to help measure the mass of such stars, we must bear in mind that binary LBVs and single LBVs may be quite distinct objects. Therefore, in order to use $`\eta `$ Car to understand some of the defining LBV characteristics such as their extremely high mass-loss rates, we first need to determine beyond all doubt that $`\eta `$ Car is in fact a binary, and then to determine the influence of the companion on the system.

Investigations over the last few years have already helped to form a basic picture of $`\eta `$ Car. The orbital parameters, although uncertain, indicate the presence of an early-type companion star, which will also have a powerful stellar wind. In such binaries, a region of hot shocked gas with temperatures in excess of 10 million K is created where the stellar winds collide (Prilutskii & Usov PU1976 (1976); Cherepashchuk C1976 (1976)). The wind-wind collision (WWC) region is expected to contribute to the observed emission from this system, particularly at X-ray and radio wavelengths. Previous X-ray observations revealed extended soft emission from the nebula and strong, hard, highly absorbed, and variable emission closest to the star (Corcoran et al. C1995 (1995); Weis et al. WDB2001 (2001)), in contrast to the emission characteristics from single stars, which are typically softer, much less absorbed, substantially weaker and relatively constant in intensity. Since 1996 Feb $`\eta `$ Car has been continuously monitored by RXTE in the 2–10 keV band (e.g. Corcoran et al. 2001a ). The lightcurve (Fig. 1) contains remarkable detail showing a slow, almost linear, rise to maximum over a period of $`1`$ yr, followed by a rapid drop to approximately 1/6 of the peak intensity for $`3`$ months, an almost as sharp rise to approximately 1/2 of the peak intensity level, and then almost constant intensity for $`3/4`$ of the proposed 5.5 yr orbital period. The drop to minimum was successfully predicted from numerical models of the WWC (Pittard et al. PSCI1998 (1998)) before being actually observed.

Small scale quasi-periodic outbursts in the X-ray lightcurve have also been detected (Corcoran et al. CISDPS1997 (1997)). Estimates of changes in the timescale between successive flares as a function of phase were made by Davidson et al. (DIC1998 (1998)) for a variety of assumed orbital elements. RXTE X-ray observations obtained after the X-ray minimum seem to show a lengthening of the flare timescale (Ishibashi et al. I1999 (1999)), which indirectly support the binary model. The latest published X-ray observation of $`\eta `$ Car is of a high resolution grating spectrum taken with the Chandra X-ray observatory (Corcoran et al. 2001b ). Preliminary analysis has revealed the presence of strong forbidden line emission, which suggests that the density of the hot gas, $`n_\mathrm{e}<10^{14}\mathrm{cm}^3`$. This can be contrasted with the newly published X-ray grating spectra of the single stars $`\theta ^1`$ Ori C (Schulz et al. SCHL2000 (2000)), $`\zeta `$ Ori (Waldron & Cassinelli WC2001 (2001)), and $`\zeta `$ Pup (Kahn et al. K2001 (2001)), all of which have weak forbidden lines (indicative of either high densities or high UV flux near the line forming region). This lends further support to a wind-wind collision model, although it is possible that some of the forbidden emission may be related to the surrounding nebula.

If $`\eta `$ Car is in fact a binary system, the orbital elements and stellar parameters are not yet tightly constrained by current observations. Ground-based observations (for which good phase coverage exist) are hampered by poor spatial resolution and thus suffer contamination from strong nebular emission. High-spatial resolution spectra have been obtained by STIS but phase coverage is currently very limited. HST STIS observations at two different phases of the 5.52 year cycle (Davidson et al. DI2000 (2000)) did not confirm the predicted variations in the radial velocity of the emission lines based on the ground-based radial velocity curve (Damineli et al. D2000 (2000)). If $`\eta `$ Car is a binary, it is vitally important to determine the stellar parameters of the companion so that the effect of the companion on observations can be understood and the correct stellar parameters of the primary can be derived.

### 1.1 X-ray Emission from Colliding Winds

The wealth of information contained in X-ray spectra of colliding wind binaries (e.g. the density, temperature, velocity, abundance, and distribution of the shocked gas in the wind collision region) has been a strong motivating force for observers and theorists alike in this field. Since the hot plasma in most colliding wind shocks is optically thin and collisionally ionized, and is generally assumed to be in collisional equilibrium, Raymond-Smith (Raymond & Smith RS1977 (1977)) or MEKAL (Mewe et al. MKL1995 (1995)) spectral models are normally fitted to such data (e.g. Zhekov & Skinner ZS2000 (2000); Rauw et al. RSPC2000 (2000); Corcoran et al. 2001b ). However, the multi-temperature, multi-density nature of the WWC region means that at best simple fits with one- or two-temperature Raymond-Smith models can only characterize the broad properties of the emission. In this way one can estimate an ‘average’ temperature of the shocked gas, and an ‘average’ absorbing column, but little is learned of the underlying stellar wind parameters. At worst the application of one- or two-temperature models to what is inherently multi-temperature emission can lead to spurious values of some of the fit parameters, e.g. abundances (cf. Strickland & Stevens SS1998 (1998)).

Complex numerical hydrodynamical models have often been applied to gain insight into colliding wind systems (e.g. Stevens et al. SBP1992 (1992), Pittard et al. PSCI1998 (1998)). However, while undoubtedly useful, their interpretation can be difficult, and to date there have been only two published papers where observed spectra are directly fitted with synthetic spectra from such models. In the pioneering work of Stevens et al. (S1996 (1996)), medium resolution ASCA spectra of the Wolf-Rayet binary $`\gamma ^2`$ Velorum were fitted against a grid of synthetic spectra. In this fashion they were able to obtain direct estimates of the mass-loss rates and terminal velocities of the individual stellar winds. As mass-loss rates obtained from measures of radio flux or spectral line fits depend on a variety of untested assumptions, the importance of a new independent method to complement estimates from free-free radio or sub-mm observations, or from H$`\alpha `$ or UV spectral line fitting, cannot be stressed enough. Rates derived by Stevens et al. (S1996 (1996)) with this new method were significantly lower than the commonly accepted estimates for $`\gamma ^2`$ Velorum based on radio observations, but an indication of the future benefits of this method was realized when both sets of estimates were later brought into agreement following a surprisingly large reduction in the distance to this star from Hipparcos data111The thermal radio flux, $`S_\nu \dot{M}^{4/3}D^2`$ (where $`D`$ is the distance to the source) whereas the X-ray flux from an adiabatic wind collision is $`F_\nu \dot{M}^2D^2`$. Therefore $`\dot{M}_{\mathrm{radio}}D^{3/2}`$ whereas $`\dot{M}_{\mathrm{xray}}D`$. If $`D`$ is revised downwards, $`\dot{M}_{\mathrm{radio}}`$ decreases faster than $`\dot{M}_{\mathrm{xray}}`$.. We note that this method can also provide insights into the values of parameters which are otherwise virtually impossible to estimate, such as the mass-loss rate of the companion, $`\dot{M}_2`$, or the characteristic ratio of the pre-shock electron and ion temperature (Zhekov & Skinner ZS2000 (2000)). The quality of recently available X-ray grating spectra now gives us access to important X-ray emission line diagnostics which should severely constrain models of the X-ray emission distribution. This means that stellar wind parameters can in principle be reliably estimated from analysis of X-ray grating spectra of colliding wind binaries.

In addition to testing the binary hypothesis, the Chandra grating spectrum of $`\eta `$ Car (Corcoran et al. 2001b ) provides the ideal opportunity to test the method developed by Stevens et al. (S1996 (1996)) against a spectrum of much higher spectral resolution, and to pin down important physical parameters of the system. In this paper we therefore fit the X-ray grating spectrum using a grid of colliding wind emission models to i) test the binary hypothesis, and ii) to attempt to obtain accurate estimates of the wind parameters of each star. The fact that we are able to obtain good fits, with sensible model parameters, gives us further confidence in the binary hypothesis. We also find that unlike the UV and optical work where fits to the primary are made difficult by significant contamination from the companion star, instead the X-ray emission arises from the shocked wind of the companion and suffers essentially zero contamination from the wind of the primary. Therefore the X-ray data uniquely samples parameters of the companion, in contrast to the optical analysis which probes the nature of the primary. In this sense our analysis is entirely complementary to the complex fits to the UV and optical HST spectrum of $`\eta `$ Car by Hillier et al. (HDIG2001 (2001)). Our analysis also provides us with a new estimate of the mass-loss rate of the primary star. For details of the Chandra observation and an initial analysis of the data the reader is referred to Corcoran et al. (2001b ). Here we only note that there is no significant contamination of the dispersed spectrum from any spatially resolved emission (i.e. in the Homunculus). In Section 2 we discuss the creation and variation of the model spectral grid; in Section 3 we describe the fit results; and in Section 4 we summarize and conclude.

## 2 The Synthetic Colliding Wind Spectra

The method applied by Stevens et al. (S1996 (1996)) as discussed in Section 1 is as follows. We first calculated a whole host of hydrodynamical models with different stellar wind parameters. We then generated a synthetic X-ray spectrum from each model. Finally, we fit the grid of synthetic spectra to the actual data. Best-fit values are then obtained for each parameter of interest (e.g. $`\dot{M}_1`$, $`\dot{M}_2`$, $`v_\mathrm{}_1`$, $`v_\mathrm{}_2`$ etc. ). Specific points in relation to this method are highlighted below.

### 2.1 Hydrodynamical Models of $`\eta `$ Carinae

The colliding wind models were calculated using VH-1 (Blondin et al. BKFT1990 (1990)), a Lagrangian-remap version of the third-order accurate Piecewise Parabolic Method (PPM; Colella & Woodward CW1984 (1984)). The stellar winds are modelled as ideal gases with adiabatic index $`\gamma =5/3`$. Radiative cooling is included via the method of operator splitting, and is calculated as an optically thin plasma in ionization equilibrium. The cooling curve for the temperature range $`4.0<\mathrm{log}T<9.0`$ was generated using the Raymond-Smith plasma code. Despite some evidence to the contrary (see Corcoran et al. 2001b ), we have used solar abundances throughout this work, since it keeps the problem as simple as possible in this first investigation. We aim to model non-solar abundances in a later paper.

The simulations were calculated assuming cylindrical symmetry of the wind collision zone - orbital effects are negligible at the phase of the Chandra observation ($`\varphi =0.60`$). Grid sizes spanned a range from $`300\mathrm{x}\mathrm{\hspace{0.33em}300}`$ to $`550\mathrm{x}\mathrm{\hspace{0.33em}450}`$ cells. Each grid cell was square and of constant size on an individual grid. The linear dimension of each cell was either $`2.67\times 10^{12}`$ or $`4.0\times 10^{12}`$ cm. The maximum distance from the axis of symmetry was $`1.2\times 10^{15}`$ cm in all cases. At the phase of the observation, the orbital separation of the stars using the ephemeris of Corcoran et al. (2001a ) is $`4\times 10^{14}`$ cm. Thus there is either 100 or 150 grid cells between the stars. The large separation means that the winds are likely to collide at very near their terminal velocities. We therefore assume that we can treat the winds as being instantaneously accelerated to their terminal velocities, and do not consider any radiative driving effects. In comparison with the work on $`\gamma ^2`$ Velorum (Stevens et al. S1996 (1996)), assuming terminal velocity winds and negligible binary rotation is more valid for $`\eta `$ Car. We further assume that the winds are spherically symmetric.

To compute a grid of synthetic spectra we initially varied four parameters ($`\dot{M}_1`$, $`\dot{M}_2`$, $`v_\mathrm{}_1`$, $`v_\mathrm{}_2`$) in the hydrodynamic models. During our investigation we also included the separation of the stars, $`D`$, as an additional free parameter. From fits at this point it was found that $`D4\times 10^{14}`$ cm with a small uncertainty, in good agreement with the expected value from our current understanding of the orbit. We therefore fixed it at this value for the rest of our analysis. However, we found that there were large uncertainties on the values of $`\dot{M}_1`$ and $`v_\mathrm{}_1`$. With the benefit of hind-sight this is not too surprising given the known slow speed of the primary wind ($`v_\mathrm{}_1500700\mathrm{km}\mathrm{s}^1`$), which means that the shocked primary wind is not a strong source of X-ray emission at hard energies. Therefore for our final grid we fixed the terminal velocity of the primary star at $`500\mathrm{km}\mathrm{s}^1`$ and adjusted $`\dot{M}_1`$ to obtain a desired value for the wind momentum ratio, $`\eta =\dot{M}_2v_\mathrm{}_2/(\dot{M}_1v_\mathrm{}_1)`$, with $`\dot{M}_2`$ and $`v_\mathrm{}_2`$ as free parameters.

The range of the free parameters used in our final grid is given in Table 1. To restrict the number of models to a manageable number the parameter steps are fairly coarse. In future papers we will use a much finer grid. The range in the value of $`\eta `$ corresponds to either the wind of the primary ($`\eta =0.1`$) or the secondary ($`\eta =5`$) dominating, or the winds having equal momentum fluxes ($`\eta =1`$). The distance of the stagnation point from the centre of the primary star is given by

$$r=\frac{1}{1+\sqrt{\eta }}D,$$
(1)
and ranges from $`1.2\times 10^{14}`$ cm to $`3.0\times 10^{14}`$ cm. In the most extreme cases, this is still of order ten to a hundred times the radius of the stars, and justifies our use of terminal velocity winds. The kinetic power of the primary wind ranged from $`4.7\times 10^{34}\mathrm{erg}\mathrm{s}^1`$ to $`7.9\times 10^{38}\mathrm{erg}\mathrm{s}^1`$, and from $`7.1\times 10^{35}\mathrm{erg}\mathrm{s}^1`$ to $`7.9\times 10^{38}\mathrm{erg}\mathrm{s}^1`$ for the secondary wind. The combined kinetic power of the winds ranges from $`7.6\times 10^{35}\mathrm{erg}\mathrm{s}^1`$ to $`1.6\times 10^{39}\mathrm{erg}\mathrm{s}^1`$. The half-opening angle of the contact discontinuity measured from the line between the secondary star and the shock apex ranges from $`50120^{}`$ (cf. Eichler & Usov EU1993 (1993)).

The effect of radiative cooling can be quantified by the parameter $`\chi `$, the ratio of the cooling timescale to the dynamical timescale of the system. For shocked gas near the local minimum in the cooling curve at $`T2\times 10^7`$ K,

$$\chi \frac{v_8^4d_{12}}{\dot{M}_7},$$
(2)
where $`v_8`$ is the wind velocity in units of $`1000\mathrm{km}\mathrm{s}^1`$, $`d_{12}`$ is the distance to the contact discontinuity in units of $`10^{12}`$ cm, and $`\dot{M}_7`$ is the mass-loss rate in units of $`10^7\mathrm{M}_{}\mathrm{yr}^1`$ (cf. Stevens et al. SBP1992 (1992)). This equation is valid for post-shock temperatures in the range $`10^7<T<10^{8.2}`$ K ($`680/\sqrt{\mu }<v<2700/\sqrt{\mu }\mathrm{km}\mathrm{s}^1`$ where $`\mu \mathrm{m}_\mathrm{H}`$ is the average particle mass in grammes). For material with solar abundance ($`\mu =0.6`$), this corresponds to $`0.9<v_8<3.5`$. For winds with slower velocities, the shocked gas lies on the negative slope of the cooling curve. In this case the dependence of $`\chi `$ on the velocity of the wind is steeper:

$$\chi \frac{v_8^{5.2}d_{12}}{\dot{M}_7}.$$
(3)
In our grid of simulations, the shocked wind of the primary star is almost always strongly radiative ($`\chi _{1,\mathrm{min}}<<1.0`$) while the secondary’s is never strongly radiative ($`\chi _{2,\mathrm{min}}=2.0`$). However, for some choices of the 3 parameters in Table 1, the primary’s wind can approach the point where it starts to become adiabatic ($`\chi _{1,\mathrm{max}}1.0`$). In Fig. 2 we show two examples of the hydrodynamic calculations used to calculate synthetic spectra. Due to the difference in terminal velocity of the two winds, the velocity shear at the wind interface always generates Kelvin-Helmholtz instabilities in our models. The rapid cooling of at least one of the two winds means that thin-shell instabilities also always occur in our models. The colliding winds region is most unstable when both shocked winds rapidly cool (cf. Stevens et al. SBP1992 (1992)).

The spectrum from each hydrodynamic model was averaged over 3 ‘snapshots’ each spaced by about 70 d. This is of order the wind dynamical timescale and serves to approximate a time-averaged spectrum. It is particularly important to adopt this approach for models where the wind collision region is very unstable since in these cases the flux can vary by greater than $`\pm 50\%`$ from one snapshot to the next. As already noted by Corcoran et al. (2001b ) no variability is seen during the $`90\mathrm{ksec}`$ exposure of the Chandra grating spectrum. This is to be expected, however, since the dynamical timescale is much longer. Two examples of ‘time-averaged’ synthetic spectra are shown in Fig. 2, next to density plots of the corresponding hydrodynamic calculation.

The synthetic spectra were calculated over the energy range 1.26–10 keV (10–1.26 $`\AA `$ respectively), and have a resolution of 0.005 $`\AA `$. The actual grating spectrum shows significant absorption below 1.5 keV and contains essentially no useful information below our lower limit of 1.26 keV. The spectral resolution of our synthetic spectra is approximately twice as high as the grating data that we model. Use of the Raymond-Smith plasma code implicitly assumes that the plasma is in thermal equilibrium. This is generally true for colliding wind binary systems (see Luo et al. LMM1990 (1990)) and is a good assumption for $`\eta `$ Car.

### 2.2 Variations with $`\eta `$, $`\dot{M}_2`$, and $`v_\mathrm{}_2`$

In Fig. 3 and the following subsections we reveal how the shape and normalization of the calculated spectra depend on the various free parameters of the grid.

#### 2.2.1 Variation with $`\eta `$

Due to the low preshock velocity of the primary wind, the X-ray emission from the WWC is almost entirely from the shocked secondary wind. A greater fraction of the secondary wind passes into the wind collision zone as the relative momentum of the primary wind is increased, so we expect the luminosity to increase as $`\eta `$ decreases, as indeed the top panel of Fig. 3 shows. An important point, however, is that the shape of the spectra do not seem dependent on the value of $`\eta `$, contrary to our expectations before starting this study.

#### 2.2.2 Variation with $`\dot{M}_2`$

As we increase the mass-loss rate of the secondary star in our models, the kinetic power and the density of its wind also increase. Therefore, for a given wind momentum ratio, the rate of conversion of kinetic energy to thermal energy increases, as does the density of the postshock gas. Thus there is more available energy to radiate and a greater efficiency of radiation. Hence we should expect the X-ray luminosity to scale strongly with increasing $`\dot{M}_2`$. This is indeed what we see in the middle panel of Fig. 3, where an order of magnitude increase in $`\dot{M}_2`$ leads to an almost two order of magnitude increase in X-ray luminosity, indicating that the calculations in this panel are in the adiabatic regime ($`L_x\dot{M}^2`$; see Stevens et al. SBP1992 (1992)). However, there is no evidence for a softening of the spectrum with increasing $`\dot{M}_2`$. This is again contrary to our initial expectations (note, however, the discussion in the next section).

#### 2.2.3 Variation with $`v_\mathrm{}_2`$

If the preshock velocity of a wind is increased, the postshock temperature increases ($`Tv^2`$). Therefore we expect the model spectra to harden as the preshock velocity of the secondary wind is increased. Because the kinetic power of the wind increases we also expect the luminosity to increase. The variation of the spectrum with $`v_\mathrm{}_2`$ is shown in the bottom panel of Fig. 3. Here we find that our expectations are only partially met. The overall luminosity does increase with $`v_\mathrm{}_2`$, and the spectrum does harden between $`v_\mathrm{}_2=15003000\mathrm{km}\mathrm{s}^1`$, but it softens between $`v_\mathrm{}_2=30005000\mathrm{km}\mathrm{s}^1`$. An investigation into this surprising finding offers the following explanation.

For fixed $`\eta `$ and $`\dot{M}_2`$, and variable $`v_\mathrm{}_2`$, the density of the wind varies as $`\rho \dot{M}/v_\mathrm{}_2`$. If $`v_\mathrm{}_2`$ is increased $`\rho `$ decreases, which reduces the radiative efficiency of the hot shocked gas near the stagnation point ($`\dot{E}n^2`$), which acts to slightly suppress the hard flux. A plot of the average density against $`T`$ for $`v_\mathrm{}_2=1500,\mathrm{\hspace{0.33em}3000},\mathrm{\hspace{0.33em}5000}\mathrm{km}\mathrm{s}^1`$ also reveals that for $`10^{6.3}<T<10^{7.4}`$ the average density is higher when $`v_\mathrm{}_2=5000\mathrm{km}\mathrm{s}^1`$ than when $`v_\mathrm{}_2=3000\mathrm{km}\mathrm{s}^1`$, which tends to keep up the soft emission relative to the hard emission. This appears to be because there is some mixing at the interface between the primary and secondary wind, and the shocked primary wind is denser when $`v_\mathrm{}_2=5000\mathrm{km}\mathrm{s}^1`$ ($`\dot{M}_1=3\times 10^4,\text{6}\times 10^4,\mathrm{\hspace{0.33em}10}^3\mathrm{M}_{}\mathrm{yr}^1`$ for $`v_\mathrm{}_2=1500,\mathrm{\hspace{0.33em}3000},\mathrm{\hspace{0.33em}5000}\mathrm{km}\mathrm{s}^1`$). The mixed gas tends to be at these intermediate temperatures.

The overall net effect is that the emission from the hot gas is suppressed relative to the cooler gas, which softens the spectrum from the model with $`v_\mathrm{}_2=5000\mathrm{km}\mathrm{s}^1`$ relative to the model with $`v_\mathrm{}_2=3000\mathrm{km}\mathrm{s}^1`$. These inferences are supported by spectra calculated with $`v_\mathrm{}_2=5000\mathrm{km}\mathrm{s}^1`$ and $`\eta =0.1`$, and variable $`\dot{M}_2`$ (Fig. 4). This shows that the spectrum hardens as $`\dot{M}_2`$ increases from $`10^6\mathrm{M}_{}\mathrm{yr}^1`$ to $`10^4\mathrm{M}_{}\mathrm{yr}^1`$ (the spectral index $`\alpha `$, where $`F_\nu \nu ^\alpha `$, increases by about 0.1). This is what we expect from our prior inferences since a higher $`\dot{M}_2`$ helps to make postshock gas with a preshock velocity of $`v_\mathrm{}_2=5000\mathrm{km}\mathrm{s}^1`$ somewhat more efficient at radiating. In summary, we find that varying $`\dot{M}_2`$ can actually have a small effect on the shape of the spectrum in some parts of parameter space, as well as the more obvious and much greater effect on the luminosity. The power of $`\dot{M}_2`$ on affecting the spectral shape is, however, much less than that of $`v_\mathrm{}_2`$, which is the primary influence.

#### 2.2.4 Summary of Spectral Dependence on $`\eta `$, $`\dot{M}_2`$, and $`v_\mathrm{}_2`$

For the models investigated in this paper, we find that the spectral shape is i) insensitive to the value of $`\eta `$, ii) can harden slightly as $`\dot{M}_2`$ increases (in other parts of parameter space, e.g. strongly radiative shocks, it may soften with increasing $`\dot{M}_2`$), and iii) generally hardens with increasing $`v_\mathrm{}_2`$, but can sometimes soften. The luminosity increases with i) decreasing $`\eta `$ (since $`\eta `$ Car is unusual in the sense that there is no discernible contribution to the X-rays from the shocked primary wind, bar mixing), ii) increasing $`\dot{M}_2`$, and iii) increasing $`v_\mathrm{}_2`$.

### 2.3 The Spectral Grid

The synthetic spectra were input into a FITS-style ‘Table’ model for use with the XSPEC data-analysis package, as specified by Arnaud (A1999 (1999)). The fitting procedure within XSPEC is able to interpolate between the parameter values on the spectral grid so that best-fit values of $`\dot{M}_2`$ etc. are obtained from fitting the Chandra data. In the $`210`$ keV band, the intrinsic luminosity on our grid ranges from $`1.8\times 10^{31}\mathrm{erg}\mathrm{s}^1`$ to $`1.6\times 10^{37}\mathrm{erg}\mathrm{s}^1`$, which brackets the observationally determined value of $`4\times 10^{34}\mathrm{erg}\mathrm{s}^1`$ (Corcoran et al. 2001b ).

## 3 Fits to the Chandra grating spectrum

We compared the model intrinsic colliding wind spectra with the observed grating spectrum allowing the mass-loss rate and wind terminal velocity of the companion to vary. We in addition included absorption from a cool overlying medium to simulate the combined absorption of the circumstellar (nebular) material in the Homunculus, the ISM, and the wind from the companion (all of which may contribute to the X-ray absorption at some level). Since the stellar wind absorption depends on the companion’s mass-loss rate (at the phase of the Chandra observation) for consistency we should model this component of the absorption directly, but this is beyond the scope of the present paper. Because of the low number of counts per spectral bin in the high-resolution grating spectrum, and because background does not contribute significantly to the observed spectrum at energies above 2 keV, we fit the gross spectrum using the C-statistic (Cash C1979 (1979)), which is appropriate for Poisson-distributed data.

The fact that the shapes of our model spectra are not dependent on the value of $`\eta `$ introduces a degeneracy into our grid as both $`\eta `$ and $`\dot{M}_2`$ primarily influence the normalization of the spectra. This means that various combinations of $`\eta `$ and $`\dot{M}_2`$ can provide similar quality fits to the data. However, it is possible to obtain an estimate for $`\eta `$ from the duration of the lightcurve minimum ($`100`$ d). To do this we have constructed a simple model of the observed X-ray emission. We assume that the intrinsic flux varies as $`1/D`$, and that the absorption from the shock apex along the line of sight is negligible when viewing through the less dense secondary wind, but total when viewing through the very dense primary wind. We use the latest orbital parameters (Corcoran et al. 2001a ) and assume that the skew angle of the shock cone from the orbital velocity, $`v_{\mathrm{orb}}`$, is $`\delta \mathrm{arctan}(v_{\mathrm{orb}}/v_\mathrm{}_1)`$. Fig. 5 shows the results, scaled so that $`L_x=1`$ at periastron. The skew of the shock breaks the symmetry of the observed emission so that the post-minimum flux is lower than the pre-minimum flux, as observed. The duration of the minimum decreases with increasing $`\eta `$ (for $`\eta =0.1,0.2,0.3`$ the duration is $`133,92,71\mathrm{d}`$), and is best matched by $`\eta 0.2`$222Incidently, the lightcurve minimum computed with the 2D hydrodynamical model in Pittard et al. PSCI1998 (1998) (redisplayed in Fig. 1 at the beginning of this paper) also matches the duration of minimum very well (although not the post-minimum flux because the skew of the shock could not be computed with a 2D calculation). Although different mass-loss rates and a slightly less eccentric orbit were used, the wind momentum ratio adopted was also $`\eta =0.2`$, which provides a certain robustness for this value.. Therefore we adopt this value for the rest of our analysis.

We attempted to fit the grid of synthetic colliding wind spectra to the observed MEG +1 order spectrum. We also constrained the spectrum normalization to lie between 0.8 and 1.2, which fixes the distance to the star between 2300 and 1900 pc, close to the canonical distance of 2100 pc. In Fig. 6 we show the best fit interpolated spectrum from our grid in the range $`1.57.5`$ keV. Fig. 7 shows a close-up of the fit to the Si and S emission lines between $`1.82.8`$ keV. Both figures demonstrate the excellent agreement found between the data and the models and highlight the progress made in both the observational and theoretical study of colliding winds over the last 10 years. The continuum shape and level from the models matches closely that of the data. All strong lines (and most weak ones) which appear in the observed spectrum are matched also in the best fit. This means that the temperature distribution seen in the grating spectrum is matched by the model. However most of the strong lines are significantly overpredicted by the model fit. This discrepancy is in the opposite sense to the fit results from simple one- and two-temperature Raymond-Smith models, and is perhaps revealing that our assumption of solar abundances needs to be modified. The best fit parameters are summarized in Table 2.

## 4 Conclusions

In this paper our aim has been to test the binary hypothesis of $`\eta `$ Car by directly fitting its X-ray spectrum using a grid of spectra calculated from hydrodynamical models of the wind-wind collision. While our analysis does not prove that it is a binary, we find that the colliding wind emission model naturally provides for the range of ionization seen in the emission line grating spectrum for reasonable values of the wind parameters. We have not shown that it is inconsistent with emission from a single star, but have noted that it is unlike any of the other single stars observed so far at high energies and dispersion.

The technique applied in this paper has only been demonstrated once before and is the first time that it has been used with a high quality grating spectrum. Due to the low velocity of the primary wind ($`500\mathrm{km}\mathrm{s}^1`$), most of the observed X-ray emission arises from the shocked wind of the companion star. We find it difficult therefore to fit both $`\dot{M}_1`$ and $`v_\mathrm{}_1`$ as free parameters. However, the duration of the observed X-ray minimum can be used to estimate the wind momentum ratio of the stars, $`\eta `$. With $`\eta `$ fixed at 0.2, and $`\dot{M}_2`$, and $`v_\mathrm{}_2`$ as free parameters, we are able to obtain a good fit to the data.

We find that the mass-loss rate of the companion is $`\dot{M}_210^5\mathrm{M}_{}\mathrm{yr}^1`$ and the terminal velocity of its wind is $`v_\mathrm{}_23000\mathrm{km}\mathrm{s}^1`$. These values suggest that the companion is probably an Of supergiant (O-stars with similar wind parameters - e.g. HD 15570 (O4If), HD 93129A (O3If), HD 93250 (O3Vf), HD 151804 (O8If), and Cyg OB2 #9 (O5If) - are listed in Howarth & Prinja HP1989 (1989)), or is possibly a WR star. The velocity of the primary wind has been determined to lie in the range $`500700\mathrm{km}\mathrm{s}^1`$ (e.g. Hillier et al. HDIG2001 (2001)). Hence our fit implies a primary mass-loss rate of $`\dot{M}_12.5\times 10^4\mathrm{M}_{}\mathrm{yr}^1`$. From the uncertainty in the value of $`\eta `$, and the interpolation on our grid, we estimate the uncertainty on our derived value for $`\dot{M}_1`$ as approximately a factor of 2.

Our best-fit estimate of $`\dot{M}_1`$ is smaller than typically inferred (cf. Davidson & Humphreys DH1997 (1997); Hillier et al. HDIG2001 (2001)) However, we note that a lower mass-loss rate for the primary star can reduce some of the problems noted by Hillier et al. (HDIG2001 (2001)) who fitted a value as high as $`10^3\mathrm{M}_{}\mathrm{yr}^1`$. In particular, the models of Hillier et al. (HDIG2001 (2001)) suffered from absorption components that were too strong and electron-scattering wings which were overestimated. Both indicate that their chosen mass-loss rate is too high. Paradoxically both the H$`\alpha `$ and H$`\beta `$ emission lines were weaker than observed, indicating that their mass-loss rate is too low. However, it is well known that the wind collision zone can be a strong source of emission lines (e.g. HD 5980, Moffat et al. M1998 (1998)), which would resolve this problem. Their inferred minimum column density is also larger than the observed X-ray value, again implying an overestimate of $`\dot{M}_1`$. To increase the mass-loss rate of the primary towards the value estimated by Hillier et al. (HDIG2001 (2001)), we require either a reduction in the wind velocity to $`100\mathrm{km}\mathrm{s}^1`$, which in turn is in conflict with previous estimates, or a reduction in the wind momentum ratio to $`\eta <0.1`$, which is in conflict with the observed duration of the X-ray minimum. Finally, it is worth noting that our results indicate a value for $`\dot{M}_1`$ which is closer to the value inferred for the Pistol star ($`\dot{M}4\times 10^4\mathrm{M}_{}\mathrm{yr}^1`$; Figer et al. F1998 (1998)), an extreme early-type star with similarities to $`\eta `$ Car.

Since current observations and theoretical modelling of the optical spectrum are unable to determine the effective temperature and the stellar radius of the primary without first determining $`\dot{M}_1`$ (cf. Hillier et al. HDIG2001 (2001)), our independent estimate may prove to be extremely useful in this regard. It will be interesting to see if the results from this paper are consistent with future X-ray observations, and whether estimates of $`\dot{M}_1`$ from observations at X-ray and other wavelengths can be reconciled. Future X-ray grating observations should also help us to fix the value of the wind momentum ratio more accurately.

As the secondary wind dominates the X-ray spectrum, and its terminal velocity appears to be high ($`3000\mathrm{km}\mathrm{s}^1`$), we should expect to see signs of Doppler broadening and shifts in the line profiles. While there is little evidence for this in the current spectrum, other orbital phases may be more favourable in this regard. The addition of Doppler effects has already been incorporated in modelled X-ray spectra (Pittard et al. in preparation), and should provide further information on wind velocities and the structure of the wind-wind collision region.

The over-prediction of the X-ray lines in our models perhaps indicates that the companion has sub-solar abundances, which favours an O-type over a WR classification, although we would need to perform a more detailed analysis to confirm this possibility. As the primary has slightly enhanced abundances of C and N compared to solar, this suggests that to date there has been no mass exchange between the stars. Lamers et al. (LLPW1998 (1998)) suggested that the star which dominates the UV GHRS spectrum is not the star which ejected the nebula since the abundances in the GHRS spectrum are not as evolved as the abundances in the nebula (which are indicative of CNO-cycle products). As the UV bright source is probably the companion (Hillier et al. HDIG2001 (2001)) this indicates that it was the primary which ejected the nebula. There are also some caveats about the analysis in Lamers et al. (LLPW1998 (1998)) since strong C lines (which Lamers et al. took to indicate normal CNO abundances in the GHRS spectrum) also appear in stars known to be deficient in C (see the discussion in Hillier et al. HDIG2001 (2001)). It is also interesting to note that the UV spectrum brightened in 1999.1 vs. 1998.2, which is suggestive of an eclipse of the UV source near periastron in 1998. In conclusion, the high energy photons (UV and X-ray) seem to be telling us about the companion.

We emphasize that contrary to the vast majority of colliding wind systems, our X-ray analysis of $`\eta `$ Car primarily probes the conditions of the shocked wind of the companion. X-ray observations of $`\eta `$ Car are therefore unique in this regard since at other wavelengths (with the possible exception of the far UV) the wind of the primary dominates the observed phenomena. While our analysis has for the first time provided a direct estimate of the wind parameters of the companion star, relating these to the stellar parameters (mass, radius, luminosity) of the companion star requires more work. It is anticipated that the continued multi-wavelength study of $`\eta `$ Car through and beyond the next periastron passage will further reveal the hidden secrets of this most enigmatic system.

###### Acknowledgements.

We would like to thank Keith Arnaud for help with constructing a table model suitable for use with XSPEC. JMP would like to thank PPARC for the funding of a PDRA position, and MFC acknowledges that support for this research was obtained through a cooperative agreement with NASA/GSFC: NCC5-356. We would also like to thank the referee Kerstin Weis whose suggestions improved this paper. This work has made use of NASA’s Astrophysics Data System Abstract Service.


format: text

text:
# Anatomy of extreme events in a complex adaptive system

## Abstract

We provide an analytic, microscopic analysis of extreme events in an adaptive population comprising competing agents (e.g. species, cells, traders, data-packets). Such large changes tend to dictate the long-term dynamical behaviour of many real-world systems in both the natural and social sciences. Our results reveal a taxonomy of extreme events, and provide a microscopic understanding as to their build-up and likely duration.

PACS: 87.23.Ge, 02.50.Le, 05.45.Tp, 87.23.Kg

Large unexpected changes or ‘extreme events’ (e.g. crashes in financial markets, or punctuated equilibria in evolution) happen infrequently, yet tend to dictate the long-term dynamical behaviour of real-world systems in disciplines as diverse as biology and economics, through to ecology and evolution. The ability to generate large internal, so-called *endogenous* changes is a defining characteristic of complex systems, and arguably of Nature and Life itself . Such changes are manifestations of subtle, short-term temporal correlations resulting from internal collective behaviour. They seem to appear out of nowhere and have long-lasting consequences. To what extent can they ever be ‘understood’? Followers of the self-organized criticality view  would claim this question is naïve because of an inherent self-similarity in Nature: any large changes are simply magnified versions of smaller changes, which are in turn magnified versions of even smaller changes, and so on. Such self-similarity is presumed to underlie the power-law scaling observed in natural, social and economic phenomena . However, there are reasons for believing that the largest changes may be ‘special’ in a microscopic sense . Power-law scaling is only approximately true, and does not apply over an infinite range of scales. Apart from being atomistic at the smallest scale, a population of competing agents cannot cause any effect larger than the population size itself: in short, the largest changes will tend to ‘scrape the barrel’ in some way. Reference  quotes Bacon from Novum Organum: “Whoever knows the ways of Nature will more easily notice her deviations; and, on the other hand, whoever knows her deviations will more accurately describe her ways”.

This paper addresses the task of understanding, and eventually controlling, the large endogenous changes arising in a complex adaptive system comprising competing agents (e.g. species, cells, traders, data-packets). Our work reveals a taxonomy of large changes, and provides a quantitative microscopic description of their build-up and duration. Our results also provide insight into how a ‘complex systems manager’ might contain or control such extreme events.

We consider a generic complex system in which a population of $`N_{tot}`$ heterogeneous agents with limited capabilities and information, repeatedly compete for a limited global resource. Our model was introduced in Ref. , and is a generalization of the El Farol bar problem and the Minority game, concerning a population of people deciding whether to attend a popular bar with limited seating . At timestep $`t`$, each agent (e.g. a bar customer, or a market trader) decides whether to enter a game where the choices are option 1 (e.g. attend the bar, or buy) and option 0 (e.g. go home, or sell): $`N_0`$ agents choose 0 while $`N_1`$ choose 1. The ‘excess demand’ $`D[t]=N_1N_0`$ (which mimics price-change in a market) and number $`V[t]=N_1+N_0`$ of active agents (which mimics volume of market orders) represent output variables. These two quantities fluctuate with time, and can be combined to construct other global quantities of interest for the complex system studied (e.g. summing the price-changes gives the current price). This model can reproduce statistical and dynamical features similar to those in a real-world complex adaptive system, namely a financial market , and exhibits the crucial feature of seemingly spontaneous large changes of variable duration . The resulting time-series appears ‘random’ yet is non-Markovian, with subtle temporal correlations which put it beyond any random-walk based description. The temporal correlations of price-changes and volume, and their cross-correlation, are of intense interest in financial markets where so-called chartists offer a wide range of rules-of-thumb  such as ‘volume goes with price trend’. Although such rules are unreliable, the intriguing question remains as to whether there could *in principle* be a ‘science of charting’.

A subset $`V[t]N_{tot}`$ of the population, who are sufficiently confident of winning, are active at each timestep. For $`N_1<N_0`$ the winning decision is 1 and vice-versa, i.e. the winning decision is given by $`H[D[t]]`$ where $`H[x]`$ is the Heaviside function. The global resource level is so limited, or equivalently the game is so competitive, that at least half the active population lose at each timestep . The only global information available to the agents is a common bit-string ‘memory’ of the $`m`$ most recent outcomes. Consider $`m=2`$; the $`P=2^m=4`$ possible history bit-strings are 00, 01, 10 and 11, which can also be represented in decimal form: $`\mu \{0,1,\mathrm{},P1\}`$. A strategy consists of a response, $`a^\mu \{1,1\}`$ to each possible bit-string $`\mu `$, $`a^\mu =1`$ option 1, and $`a^\mu =1`$ option 0. Hence there are $`2^P=16`$ possible strategies. The heterogeneous agents randomly pick $`s`$ strategies each at the outset, and update the scores of their strategies after each timestep with the reward function $`\chi [D]=sgn[D]`$, i.e. $`+1`$ for choosing the minority action, $`1`$ for choosing the majority action. Agents have a time horizon $`T`$ over which strategy points are collected, and a threshold level $`r`$ which mimics a ‘confidence’. Only strategies having $`r`$ points are used, with agents playing their highest scoring strategy. Agents with no such strategy become temporarily inactive . We focus on the regime where the number of strategies in play is comparable to the total number available, since this yields seemingly random dynamics with occasional large movements . The coin-tosses used to resolve ties in decisions (i.e. $`N_0=N_1`$) and active-strategy scores, inject stochasticity into the game’s evolution. Reference  showed that a simplified version of this system in the limit $`r\mathrm{}`$ and $`T\mathrm{}`$, can be usefully described as a stochastically disturbed deterministic system. We are interested in the dynamics of large changes, and adopt the approach and terminology of Ref. . Averaging over our model’s stochasticity yields a description of the game’s deterministic dynamics via mapping equations for the strategy score vector $`\underset{¯}{S}[t]`$ and global information $`\mu [t]`$. For $`s=2`$ the deterministic dynamics are given exactly by the following equations:

$`\underset{¯}{S}[t]`$ $`=\underset{¯}{S}[0]{\displaystyle \underset{i=tT}{\overset{t1}{}}}\underset{¯}{a}^{\mu [i]}sgn\left[D[i]\right],`$ (1)
$`\mu [t]`$ $`=2\mu [t1]PH\left[\mu [t1]P/2\right]+H\left[D[t1]\right].`$
The corresponding demand function is given by

$$D[t]=\underset{R=1}{\overset{2P}{}}a_R^{\mu [t]}H[S_Rr]\underset{R^{}=1}{\overset{2P}{}}\left(1+sgn\left[S_R[t]S_R^{}[t]\right]\right)\mathrm{\Psi }_{R,R^{}},$$

where $`\underset{¯}{\underset{¯}{\mathrm{\Psi }}}`$ is the symmetrized strategy allocation matrix which constitutes the *quenched disorder* present during the system’s evolution . Elements $`\mathrm{\Psi }_{R,R^{}}`$ enumerate the number of agents holding both strategy $`R`$ and $`R^{}`$. The volume $`V[t]`$ is given by the same expression as $`D[t]`$ replacing $`a_R^{\mu [t]}`$ by unity.

Large changes such as financial market crashes, seem to exhibit a wide range of possible durations and magnitudes making them difficult to capture using traditional statistical techniques based on one or two-point probability distributions . A common feature, however, is an obvious trend (i.e. to the eye) in one direction over a reasonably short time window: we use this as a working definition of a large change. In fact, all the large changes discussed here represent $`>3\sigma `$ events. In both our model and the real-world system, these large changes arise more frequently than would be expected from a random-walk model . Our model’s dynamics can be described by trajectories on a de Bruijn graph : see Fig. 1 for $`m=3`$, with a transition incurring an increment to the score vector $`\underset{¯}{S}`$. There are $`P`$ orthogonal increment vectors $`\underset{¯}{a}^\mu `$, one for each node $`\mu `$. Setting the initial scores $`\underset{¯}{S}[0]=\underset{¯}{0}`$, the strategy score vector in Eq. (1) can be written exactly as:

$$\underset{¯}{S}[t]=c_0\underset{¯}{a}^0+c_1\underset{¯}{a}^1+\mathrm{}+c_{P1}\underset{¯}{a}^{P1}=\underset{j=0}{\overset{P1}{}}c_j\underset{¯}{a}^j$$

where $`c_j`$ represents the *nodal weights* for history node $`\mu =j`$. The nodal weights enumerate the number of negative return transitions from node $`\mu `$ minus the number of positive return transitions, in the time window $`tTt1`$. High absolute nodal weight implies persistence in transitions from that node i.e. persistence in $`D|\mu `$. Large changes will occur when connected nodes become persistent. The simplest type of large movement exhibiting perfect nodal persistence would be $`\mu =0,0,0,0,\mathrm{}`$ in which all successive price changes are in the *same* direction. We call this a ‘fixed-node crash’ (or rally). However, there are many other possibilities reflecting the wide range of forms and durations of the large change. For example, on the $`m=3`$ de Bruijn graph in Fig. 1 the cycle $`\mu =0,0,1,2,4,0,\mathrm{}`$ has four out of the five transitions producing price-changes of the same sign (it is persistent on nodes 1, 2, 4 and antipersistent on node 0). We call this a ‘cyclic-node crash’ (or rally). Figure 2 illustrates a large change which starts as a fixed-node crash then subsequently becomes a cyclic-node crash. Cyclic-node crashes can be treated simply as interlocking fixed-node crashes, hence for clarity we focus here on a single fixed-node crash (or rally). For the parameter ranges of interest, the choice about whether a strategy is played by an agent is more determined by whether that strategy’s score is above the threshold, than whether it is their highest-scoring strategy . This is because agents are only likely to have at most one strategy whose score lies above the threshold for confidence levels $`r0`$. Making the additional numerically-justified approximation of small quenched disorder (i.e. the variance of the entries in the strategy allocation matrix $`\underset{¯}{\underset{¯}{\mathrm{\Psi }}}`$ is smaller than their mean for the parameter range of interest ), the demand and volume become:

$`D[t]`$ $`={\displaystyle \frac{N}{4P}}{\displaystyle \underset{R=1}{\overset{2P}{}}}a_R^{\mu [t]}sgn\left[S_R[t]r\right],`$ (2)
$`V[t]`$ $`={\displaystyle \frac{N}{2}}+{\displaystyle \frac{N}{4P}}{\displaystyle \underset{R=1}{\overset{2P}{}}}sgn\left[S_R[t]r\right].`$ (3)
Suppose persistence on node $`\mu =0`$ starts at time $`t_0`$. How long will the resulting crash last? To answer this, we decompose Eq. (2) into strategies which predict 1 at $`\mu =0`$, and those that predict 0. We first consider the particular case where the node $`\mu =0`$ was *not* visited during the previous $`T`$ timesteps, hence the loss of score increment from time-step $`tT`$ will not affect $`\underset{¯}{S}[t]`$ on average. At any later time $`t_0+\tau `$ during the crash, (i.e. $`\mu =0`$) Eqs. (2) and (3) are hence given by:

$`D[t_0+\tau ]`$ $`={\displaystyle \frac{N}{4P}}\left\{{\displaystyle \underset{Ra_R^\mu =1}{}}sgn\left[S_R[t_0]r\tau \right]{\displaystyle \underset{Ra_R^\mu =1}{}}sgn\left[S_R[t_0]r+\tau \right]\right\},`$ (4)
$`V[t_0+\tau ]`$ $`={\displaystyle \frac{N}{2}}+{\displaystyle \frac{N}{4P}}\left\{{\displaystyle \underset{Ra_R^\mu =1}{}}sgn\left[S_R[t_0]r\tau \right]+{\displaystyle \underset{Ra_R^\mu =1}{}}sgn\left[S_R[t_0]r+\tau \right]\right\}.`$
$`|D[t_0+\tau ]|`$ decreases as the persistence time $`\tau `$ increases, and hence the crash ends at time $`t_0+\tau _c`$ when the right-hand side of Eq. (4) changes sign. The persistence time or ‘crash-length’ $`\tau _c`$ is thus given by the mean of the scores of the strategies predicting 0, i.e. $`\tau _c=\overline{S}_{Ra_R^\mu =1}[t_0]=c_0[t_0]`$. In the more general case where the node $`\mu =0`$ *was* visited during the previous $`T`$ timesteps, $`\tau _c`$ is given by the largest $`\tau `$ value which satisfies:

$$\tau =\left(c_0[t_0]+\underset{\{t^{}\}}{}sgn\left[D[t^{}]\right]\right)$$

where $`\{t^{}\}(\mu [t^{}]=0t_0Tt^{}t_0+\tau T)`$. Assume that the scores have a near-Normal distribution, i.e. $`S_{Ra_R^\mu =1}[t_0]N[\overline{S}_1,\sigma ]`$ as in Fig. 3a. For each strategy $`R`$ there exists an anticorrelated strategy $`\overline{R}`$ and hence $`S_R[t]=S_{\overline{R}}[t]`$ for all $`t`$. Consequently, prior to a crash, the score distribution tends to split into two halves as indicated schematically in Fig. 3a. The expected demand (and volume) during the crash are then:

$`<D[t_0+\tau ]>`$ $`\left(\mathrm{erf}\left[{\displaystyle \frac{c_0[t_0]+r+\tau }{\sqrt{2}\sigma }}\right]\mathrm{erf}\left[{\displaystyle \frac{c_0[t_0]+r\tau }{\sqrt{2}\sigma }}\right]\right)`$
$`<V[t_0+\tau ]>`$ $`\left(2\mathrm{erf}\left[{\displaystyle \frac{c_0[t_0]+r+\tau }{\sqrt{2}\sigma }}\right]\mathrm{erf}\left[{\displaystyle \frac{c_0[t_0]+r\tau }{\sqrt{2}\sigma }}\right]\right)`$
These forms are illustrated in Figure 3b. As the spread in the strategy score distribution is increased, the dependence of $`<D>`$ and $`<V>`$ on the parameters $`\tau `$ and $`r`$ becomes weaker and the surfaces flatten out leading to a smoother drawdown, as opposed to a sudden severe crash. As the parameters $`\overline{S}_1,\sigma ,r`$ are varied, it can be seen that the behaviour of the demand and volume during the crash can exhibit markedly different qualitative forms , yielding a *taxonomy* of *different species* of large change *even within the same single-node family*. This result could explain why financial market chartists’ rules-of-thumb , such as ‘volume goes with price trend’, are far too simplistic.

We now turn to the important practical question of whether history will repeat itself, i.e. given that a crash has recently happened, is it likely to happen again? If so, is it likely to be even bigger? Suppose the system has built up a negative nodal weight for $`\mu =0`$ at some point in the game (see Fig. 4a). It then hits node $`\mu =0`$ at time $`t_0`$ producing a crash (Fig. 4b). The nodal weight $`c_0`$ is hence restored to zero (Fig. 4c). In this model the previous build-up is then forgotten because of the finite $`T`$ score window, hence $`c_0`$ becomes positive (Fig. 4d). The system then corrects this imbalance (Fig. 4e), restoring $`c_0`$ to 0. The crash is then forgotten, hence $`c_0`$ becomes negative (Fig. 4f). The system should therefore crash again - however, a crash will *only* re-appear if the system’s trajectory subsequently returns to node $`\mu =0`$. Interestingly, we find that the *disorder* in the initial distribution of strategies among agents (i.e. the quenched disorder in $`\underset{¯}{\underset{¯}{\mathrm{\Psi }}}`$) can play a deciding role in the issue of crash ‘births and revivals’ since it leads to a slight bias in the outcome, and hence the subsequent transition, at each node. When $`c_{\mu [t]}=0`$ (see Fig. 4c), it follows that $`sgn[D[t]]`$ is more likely to be equal to $`sgn\left[\underset{¯}{a}^{\mu \left[t\right]}\underset{¯}{x}\right]`$ where $`\underset{¯}{x}=_{R^{}=1}^{2P}\underset{¯}{\underset{¯}{\mathrm{\Psi }}}_R^{}`$ is a strategy weight vector with $`x_R`$ corresponding to the number of agents who hold strategy $`R`$ . The quenched disorder therefore provides a crucial bias for determining the future trajectory on the de Bruijn graph when the nodal weight is small, and hence can decide whether a given crash recurs or simply disappears. The quenched disorder also provides a *catalyst* for building up a very large crash.

Our work opens up the study of how a ‘complex-systems-manager’ might use this information to control the long-term evolution of a complex system by introducing, or manipulating, such large changes. As an example, we give a quick three-step solution to prevent large changes: (1) use the past history of outcomes to build up an estimate of the score vector $`\underset{¯}{S}[t]`$ and the nodal weights $`\{c_{\mu [t]}\}`$ on the various critical nodes, such as $`\mu =0`$ in the case of the fixed-node crash. (2) Monitor these weights to check for any large build-up. (3) If such a build-up occurs, step in to prevent the system hitting that node until the weights have decreased.

Figure Captions

Figure 1: Dynamical behaviour of the global information is described by transitions on the de Bruijn graph. Graph for population of $`m=3`$ agents. Blue transitions represent positive demand $`D`$, red transitions represent negative demand.

Figure 2: Dynamical behaviour of complex system (e.g. price $`P\left[t\right]`$ in financial market) described by evolution of nodal weights $`c_\mu `$. History at each timestep indicated by black square. Large change preceded by abnormally high nodal weight. Large change incorporates fixed-node and cyclic node crashes

Figure 3: (a) Schematic representation of strategy score distribution prior to crash. Arrows indicate subsequent motion during crash period. (b) Plots of expected demand and volume during crash period showing range of different possible behaviour as system parameters are varied.

Figure 4: Representation of how large changes can recur due to finite memory of agents. Grey area shows history period outside agents’ memory. Example shows recurring fixed-node crash at node $`\mu =0`$.


text:
# Interactions between Pulsars, Pulsar Nebulae and Supernova Remnants

## 1. Introduction

Pulsars are thought to be born in the core collapse and explosion of massive stars, so that the initial evolution and observability of the PWNe (pulsar wind nebulae) depends on the properties of the supernova and the supernova remnant evolution. Several phases of evolution are expected (Pacini & Salvati 1973; Reynolds & Chevalier 1984; Chevalier 1998). Initially the radiation from the PWN is absorbed by the surrounding supernova gas; high energy photons may be reprocessed and appear at a lower energy. After $`>`$10’s of years, the supernova becomes optically thin and the PWN radiation directly escapes. The Crab Nebula may be in this phase of evolution. After $`10^3`$ years for a pulsar like that in the Crab, the spindown power from the pulsar decreases and the PWN fades. After $`10^4`$ years, the reverse shock wave from the supernova remnant moves to the center of the remnant, crushing the PWN and leading to an increase in brightness of the nebula.

The discussion of PWNe has typically centered on objects like the Crab and its pulsar, which is estimated to have an initial spin period $`P_019`$ ms and a magnetic field of $`5\times 10^{12}`$ G. It is thought to belong to the main (not binary related) group of radio pulsars, which appear to have a fairly narrow distribution of initial magnetic field strengths $`3\times 10^{12}`$ G (Stollman 1987). The initial rotation rates of these pulsars remain uncertain; studies of the pulsar population are influenced by selection effects (e.g., Emmering & Chevalier 1989). In recent years there has been increasing evidence for a class of strongly magnetized neutron stars, or magnetars, with $`B>10^{14}`$ G (Thompson & Duncan 1995). Magnetars are thought to be associated with soft gamma-ray repeaters and possibly with AXPs (anomalous X-ray pulsars), both of which appear to contain slowly rotating neutron stars. Although the youth of AXPs is indicated by their association with supernova remnants, there is no evidence for PWNe. Finally there is a class of “quiet” neutron stars that have been found as weak X-ray sources in young supernova remnants, such as Cas A (Tananbaum 1999) and Puppis A (Petre, Becker, & Winkler 1996). Although the relative formation rates of the normal radio pulsars, magnetars, and quiet neutron stars is not known, the likely youth of the objects in the latter two categories indicates that all three groups make a significant contribution to the total rate of massive star core collapses.

Here, I discuss four topics related to neutron stars and their supernova remnants. The first is the quest for young PWNe, considerably younger than the Crab Nebula. Next, I consider the status of the Crab Nebula and its surroundings. The third topic concerns the possible relation between a neutron star and its surrounding supernova, given the growing number of such associations. Finally, I review recent investigations of the effect of a reverse shock front on a PWN.

## 2. Searching for Young Pulsar Nebulae

The Crab pulsar, with an age of 947 years, was until recently the youngest known pulsar. A newly discovered pulsar in the remnant Kes 75 has a characteristic age of only 723 years (Gotthelf et al. 2000) and so may be younger. However, the search for very young pulsars (age $`<100`$ years) continues. The discovery of such objects is important for understanding the early evolution of pulsars. In their early phases, pulsars are expected to interact with their surrounding supernovae. This occurs both through the dynamical interaction of the pulsar bubble with the supernova and the interaction of the radiation with the surrounding supernova gas. If the pulsar injects energy at a constant rate $`\dot{E}`$ in a uniformly expanding central region of approximately constant density $`\rho _c`$, the bubble radius increases as (Chevalier 1977)

$$R=\left(\frac{125}{132\pi }\frac{\dot{E}}{\rho _ct^3}\right)^{1/5}t^{6/5}$$
(1)
where $`\rho _ct^3`$ is a constant because of the uniform expansion. The shock front generated by the pulsar bubble is expected to initially be radiative, with a total power $`0.015\dot{E}`$. The initial spin-down power of a pulsar in the vacuum dipole approximation can be expressed as

$$\dot{E}_0=1.6\times 10^{39}\left(\frac{B_p}{3\times 10^{12}\mathrm{G}}\right)^2\left(\frac{R_{ns}}{12\mathrm{km}}\right)^6\left(\frac{P_0}{20\mathrm{ms}}\right)^4\mathrm{erg}\mathrm{s}^1,$$
(2)
where $`B_p`$ is the dipole magnetic field at the pole of the star, $`R_{ns}`$ is the neutron star radius, and $`P_0`$ is the initial period. A dipole magnetic field perpendicular to the rotation axis has been assumed. The velocity of the shell at the bubble boundary is given by

$$V=470\left(\frac{\dot{E}}{10^{39}\mathrm{erg}\mathrm{s}^1}\right)^{1/5}\left(\frac{\rho _ct^3}{10^9\mathrm{g}\mathrm{s}^3\mathrm{cm}^3}\right)^{1/5}\left(\frac{t}{10\mathrm{yr}}\right)^{1/5}\mathrm{km}\mathrm{s}^1,$$
(3)
where the reference value of $`\rho _ct^3`$ is based on models for SN 1987A.

The radiation can be supplemented by synchrotron radiation from the pulsar bubble. The radiation from a pulsar bubble like the Crab Nebula appears to be from relativistic electrons and positrons that have been accelerated in a shock front which begins the deceleration of the pulsar wind (Rees & Gunn 1974; Kennel & Coroniti 1984a, 1984b). Kennel & Coroniti (1984b) find that the magnetization parameter $`\sigma `$ of the pulsar wind (ratio of Poynting flux to kinetic energy flux) must be $`0.003`$ to fit the Crab parameters. However, when $`\sigma <v/c`$ flow solutions with a shock front are no longer possible; the shock wave moves in toward the position of the pulsar (Kennel & Coroniti 1984a; Emmering & Chevalier 1987). In the case of the Crab, $`v2000\mathrm{km}\mathrm{s}^1`$, so that $`\sigma <6.7\times 10^3`$ is needed. The value indicated for the Crab is only slightly smaller than the limit. According to eqn. (1) with $`Vt^{1/5}`$, the current value of $`\sigma `$ was below the limit at an age of $`20`$ years, which shows that a young pulsar nebula may have different properties than the Crab.

As $`\sigma `$ approaches $`v/c`$ in nebular models, the fact that the shock front moves in to the central pulsar suggests that the rotating pulsar is in direct causal contact with the surrounding medium (Emmering & Chevalier 1987; Lyutikov 2001). There is no supersonic wind zone. The pulsar spindown law is then not given by eqn. (2), but is determined by the torque acting on the pulsar from the surrounding supernova gas through the magnetic field. The pulsar may do work on the surrounding medium, leading to some radiation, but the normal nebular synchrotron radiation does not occur.

If a pulsar wind has a low magnetization parameter or the pulsar is powerful so it can drive a fast surrounding shell, it is possible that even a young pulsar can have a surrounding pulsar nebula with synchrotron emission from shock accelerated particles. Chevalier & Fransson (1992) investigated the emission from a supernova that contains such a pulsar nebula. During the first 10 years, the emission is absorbed by the core gas around the pulsar nebula and re-radiated at optical and ultraviolet wavelengths. The line widths depend on the pulsar power, but a typical width is $`1000\mathrm{km}\mathrm{s}^1`$. If the core material is oxygen rich, the expected strong lines are \[O III\] $`\lambda \lambda `$4959,5007, O I $`\lambda `$7774, Mg I\] $`\lambda `$4571, and various ultraviolet O lines. A characteristic property of the line widths is that they should broaden with time because as the supernova density decreases the ionized region becomes broader.

The best tests for the presence of a young pulsar are in the nearest core collapse supernovae, SN 1987A in the LMC (Large Magellanic Cloud) and SN 1993J in M81. The ultraviolet/optical spectrum of SN 1987A is well explained by power input from $`10^4M_{}`$ of 44Ti; any other power source is at a level $`<10^{36}\mathrm{erg}\mathrm{s}^1`$ (Chugai et al. 1997). If a quiet neutron star is present in SN 1987A, accretion onto the neutron star is expected, giving an Eddington limit luminosity. An accretion rate of at least the electron scattering Eddington rate is expected for $`200`$ years (Chevalier 1989), but the observed luminosity limit is considerably lower than the expected $`2\times 10^{38}\mathrm{erg}\mathrm{s}^1`$. One possibility is that a high Fe opacity effectively decreases the Eddington limit; a photon driven wind can then clear the region around a central neutron star on a timescale of $`1`$ yr (Fryer, Colgate, & Pinto 1999). Better limits on a central source of emission will become available because SN 1987A can now be spatially resolved at optical (Pun et al. 2001), X-ray (Park et al. 2001a), and radio (Manchester et al. 2001) wavelengths; in no case is there any sign of central emission.

The case of SN 1993J is complicated by the fact that it is undergoing strong circumstellar interaction (Matheson et al. 2000). The lines that are mentioned above as resulting from a pulsar nebula can also be produced by circumstellar interaction and are observed in SN 1993J. Methods of distinguishing between the sources of power are that the circumstellar interaction lines tend to be initially broader and to narrow with time because of deceleration, whereas the pulsar nebula lines are initially narrower and broaden with time, as discussed above.

While extragalactic supernovae are excellent sites to search for very young pulsar nebulae, somewhat older ones might be found by searching for extragalactic objects like the Crab Nebula. In the LMC, the remnant 0540-69 has a pulsar nebula very much like the Crab and may the penultimate core collapse supernova in the galaxy before SN 1987A (Kirshner et al. 1989). Beyond the LMC, however, pulsar nebulae have not been identified. At radio wavelengths, Reynolds & Fix (1987) failed to find PWNe in the galaxy M33. At X-ray wavelengths, deep surveys of star forming galaxies are becoming available and it is likely that PWNe are being detected. For example, Pence et al. (2001) have detected 110 sources in an image of M101 down to a level of $`10^{36}\mathrm{erg}\mathrm{s}^1`$ in the 0.5–10 keV band. The problem is to distinguish the PWNe from binary and other sources. A strongly absorbed spectrum is possible if absorption by the surrounding supernova gas is still important.

## 3. The Crab and its Surroundings

A long-standing mystery concerning the Crab Nebula is the nature of the medium into which it is expanding and which contains the pulsar nebula. There has not yet been an observation of an extended interaction region around the Crab. One possibility is that it is expanding into the freely expanding ejecta of SN 1054 (e.g., Chevalier 1977). Another is that the surroundings are a relatively slowly moving stellar wind from mass loss before the supernova (e.g., Nomoto et al. 1982).

There are a number of arguments that appear to favor the ejecta as the surrounding medium. In the ejecta case, the outer edge of the Crab is expected to be slowly accelerating or at least moving at a constant velocity. The observations of proper motion are consistent with free expansion, even in the extended “jet” on the north side of the Crab (Fesen & Staker 1993). If the Crab were interacting with a slow wind, deceleration at the outer edge would be expected.

In the ejecta interaction model, a radiative shock wave with velocity of $`230\mathrm{km}\mathrm{s}^1`$ is driven into ejecta over much of the lifetime of the Crab (Chevalier & Fransson 1992). The shocked ejecta cool and accumulate in the observed optical filaments. The filamentary structure can be attributed to the action of the Rayleigh-Taylor instability. Jun (1998) has shown in numerical simulations how the instability can give rise to structures like those observed in the Crab. Sankrit & Hester (1997) have found observational evidence that the radiative shock wave is still present in the Crab, based on the strength of \[O III\] lines at the boundary of the Crab. These lines are difficult to form by photoionization and are more likely to be due to a shock front. In their picture, the shock front is radiative only over the middle part of the Crab, so it is currently making a transition from the radiative to the nonradiative phase.

A recent observation is the ultraviolet spectrum of the Crab pulsar in a search for absorption by freely expanding ejecta gas (Sollerman et al. 2000). An absorption feature in the C IV $`\lambda `$1550 line was found with a negative velocity extent to $`2500\mathrm{km}\mathrm{s}^1`$. This velocity is somewhat beyond the velocity of the filaments in the Crab Nebula, but not by much. It does not extend to the high velocities that would be needed if an explosion energy of $`10^{51}`$ ergs were in the freely expanding ejecta with a plausible density profile. One of the initial arguments for expansion into free ejecta was that the ejecta could be carrying a normal supernova energy, as opposed to the low energy that is directly observed in the Crab Nebula (Chevalier 1977).

The overall result then is that the case for interaction with freely expanding ejecta is good, but there is no evidence for a normal supernova energy in the ejecta. The lower limit found by Sollerman et al. (2000) is only $`1.5\times 10^{49}`$ ergs. H$`\alpha `$ emission should provide a means of detecting more extended, cool ejecta, but only upper limits have been obtained so far (Fesen, Shull, & Hurford 1997). Type II supernovae with low expansion velocities have been found, such as SN 1994W (Sollerman, Cumming, & Lundqvist 1998) and SN 1997D (Turatto et al. 1998), but it is not clear whether these supernovae have unusually low energies.

## 4. Neutron Stars and Their Supernovae

The increasing number of neutron stars discovered in association with young supernova remnants gives the possibility of searching for correlations between core collapse supernovae and their central compact objects. Table 1 contains a list of core collapse supernovae with ages $`<3,000`$ years. The ages of 3C 58, G11.2–0.3, and MSH 15–52 have been put in parentheses because uncertain identifications with SN 1181, SN 386, and SN 185, respectively, have been assumed. The first 10 objects are in the Galaxy, the next 3 are in the LMC, and the last one is in the SMC. If there is not direct evidence for a compact object, the determination of a supernova being of the core collapse variety is made based on element overabundances in the ejecta or in the circumstellar medium (CSM). A controversial object is Kepler’s remnant, which has been suggested to have a massive star origin based on the circumstellar material (Bandiera 1987). Gerardy & Fesen (2001) note the similarity between the spectrum of a circumstellar knot in Kepler and the quasi-stationary flocculi in Cas A.

The nature of the compact source and the X-ray luminosity of the PWN (column 3 and 4) are related in that normal pulsars appear to be capable of creating synchrotron wind nebulae. The quoted X-ray luminosity, $`L_x`$, is in the 0.2–4 keV range. The pulsar in 3C 58 was discovered only recently (Murray et al. 2001). The compact source in G292.0+1.8 is likely to be a normal pulsar, given the other correspondences with PWNe.

The weak X-ray sources in Cas A, RCW 103, and Puppis A appear to be qualitatively different from the pulsars. In the case of Puppis A, Gaensler et al. (2001) have set strong limits on the presence of a radio pulsar nebula. The other 4 objects have neither an observed compact object nor a PWN. Three of these are sufficiently distant that a Cas A type source or a weak pulsar would not be observable. Kepler’s remnant is closer and is worthy of a careful search for a compact object.

The overall impression from Table 1 is that remnants with normal pulsars and nebulae tend not to have O-rich ejecta, whereas the O-rich young remnants tend not to have pulsars. One possible exception is 0540–69 in the LMC. This remnant has weak H emission, but lines of O, S, Ar, Ni, and Fe have relative strengths that are comparable to those observed in the Crab Nebula (Kirshner et al. 1989). This remnant is worthy of detailed abundance studies. A clearer exception is G292.0+1.8.

The implication is that the pulsars are born in stars with a range of masses. The ejecta in the Crab Nebula are not O-rich and are consistent with a low mass star, $`810M_{}`$ (Nomoto et al. 1982). The presence of O-rich core material in G292.0+1.8 implies an initial mass $`>13M_{}`$. The tendency of the pulsars to avoid O-rich remnants suggests that they typically have lower mass progenitors, but another possibility is that the O-rich remnants are ones with a dense circumstellar medium, so that the ejecta are observable through their interaction with the surroundings. Blair et al. (2000) suggested that N132D and 1E 0102 had Wolf-Rayet star progenitors, and this suggestion has long been made for Cas A (Chevalier 1976). However, G292.0+1.8 is a case with a pulsar where strong circumstellar interaction appears to be taking place (Park et al. 2001b). There is thus not a clear link between neutron stars and their supernovae, but the data are reaching the point where interesting questions can be asked.

## 5. The Reverse Shock Wave Interaction

In the later stage of evolution of a supernova remnant, the reverse shock front in the ejecta moves back to the center of the remnant, crushing a central PWN. If the surrounding density is $`1`$ H atom cm-3, this occurs at an age $`10^4`$ years. Reynolds & Chevalier (1984) discussed the effects of the crushing on increasing the synchrotron luminosity of a PWN. There have recently been hydrodynamic studies of the crushing process.

Van der Swaluw et al. (2001) et al. carried out one dimensional simulations of the process using a one fluid code. The expansion of the PWN inside the supernova was approximately as described by the analytical model of Reynolds & Chevalier (1984) during the pre- and post-reverse shock phases, but there were considerable transient waves during the transition. During the late phase, the radius of the PWN increases as $`t^{0.3}`$, which is close to the expansion rate of an adiabatic blast wave $`t^{0.4}`$. The ratio of the radius of a PWN to the radius of the outer shock front can thus be related to the energy injected by the pulsar, with only a weak dependence on age. Van der Swaluw & Wu (2001) use this property to obtain estimates of the total energy deposited into the PWN nebula for a number of observed cases and thus deduce the initial periods of the pulsars. The result is that the deduced periods are $`>10`$’s of seconds, i.e. relatively long. However, radiative losses from the PWN can allow a larger initial energy. In addition, a large energy deposited early in the evolution is primarily lost to adiabatic expansion, so that this method is not very sensitive to the initial period (Blondin, Chevalier, & Frierson 2001).

Blondin et al. (2001) carried out two dimensional, two fluid simulations of the late evolution of a PWN and its surrounding supernova remnant. The PWN was treated as a $`\gamma =4/3`$ gas with a high sound speed. These simulations showed that the crushing and re-expansion of the PWN are subject to hydrodynamic instabilities so that the thermal and nonthermal gases become mixed. The simulations did not include a magnetic field; the formation of magnetic filaments is likely to occur. In addition, the reverse shock front is unlikely to be spherically symmetric, so that the PWN can be pushed to the side relative to its parent pulsar. This effect is most likely to be observed at radio wavelengths because of the large ages of the radio emitting electrons, and can explain the off-center position of the Vela pulsar relative to the Vela X nebula. Off center radio pulsar nebulae have also been found in a number of other sources.

## 6. Discussion

Of the various phases expected for pulsar nebula evolution, the earliest phase when nebular emission is absorbed by the surrounding supernova is least understood because of a lack of observations. The next phase when the PWN can be observed directly and is expanding into freely expanding supernova ejecta is presumably observed in many cases, although clear proof of this phase has been difficult to obtain even in the case of the Crab Nebula. Most of the remnants listed in Table 1, with ages $`<3000`$ yr, are likely to be in this phase. Continued study of these sources, especially analysis of the composition from X-ray emission, is warranted in order to search for correlations between neutron stars and their progenitor stars. In the late, post-reverse shock phase, a pulsar nebula is easily distorted and displaced by the thermal gas. Careful studies of the boundaries of these nebulae should reveal the nature of the interaction.

### Acknowledgments.

This work was supported in part by NASA grant NAG5-8130.

## References

Bandiera, R. 1987, ApJ, 319, 885

Blair, W. P., et al. 2000, ApJ, 537, 667

Blanton, E. L., & Helfand, D. J. 1996, ApJ, 470, 961

Blondin, J. M., Chevalier, R. A., & Frierson, D. M. 2001, ApJ, 563, 806

Canizares, C. R., & Winkler, P. F. 1981, ApJ, 246, L33

Chevalier, R. A. 1976, ApJ, 208, 826

Chevalier, R. A. 1977, in Supernovae, ed. D. N. Schramm (Dordrecht: Reidel), 53

Chevalier, R. A. 1989, ApJ, 346, 847

Chevalier, R. A. 1998, Mem. Soc. Astron. Italia, 69, 977

Chevalier, R. A., & Fransson, C. 1992, ApJ, 511, 798

Chevalier, R. A., & Kirshner, R. P. 1978, ApJ, 219, 931

Chugai, N. N., Chevalier, R. A., Kirshner, R. P., & Challis, P. M. 1997, ApJ, 483, 925

Davidson, K., & Fesen, R. A. 1985, ARA&A, 23, 119

Emmering, R. T., & Chevalier, R. A. 1987, ApJ, 321, 334

Emmering, R. T., & Chevalier, R. A. 1989, ApJ, 345, 931

Fesen, R. A., & Staker, B. 1993, MNRAS, 263, 69

Fesen, R. A., Kirshner, R. P., & Becker, R. H. 1988, in Supernova Remnants and the Interstellar Medium, ed. R. S. Roger & T. L. Landecker (Cambridge: CUP), 55

Fesen, R. A., Shull, J. M., & Hurford, A. P. 1997, AJ, 113, 354

Fryer, C. L., Colgate, S. A., & Pinto, P. A. 1999, ApJ, 511, 885

Gaensler, B. M., Bock, D. C.-J., & Stappers, B. W. 2001, ApJ, 537, L35

Gerardy, C. L., & Fesen, R. A. 2001, 121, 2781

Gotthelf, E. V., Vasisht, G., Boylan-Kolchin, M., & Torii, K. 2000, ApJ, 542, L37

Hughes, J. P., Rakowski, C. E., Burrows, D. N., & Slane, P. O. 2000, ApJ, 528, L109

Hughes, J. P., et al. 2001, ApJ, 559, L153

Jun, B.-I. 1998, ApJ, 499, 282

Kennel, C. F., & Coroniti, F. V. 1984a, ApJ, 283, 694

Kennel, C. F., & Coroniti, F. V. 1984b, ApJ, 283, 710

Kirshner, R. P., Morse, J. A., Winkler, P. F., & Blair, W. P. 1989, ApJ, 342, 260

Kozma, C., & Fransson, C. 1998, 497, 431

Leibowitz, E. M., & Danziger, I. J. 1983, MNRAS, 204, 273

Lyutikov, M. 2001, talk at this meeting

Manchester, R. N., et al. 2001, PASA, in press (astro-ph/0110693)

Matheson, T., Filippenko, A. V., Ho, L. C., Barth, A. J., & Leonard, D. C. 2000, AJ, 120, 1499

Murdin, P., & Clark, D. H. 1979, MNRAS, 189, 501

Murray, S. S., Slane, P. O, Seward, F. D., Ransom, S. M., & Gaensler, B. M. 2001, preprint (astro-ph/0108502)

Nomoto, K., et al. 1982, Nature, 299, 803

Pacini F., & Salvati, M. 1973, ApJ, 186, 249

Park, S. et al. 2001a, ApJ, in press (astro-ph/0111116)

Park, S. et al. 2001b, ApJ, 564, L39

Pence, W. D., Snowden, S. L., Mukai, K., & Kuntz, K. D. 2001, ApJ, 561, 189

Petre, R., Becker, C. M., & Winkler, P. F. 1996, ApJ, 465, L43

Pun, C. S. J., et al. 2001, in preparation

Rees, M. J., & Gunn, J. E. 1973, MNRAS, 167, 1

Reynolds, S. P., & Chevalier, R. A. 1984, ApJ, 278, 630

Reynolds, S. P., & Fix, J. D. 1987, ApJ, 322, 673

Sankrit, R., & Hester, J. J. 1997, ApJ, 491, 796

Seward, F. D., Harnden, F. R., Jr., & Helfand, D. J. 1984a, ApJ, 287, L19

Seward, F. D., Harnden, F. R., Jr., Symkowiak, A., & Swank, J. 1984b, ApJ, 281, 650

Sollerman, J., Cumming, R. J., & Lundqvist, P. 1998, ApJ, 493, 933

Sollerman, J., et al. 2000, ApJ, 537, 861

Stollman, G. M. 1987, A&A, 178, 143

Tananbaum, H. 1999, IAU Circ. 7246

Thompson, C., & Duncan, R. C. 1995, MNRAS, 275, 255

Torii, K., Tsunemi, H., Dotani, T., & Mitsuda, K. 1997, ApJ, 489, L145

Tuohy, I., & Garmire, G. 1980, ApJ, 239, L107

Turatto, M., et al. 1998, ApJ, 498, L129

Winkler, P. F., & Kirshner, R. P. 1985, ApJ, 299, 981

van der Swaluw, E., Achterberg, A., Gallant, Y. A., & Tóth, G. 2001, A&A, 380, 309

van der Swaluw, E., & Wu, Y. 2001, ApJ, 555, L49