text:
# Simulating Chandra observations of galaxy clusters

## 1 Introduction

Theoretical studies of the dynamical processes underlying the formation and evolution of galaxies and galaxy clusters are now mainly based on the results of numerical hydro-N-body simulations. Since the pioneering attempts of solving simple N-body problems in the early ’70s, nowadays simulation codes have dramatically improved. Besides the gas hydrodynamics, their most recent versions may account for some complex physical processes, that include, but are not limited to, gas cooling and heating, star formation and feedback, thermal conduction, magnetic fields (see, e.g., Lewis et al. 2000; Yoshida et al. 2002; Muanwong et al. 2002; Dolag, Bartelmann & Lesch 2002; Marri & White 2003; Kay, Thomas & Theuns 2003; Springel & Hernquist 2003; Tornatore et al. 2003; Borgani et al. 2003). Simulation results clearly show that, during their evolution, clusters of galaxies experience violent events that release an enormous amount of energy in the intracluster medium (ICM). These events induce strong, but transient, variation of both ICM density and temperature, so that their distribution is far to be smooth and spherical symmetric.

From the observational side, thanks to the superb angular and spectral resolution of the latest generation X-ray satellites, *Chandra* and XMM-Newton, we now know that many galaxy clusters, including the ones previously identified as relaxed, actually present a great deal of spatial features and have a rather complex thermal structure. Some of these include cold fronts (see, e.g., Abell 2142, Markevitch et al. 2000; Abell 3667, Vikhlinin, Markevitch & Murray 2001; RX J1720+26, Mazzotta et al. 2001; A1795, Markevitch, Vikhlinin & Mazzotta 2001; 2A 0335+096, Mazzotta, Edge & Markevitch 2003), X-ray cavities (i.e. Hydra A, McNamara et al. 2000; Perseus, Fabian et al. 2000; Abell 2052, Blanton et al. 2001; Abell 2597, McNamara et al. 2001; MKW3s, Mazzotta et al. 2002b; RBS797, Schindler et al. 2001; Abell 2199, Johnstone et al. 2002; Abell 4059, Heinz et al. 2002; Virgo, Young, Wilson & Mundell 2002; Centaurus, Sanders & Fabian 2002; Cygnus A, Smith et al. 2002), X-ray blobs and/or filaments (see, e.g., Abell 1795, Fabian et al. 2001; Abell 3667, Mazzotta, Fusco-Femiano & Vikhlinin 2002a; 2A 0335+096, Mazzotta et al. 2003). Furthermore, recent observations indicate that also the gas metal content may have a rather complex distribution that may lead to what we observe as an off-centre peaked metallicity profile (see, e.g., Persueus, Schmidt, Fabian & Sanders 2002, and Churazov et al. 2003; Centaurus, Sanders & Fabian 2002; 2A 0335+096, Mazzotta et al. 2003). Due to their complex nature and the impossibility of using simple deprojection techniques, most of these observed features can be quantitatively studied only through a direct comparison to numerical simulations. Ideally, to make these comparisons straightforward one needs to re-process the simulations themselves through a sort of virtual observatory in such a way that the information provided is as much as possible similar to what an observer can obtain through real X-ray observations of clusters.

In this paper we present X-ray MAp Simulator (X-MAS), a software package we developed to simulate X-ray observations of galaxy clusters obtained from hydro-N-body simulations. The main characteristic of our code is that, giving as input any hydro-N-body simulation, it produces as output an event file which is completely similar to what an X-ray observer would obtain from a real observation. This means that the simulated data can be analysed in the same way and using the same tools of real observations. For the moment our software package simulates ACIS-S3 *Chandra* observations only. In the future we will extend the code to simulate *Chandra* in the ACIS-I mode and XMM-Newton observations with both EPIC and MOS detectors.

The outline of the paper is as follow. In Section 2 we present the general characteristics of X-MAS. In Section 3 we show the results of a simulation of an ACIS-S3 observation of a galaxy cluster. In Section 4 we discuss the possible discrepancy between the projected temperature derived from the spectral analysis of the observation and the emission-weighted temperature directly obtained from the hydro-N-body simulation. In Section 5 we show some applications of X-MAS, namely temperature profiles and maps. Finally in Section 6 we give our conclusions.

## 2 X-ray MAp Simulator: The method

In this section we describe how the X-MAS package works. The package can be divided into two main units. The first unit is quite general and does not depend on the specific characteristics of the X-ray telescope. For each considered energy channel, it generates a corresponding map of the differential flux obtained by projecting the specific emission of each particle along the line of sight. The resulting information for the angular position and energy is stored in a three-dimensional array. An equivalent way to describe the task of this first unit is to say that, for each line of sight within a defined field of view, it calculates and stores the corresponding projected mass-weighted spectrum. The second unit takes each spectrum calculated by the first one and simulates the data relevant to an observation with a specific X-ray telescope and detector for a defined amount of time. Of course, this second unit strongly depends on the characteristics of the X-ray telescope and detector we consider. At present our software package simulates *Chandra* observations in ACIS-S3 configuration only. The application of our simulation package to simulate *Chandra* in the ACIS-I mode and XMM-Newton observations with both EPIC and MOS detectors requires an adaptation of this second unit only. In the following we describe in details how the two package units work.

### 2.1 First Unit: generating differential flux maps

The first unit of our package X-MAS requires as input the output of an hydro-N-body simulation. After selecting the direction and the depth of the galaxy cluster for which we want to simulate the observation, the program generates the cluster projected spectra corresponding to all the lines of sight in a defined field of view. This is simply done by considering an energy interval $`[E_{\mathrm{min}},E_{\mathrm{max}}]`$, which we divide in $`n__E`$ regular energy channels with energy width $`\mathrm{\Delta }E=(E_{\mathrm{max}}E_{\mathrm{min}})/n__E`$. The energy interval and the energy resolution selected for unit one of the program need to be higher than or equal to the energy response and the energy resolution of the instrument we intend to simulate later with unit two, respectively. We use $`E_{\mathrm{min}}=0.1`$ keV, $`E_{\mathrm{max}}=10`$ keV and $`n__E=495`$ (or equivalently $`\mathrm{\Delta }E=20`$ eV111 The energy response and resolution of *Chandra* are \[0.1,10.0\] keV and $`\mathrm{\Delta }E100`$ eV, respectively (see the *Chandra* Proposers’ Observatory Guide; http://asc.harvard.edu/udocs/docs/docs.html)).

For each of these channels we calculate the two-dimensional map of the corresponding differential flux produced by the simulated galaxy cluster simply by projecting on a regular grid (1024 pixels $`\times `$ 1024 pixels) the flux $`F_i^\nu `$ of each particle and summing over all particles. The differential flux for the cluster is finally stored in a three-dimensional array in which two dimensions represent the angular coordinates and the third one is for the energy.

In order to calculate the flux associated with each particle we follow a procedure very similar to the one described in Mathiesen & Evrard (2001). Starting from its three-dimensional position $`𝐱__i`$, mass $`m_i`$ and density $`\rho _i`$, we assign to the $`i`$-th gas particle in the simulation an effective volume $`V_i=m_i/\rho _i`$, which is assumed for simplicity to be cubic and centred on $`𝐱__i`$.

For cosmological sources, the standard relation between flux $`F`$ and luminosity $`L`$ holds: $`F=L/4\pi d_L^2`$, where $`d_L`$ is the luminosity distance (depending on cosmology and on the source redshift $`z`$). Considering differential quantities, the previous relation becomes $`F_\nu =[(1+z)L_{\nu (1+z)}]/4\pi d_L^2`$, where the $`\nu `$ pedex accounts for the dependence on the energy band. If the flux is expressed in terms of incoming photons instead of energy, we can introduce the quantity $`F_\nu ^\gamma F_\nu /h\nu `$ (in photon/s/cm2/keV), so that $`F_\nu ^\gamma d(h\nu )`$ represents the flux of the incoming photons with energy in the range $`[h\nu ,h(\nu +d\nu )]`$ keV. In general, we prefer to give quantities like luminosity or flux in terms of incoming photons or counts (instead of energy) because they can be more easily related to real observational data. Finally the relation between flux and luminosity simply becomes

$$F_\nu ^\gamma =\frac{(1+z)^2L_{\nu (1+z)}^\gamma }{4\pi d_L^2}.$$
(1)
The emissivity per unit of frequency $`ϵ_\nu `$ in a region of volume $`V`$ is related to the luminosity as $`L_\nu =_Vϵ_\nu 𝑑V^{}`$; a similar relationship holds for differential and photon quantities.

For a galaxy cluster the emissivity from a small enough region of plasma is given by a single temperature thermal model. Given the electron and hydrogen densities ($`n_e`$ and $`n__H`$, respectively), the emissivity can be written as $`ϵ=n_en__HP(T,Z)`$. The quantity $`P(T,Z)`$ is usually called power coefficient and depends only on the temperature $`T`$ and metallicity $`Z`$ of the gas. Using the relation above, the photon luminosity can be written as $`L_\nu ^\gamma =_Vϵ_\nu ^\gamma 𝑑V^{}=P_\nu ^\gamma EM`$; the quantity $`EM_Vn_en__H𝑑V^{}`$ is often referred to as emission measure. The final relation between flux and power coefficients is then

$$F_\nu ^\gamma =\frac{(1+z)^2}{4\pi d_L^2}EMP_{\nu (1+z)}^\gamma (T,Z).$$
(2)
Using the temperature and the metallicity of each particle in the simulation, we calculate $`P(T_i,Z_i)`$ using the single temperature thermal model mekal (see, e.g., Kaastra & Mewe 1993; Liedahl, Osterheld & Goldstein 1995, and references therein) implemented in the utility XSPEC (Arnaud 1996).

To make our simulator of X-ray observation more complete, we allow to include the effects on the spectra induced by the Galactic HI absorption. This is done, at the end of the procedure, by multiplying the flux in each energy channel by an absorption coefficient given by the wabs model (Morrison & McCammon 1983), once a value for the column density $`N__H`$ is assumed.

### 2.2 Second Unit: Simulating Chandra ACIS-S3 observations

The second unit of our simulation package X-MAS takes as input the projected spectra produced by the first unit and, after convolving them with the technical characteristics of a specific instrument, generates an event file similar to the one obtained from a real observation. To do that we use the data simulation command fakeit in the utility XSPEC (see, e.g., Xspec User’s Guide version 11.2.x; Dorman & Arnaud 2001222http://legacy.gsfc.nasa.gov/docs/xanadu/xspec/manual). The above command creates simulated data from the input spectral model by convolving it with the ancillary response files (ARF) and the redistribution matrix files (RMF), which fully define the response of the considered instrument, and by adding noise appropriate to the specified integration time. Once the data of all spectra have been simulated, we generate a photon event file satisfying the standards defined for real observations. This is quite important because it allows our mock observations to be analysed by using the same tools and procedures used for the real ones.

In the following we describe how we use the second unit of our software package to simulate X-ray observations performed using *Chandra* with the back illuminated CCD ACIS-S3.

The ACIS-S3 detector is a square with a grid of 1024 pixels $`\times `$ 1024 pixels, and a field of view of about 8.3 arcmin by side. The nominal angular resolution is about 0.5 arcsec/pixel. Although the instrumental response is position-dependent, it is quite constant over CCD subregions of $`32\times 32`$ pixels, which for convenience we call detector tiles (or tiles for short). To fully account for the instrument response the *Chandra* calibration team produced 1024 ARF and 1024 RMF, one for each of the detector tile above.

Detection events have to preserve the spatial and spectral information. As explained above, their number and energy are obtained by executing the command fakeit of the utility XSPEC using the appropriate ARF and RMF. To account for the instrumental background we provide to the fakeit command an appropriate background file extracted from the blank-sky background dataset (Markevitch 2001333http://asc.harvard.edu/ “Instruments and Calibration”, “ACIS”, “ACIS Background”) in the same tile subregion corresponding to the spectrum to be simulated.

To speed up the simulation process we use an adaptive algorithm: before simulating the data of each spectrum, we estimate the expected number counts associated with each detector tile. If this number is lower than a given threshold we generate the events, otherwise we iteratively subdivide the region in four squares, until the threshold is reached. The spatial position of the events obtained in this way is then reconstructed by randomly distributing the simulated photons inside the region by using a weight proportional to the original flux.

## 3 ACIS-S3 observation of a simulated cluster

As a first example of possible applications of our software package, we generate a 300 ks *Chandra* ACIS-S3 observation of a high-resolution hydro-N-body simulation of a galaxy cluster. The object was selected from a sample of 17 objects obtained using the technique of re-simulating at higher resolution a patch of a pre-existing cosmological simulation. The assumed cosmological framework is a cold dark matter model in a flat universe with a present matter density parameter $`\mathrm{\Omega }_m=0.3`$ and a contribution to the density due to the cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$; the baryon content corresponds to $`\mathrm{\Omega }_B=0.03`$; the value of the Hubble constant (in units of 100 km/s/Mpc) is $`h=0.7`$, and the power spectrum normalization is given by $`\sigma _8=0.9`$. The re-simulation method, called ZIC (for Zoomed Initial Conditions), is described in detail in Tormen, Bouchet & White (1997), while an extended discussion of the properties of the whole sample of these simulated clusters is presented elsewhere (Tormen, Moscardini & Yoshida 2003; Rasia, Tormen & Moscardini 2003). Here we remind only some of the characteristics of the cluster used in this paper. It has been obtained by using the publicly available code GADGET (Springel, Yoshida & White 2001); during the run, starting at redshift $`z_{\mathrm{in}}=35`$, we took 51 snapshots equally spaced in $`\mathrm{log}(1+z)`$, from $`z=10`$ to $`z=0`$. Its virial mass at $`z=0`$ is $`1.46\times 10^{15}h^1M_{}`$, corresponding to a virial radius of $`2.3h^1`$ Mpc; the mass resolution is $`4.5\times 10^9h^1M_{}`$ per dark particles and $`5\times 10^8h^1M_{}`$ per gas particles; the total number of particles found inside the virial radius is 566,374, 48 per cent of which are gas particles. The gravitational softening is given by a $`5h^1`$ kpc cubic spline smoothing. Since this particular simulation does not provide information on the cluster metallicity, we fixed its value to $`Z=0.3Z_{}`$. Furthermore, we assumed a galactic equivalent column density of $`N_H=5\times 10^{20}`$ cm-2. Among the available snapshots at different redshifts we chose to observe the one at $`z=0.21`$. At this redshift the cluster, having a virial mass of $`1.05\times 10^{15}h^1M_{}`$ and a virial radius of $`2.2h^1`$ Mpc, is undergoing several merger events, so its structure is quite complex.

The cluster flux map, obtained with the first unit of our X-MAS, is presented in Fig. 1. The figure shows the flux in the \[0.1,10.0\] keV energy interval. The displayed region is 8.3 arcmin, which corresponds to approximately 1.7 (proper) Mpc at the cluster redshift. The superimposed isocontours are obtained from the image after a Gaussian smoothing with $`\sigma =8^{\prime \prime }`$. Levels are spaced by a factor of $`2`$ with the highest level corresponding to $`1.8\times 10^{15}`$ erg cm-2 s-1. In the external part of the flux map we find the presence of a number of merging subclumps, confirming the highly perturbed dynamical phase of the cluster. Notice that the flux map also shows an orange-skin-like texture induced by angular structures on scales of the order of few arcsec. This small-scale structure is actually an artifact that we deliberately introduced by considering the particle nature of the Smoothed Particle Hydrodynamics simulation used. For the purpose of this paper this artifact turns useful to test the imaging performance of our simulator on very small scales.

In Fig. 2 we show the photon image of the 300 ks *Chandra* ACIS-S3 observation of the simulated cluster shown in Fig. 1, as obtained by applying the second unit of our software X-MAS. The image, extracted from the event file in the \[0.3,9.0\] keV energy band, is background-subtracted, vignetting-corrected, and binned to $`1\mathrm{}`$ pixels. We notice that, as expected, after being observed with *Chandra*, the spatial features present in the simulation, but fainter than the instrument background, are no longer detected. For example, the three faint subclumps on the North, North-West and South-West, clearly visible in Fig. 1, have been washed out in Fig. 2. Conversely all the brighter features appear to be well reproduced by our simulator. We find that this is true even on scales as small as few arcsecs: the previously mentioned angular artifacts of the simulation on these scales are well visible in the photon image, indeed.

## 4 Spectral analysis: preliminary tests

The spectral analysis of *Chandra* observations of the simulated clusters is done by applying standard procedures and tools used for real observations. In particular, in order to extract spectra, we use the dmextract tool of the CIAO software (see http://cxc.harvard.edu/ciao/). Spectra are extracted in the \[0.6,9.0\] keV band in PI channels, re-binned to have a minimum of 10 counts per bin and fitted using the XSPEC package (Arnaud 1996). As background we use the spectrum extracted from the publicly available background dataset described in Section 2.2. The position-dependent RMFs and ARFs are computed and weighted by the X-ray brightness over the corresponding image region using the “calcrmf” and “calcarf” tools444A. Vikhlinin 2000 (http://asc.harvard.edu/ “Software Exchange”, “Contributed Software”).. Spectra are fitted with a single temperature absorbed mekal model using the C-statistics (Cash 1979; Arnaud 1996) and the Levenberg-Marquardt minimization method (i.e. the XSPEC default minimization method). In the fit procedure we fix the cluster redshift, metallicity, and hydrogen column density to the values used as inputs to compute the *Chandra* observation. As result of the fit we obtain what, from now on, we call projected spectroscopic temperature and its 68 per cent confidence level error for one interesting parameter, $`T_{\mathrm{spec}}`$ and $`\sigma _{\mathrm{spec}}`$, respectively.

Our goal is to use our *Chandra* simulator to compare the spectroscopic temperature $`T_{\mathrm{spec}}`$ to one of the possible temperature estimator adopted in the analysis of the results of hydro-N-body simulations. In particular we use the emission-weighted temperature $`T_{\mathrm{sim}}`$, which is the one most commonly adopted. This is defined as:

$$T_{\mathrm{sim}}\frac{WT𝑑V}{W𝑑V}.$$
(3)
In the previous equation, $`T`$ is the cluster gas temperature, $`dV`$ is the volume along the line of sight, while the characteristics of the emissivity are taken into account by the weight $`W`$, which is usually defined as $`W=\mathrm{\Lambda }(T)\rho ^2`$, where $`\mathrm{\Lambda }(T)`$ is the cooling function and $`\rho `$ the gas density (see, e.g., Navarro, Frenk & White 1995). The emission-weighted temperature for a specific cluster region is derived directly from the original hydro-N-body simulation by:

$$T_{\mathrm{sim}}=\frac{_im_i\rho _i\mathrm{\Lambda }(T_i)T_i}{_im_i\rho _i\mathrm{\Lambda }(T_i)},$$
(4)
where the sums are extended to all the particles which are projected inside the considered region. To measure the deviation from a completely isothermal distribution for the gas within the projected region, we use the relative emission-weighted temperature dispersion, defined as

$$\frac{\sigma __T}{T_{\mathrm{sim}}}=\frac{1}{T_{\mathrm{sim}}}\left[\frac{WT^2𝑑V}{W𝑑V}T_{\mathrm{sim}}^2\right]^{1/2}:$$
(5)
quasi-isothermal and highly perturbed regions will have low and high values of $`\sigma __T/T_{\mathrm{sim}}`$, respectively.

In the following subsections, we will compare $`T_{\mathrm{spec}}`$ to $`T_{\mathrm{sim}}`$ in two specific cases: i) an ideal toy-model corresponding to a perfect isothermal cluster, and ii) a more realistic cluster with a complex thermal structure.

### 4.1 Toy isothermal cluster

Here we discuss how the spectroscopically inferred temperature $`T_{\mathrm{spec}}`$ is affected by the total number of detected photons. At this goal we take our simulated cluster and set the temperature of all gas particles to an arbitrarily chosen constant value, $`T_{\mathrm{sim}}`$, ranging from 3 to 12 keV. After producing a 300 ks *Chandra* observation of this isothermal cluster, we extract all the spectra from each of the $`32^2`$ pixels tile regions, defined in Section 2.2. As the total number of photons in each spectrum is proportional to the total cluster emissivity in the tile region where it was extracted from, this simple procedure allows us to produce a distribution of spectra with the same temperature, but different total photon counts. Among the extracted spectra we select only the ones for which the source flux is higher than the background. For the exposure time and the fixed region size here considered, this is equivalent to selecting spectra with net total counts $`N_\gamma >250`$.

Consistently with what already discussed in literature (see, e.g., Nousek & Shue 1989 and references therein), we find a small bias between the input temperature and the best fit value of the temperature obtained from spectra with small total number counts ($`N_\gamma <10^3`$). This bias increases if the total number counts gets smaller and/or the temperature is higher. Nevertheless, we find that the input temperature value is always consistent within the 1$`\sigma `$ errors associated with the spectroscopic temperature.

### 4.2 Realistic cluster with a complex thermal structure

To study the effect of temperature inhomogeneities on the final projected temperature estimates we repeat the analysis of the previous subsection using the actual photon temperatures obtained from the simulation, instead of forcing them to be isothermal. As in § 4.1 we selected only spectra for which the source flux is higher than the background (net total counts $`N_\gamma >250`$). We estimate $`T_{\mathrm{spec}}`$ and its error $`\sigma _{\mathrm{spec}}`$ by fitting the data with an absorbed single temperature thermal model.

The result of our analysis is shown in Fig. 3. All the considered spectra are divided in three bins with increasing relative emission-weighted temperature dispersion $`\sigma __T/T_{\mathrm{sim}}`$ (see equation 5). The three bins, containing the same number of spectra, correspond to regions that, from the point of view of the temperature distribution, are highly homogeneous (i.e. almost isothermal), mildly inhomogeneous and highly inhomogeneous, respectively. For each spectrum we calculate the temperature discrepancy between the spectroscopic and the emission-weighted estimates, $`(T_{\mathrm{spec}}T_{\mathrm{sim}})/T_{\mathrm{sim}}`$. The points in the figure represent the median discrepancy in each bin and are located at the median value of $`\sigma __T/T_{\mathrm{sim}}`$. The vertical error bars represent the 25 and 75 percentiles of the distribution of the discrepancy $`(T_{\mathrm{spec}}T_{\mathrm{sim}})/T_{\mathrm{sim}}`$ inside each bin. It is interesting to notice that for regions where the temperature structure is not highly inhomogeneous (i.e. $`\sigma __T/T_{\mathrm{sim}}<0.5`$), on average we find that $`T_{\mathrm{spec}}`$ and $`T_{\mathrm{sim}}`$ are consistent, although it is still possible to obtain discrepancies of the order of 10-15 per cent. Conversely, for highly inhomogeneous regions (i.e. $`\sigma __T/T_{\mathrm{sim}}>0.5`$), we find a larger discrepancy between the spectroscopic and the emission-weighted temperatures, with on average $`T_{\mathrm{spec}}`$ systematically lower than $`T_{\mathrm{sim}}`$.

Concluding this section, we stress that our results clearly show that, when strong temperature inhomogeneities are present, the spectroscopic and the emission-weighted temperature measurements are likely to be inconsistent with each other. Consequently, much attention must be paid in the theoretical interpretation of observational temperatures. Similar conclusions have been reached by an equivalent analysis of Mathiesen & Evrard (2001) in which they compare the overall spectroscopic and emission-weighted temperatures of a ensemble of 24 simulated clusters of galaxies. Although with a large scatter, they claim that the spectroscopic temperature is 20 per cent lower than the emission-weighted one.

## 5 Examples of applications

In this section we discuss, as examples, two different applications of our simulation method: the computation of temperature profiles and the production of projected temperature maps of galaxy clusters.

### 5.1 Temperature profiles

One of the standard ingredients of the theoretical modeling of galaxy clusters is the assumption of an isothermal distribution of gas. For example, this is often used to obtain estimates of their mass by using the equation of hydrostatic equilibrium. Another standard application is to directly relate the observed cluster temperature function to the theoretically estimated mass function to obtain constraints on the main cosmological parameters, as the matter density parameter and the normalization and shape of the power spectrum of primordial fluctuations.

However, the results of high-resolution hydro-N-body simulations indicate that the temperature shows radial gradients which cannot be neglected in a dynamical analysis (see, e.g., Evrard, Metzler & Navarro 1996; Eke, Navarro & Frenk 1998; Rasia, Tormen & Moscardini 2003; Borgani et al. 2003). This is confirmed by recent observational data, mainly obtained using the X-ray observatories ASCA, BeppoSAX, and, more recently, *Chandra* and XMM-Newton (see, e.g., Markevitch et al. 1998; De Grandi & Molendi 2002; Pratt & Arnaud 2002).

In this subsection we use our simulator X-MAS to derive the projected temperature profiles of both a major merging system and a quasi relaxed one. These profiles are compared to the emission-weighted temperature profiles obtained directly from the hydro-N-body simulations. For the merger system we consider the output of the hydro-N-body simulation at redshift $`z=0.21`$, already described in Section 3 (see Fig. 2). For the relaxed system we use the output of the same simulation at $`z0.33`$, when the galaxy cluster is not undergoing any strong merging events. At this redshift the cluster has a virial mass of $`3.87\times 10^{14}h^1M_{}`$ and a virial radius of $`1.7h^1`$ Mpc.

To measure the projected temperature profile we extract the spectra from circular annuli centred on the cluster X-ray peak. The size of the bin has been chosen in order to have approximately the same number of photons inside each annulus. Again, spectra are extracted in the \[0.3,9.0\] keV energy band and fitted with a single temperature absorbed mekal model with the values for $`N__H`$, metallicity, and redshift fixed at the simulated values. The spectroscopic temperature profiles $`T_{\mathrm{spec}}`$, together with their relative 68 per cent confidence level errors $`\sigma _{\mathrm{spec}}`$, are shown as filled circles in Fig. 4. Left and right panels refer to the relaxed and merging systems, respectively. In the same figure we show the emission-weighted temperature $`T_{\mathrm{sim}}`$ and the corresponding value of $`\sigma __T`$, computed as in equation 5. As discussed in the previous section, the value of $`\sigma __T`$ measures the degree of thermal inhomogeneity of the considered cluster region. From Fig. 4 we see that, as expected, the projected radial thermal structure of the cluster in its relaxed phase is far more homogeneous than the structure in the perturbed one. This different degree of thermal homogeneity has strong implications on the temperature profiles. In fact, we notice that while for the relaxed phase the spectral and the emission weighted temperature profiles are in good agreement, this is not longer true for the perturbed phase. Furthermore, we confirm the systematic trend previously discussed: the spectral temperatures are lower than the emission-weighted temperatures.

### 5.2 Temperature maps

As a further example of possible applications of our X-ray observatory simulation package, we now compare the emission-weighted temperature map determined from the cluster simulation to the projected temperature map obtained from the spectral analysis of its corresponding *Chandra* observation. In this section we focus only on the output of the hydro-N-body simulation at $`z0.21`$, i.e. the perturbed phase considered in Section 3.

First, we calculate the emission-weighted temperature map of the simulated clusters. As the simulation provides us with the density and the temperature of each particle, the emission-weighted temperature map for the simulated cluster can be obtained with the same spatial resolution of its X-ray image. The result is shown in Fig. 5. For reference, in the same figure we superimpose the isocontours corresponding to the cluster flux distribution shown in Fig. 1. As already evident from the temperature profile (see Fig. 4), the internal region is far from being isothermal. Between the two colder central subclumps with $`T56`$ keV, we notice a region with higher gas temperature, $`T10`$ keV, which corresponds to the gas compressed by the merging blobs. Moreover we notice that the more external merging subclumps are also significantly colder than the cluster ambient gas. Of particular interest are the two subclumps on the lower-left corner of the image: they seem to form a single structure of cold gas. In addition, we notice the presence of a shock front with a post-shock gas temperature of approximately 20 keV which is produced by the motion toward the cluster centre of the most internal of these two clumps.

To calculate the spectroscopic projected temperature map we use the 300 ks *Chandra* observation of the simulated cluster described in Section 3. We extract spectra using adjacent square regions. Over most of the map the extraction regions coincide with the tile regions defined in Section 2.2. However, in the outskirt of the cluster, where the surface brightness is lower, we use larger extraction regions which are obtained by combining two or more tile regions in such a way that the total net number of photons per spectrum is $`N_\gamma >250`$. Each spectrum is then fitted with an absorbed single temperature mekal model with $`N__H`$, redshift, metallicity fixed to the input values.

The projected spectroscopic temperature map $`T_{\mathrm{spec}}`$ is shown in the left panel of Fig. 6. In order to make a direct comparison of the spectroscopic temperature map to the emission-weighted one, we degenerate the resolution of the latter to match the former. Thus, in the right panel of Fig. 6 we report the same emission-weighted temperature map shown in Fig. 5, but re-binned as the spectroscopic temperature map of the left panel. It is worth to notice that, although with a lower spatial resolution, the re-binned emission-weighted temperature map presents all the main temperature structures described before. In particular we see that the two central blobs are cold and we can easily identify the compression-heated gas between them. At the same way we can recognize the two cold subclumps on the left-bottom corner, as well as the presence of the shock-heated gas in front of the innermost of the two which is moving toward the cluster centre.

If we now compare the emission-weighted temperature map on the right to the spectral temperature map on the left, we notice that, although qualitatively similar, they show a number of important differences. Among others, we point our attention on the fact that: i) the central merging blobs appear to be colder in the observed spectroscopic temperature map than in the emission-weighted one; ii) the shock front produced by the motion of the innermost subclump in the lower-left corner, clearly visible in the emission-weighted map, is no longer detected in the spectroscopic temperature map.

To better visualize the temperature differences among the two previous maps, in the left panel of Fig. 7 we show the spatial distribution of the difference $`T_{\mathrm{sim}}T_{\mathrm{spec}}`$. We find temperature differences $`|T_{\mathrm{sim}}T_{\mathrm{spec}}|>1`$ keV for $`50`$ per cent of the pixels. In particular we notice that most of these differences are such that $`T_{\mathrm{sim}}>T_{\mathrm{spec}}`$. To quantify the significance of this discrepancy, in the right panel of Fig. 7 we present the map of $`(T_{\mathrm{sim}}T_{\mathrm{spec}})/\sigma _{\mathrm{spec}}`$ (being $`\sigma _{\mathrm{spec}}`$ the 68 per cent confidence level error associated with the spectroscopic temperature measurement). From this plot we see that most of the temperature differences are significant at $`3\sigma `$ confidence level. In particular the temperature discrepancies previously noted for the central subclumps and for the shock regions are significant at $`9\sigma `$.

## 6 Discussion and Conclusions

In this paper we presented a numerical technique and the relative software package X-MAS devoted to simulate X-ray observations of galaxy clusters obtained from hydro-N-body simulations. As specific application, we used this technique to simulate *Chandra* ACIS-S3 observations.

We stress that one of the main features of our code is that it generates event files following the same standards used for real observations. This is extremely important as it implies that our simulated observations can be analysed in the same way and with the same tools of real observations.

Using the results of high-resolution hydro-N-body simulations, we generated the *Chandra* observations of a number of simulated clusters. By performing a standard spectral analysis, we derived the projected spectral temperature in specific cluster regions. The spectral temperature has been finally compared to the emission-weighted temperature commonly used to describe to cluster gas thermal properties in numerical works.

Our main finding is that the two temperature estimates are likely to show a significantly large discrepancy, the spectroscopic temperature being lower than the emission-weighted temperature. This effect is more evident if the thermal structure of the cluster within a particular projected region is relatively complex (i.e. the region is thermally highly inhomogeneous). We point out that the data analysis procedure used in this paper assumes that both the background level and spectral shape are well known. In fact the background used to produce the *Chandra* observations is the same used in the spectral analysis. Uncertainties in the background level and spectral shape, which are quite likely in real observations, would inevitably result in an increase of the estimated errors and in discrepancies between the spectroscopic and emission-weighted temperatures even larger than what has been shown in this paper. Regardless of the background, the main reason behind the observed temperature discrepancy can be easily explained if one considers that spectroscopically the temperature is determined by fitting a thermal model (in particular its bremsstrahlung component) to the observed spectrum. The point is that the sum of two bremsstrahlung spectra with similar emission but different temperatures $`T_1`$ and $`T_2`$ is no longer a bremsstrahlung with a given temperature $`T_3`$. In fact,

$$\mathrm{exp}(E/T_1)+\mathrm{exp}(E/T_2)A\mathrm{exp}(E/T_3),$$
(6)
unless $`T_1=T_2`$. When such a combined spectrum is fitted with a single bremsstrahlung spectral model, we obtain a temperature which is not exactly the mean value of $`T_1`$ and $`T_2`$, though it will be an intermediate value between the two. The larger is the difference $`|T_2T_1|`$, the bigger will be the discrepancy between the spectroscopic and the mean temperature. This result is important because it implies that much attention must be paid in the theoretical interpretation of observational temperatures. Just to give a more detailed idea of the problem, in Section 5 we discussed the implications for two specific cases in which such comparisons have been done in the past: temperature profiles and temperature maps.

Hydro-N-body simulations have been (and still are) largely used to derive universal temperature profiles for galaxy clusters. These profiles are directly compared to the ones obtained from real observations. In Section 5.1 we used the *Chandra* simulator to compare spectroscopic and emission-weighted temperature profiles of two extreme phases of the cluster evolution: an almost relaxed system and a highly perturbed object. In agreement with what said before, we find that if the cluster is relaxed, the emission-weighted temperature profile agrees with the spectroscopic one. However, if the system is highly perturbed this is not longer true. This result indicate that direct comparisons of observed and “simulated” temperature profiles are not fully justified. In Section 5.2 we show that a similar problem applies also to the temperature maps.

In a forthcoming paper we intend to use X-MAS to investigate in detail a number of problems related to the actual observations of X-ray galaxy clusters. In particular we intend to study the complex problem of spectral deprojection and relative cluster mass determination. Moreover we are planning to use X-MAS to perform observations of simulated clusters that account for the processes of the metal enrichment of the ICM. This will allow us to address the even more complex problem of determining both the mean and the single-element metallicity structure of real galaxy clusters.

To conclude we stress that the possible applications of our simulations of observations are quite vast. Besides the here discussed problem of the comparisons of outputs of numerical simulations to real observed galaxy clusters, it can be very useful to better plan observations with existing X-ray telescopes and even more to verify what will be the real capabilities of futures ones. At this goal, results obtained by using X-MAS will be soon publicly available on the web.

Concerning our software package, we remind that at present we only completed the software module that simulates ACIS-S3 *Chandra* observations. Work is in progress to include further modules to simulate *Chandra* observations in the ACIS-I mode, and XMM-Newton observations with both EPIC and MOS detectors.

## Acknowledgments

This work has been partially supported by Italian MIUR (Grant 2001, prot. 2001028932, “Clusters and groups of galaxies: the interplay of dark and baryonic matter”) and ASI. We thank also K.A. Arnaud, S. Borgani, J. McDowell, M. Meneghetti and M. Wise, for clarifying discussions. PM, GT and LM are grateful to the Aspen Center for Physics, where the paper has been discussed and partially written up.


format: text

text:
# An L′-band survey for circumstellar disks around low-mass stars in the young 𝜎 Orionis cluster

## 1 Introduction

Disk-like structures are believed to be ubiquitous around young protostars. These disks are dissipated very early in pre-main-sequence (PMS) evolution, perhaps by powerful stellar jets/outflows or photodissociation by the far-ultraviolet flux from nearby massive OB stars. Despite their short lives, the timescales and mass dependence of disk dissipation have far reaching consequences in astrophysics: the efficiency of disk depletion could be the strongest factor in determining the timescales on which planets form in a particular stellar system (Haisch, Lada & Lada, 2001a), or whether they form at all (Brandner et al., 2000). Disks probably play a significant role in early angular momentum regulation and the dissipation timescale is thought to control the spread in rotation rates of young stars (Sills, Pinsonneault & Terndrup, 2000). Stars may accrete a significant fraction of their final mass from a circumstellar disk, so the timescale and mass dependence of that accretion influences PMS evolution and thus attempts to estimate ages and masses from evolutionary PMS models (Comerón et al., 2003). The mass dependence of disk frequencies can provide a stern test for low-mass stellar and brown dwarf formation theories. For instance, models involving competitive accretion and subsequent ejection of brown dwarfs from protostellar aggregates (Reipurth & Clarke, 2000; Bate et al., 2003) may imply shorter disk dissipation times for the lower mass fragments.

Observed disk frequencies in samples of young stars with different ages, masses and environments provide an empirical determination of disk lifetimes. Judging from L-band excesses, young clusters exhibit high disk frequencies ($``$ 80%, e.g. the Trapezium cluster: Lada et al. 2000) up to ages of $``$ 1.5 Myr, which then decrease rapidly with age: at $``$ 3 Myr, 50% of disks have been dissipated, and the timescale for all cluster members to lose their disks may be as short as $``$ 6 Myr (Haisch et al., 2001a). Such timescales have been questioned by a high disk frequency in the 9 Myr $`\eta `$ Chamaeleontis cluster, a sparsely populated cluster with no massive stars (Lyo, Lawson & Mamajek, 2003).

$`\sigma `$ Orionis is a Trapezium-like system with an O9.5 V primary. The population of low-mass stars spatially clustered around this system was discovered as bright X-ray sources in ROSAT images, and follow-up optical spectroscopy confirmed most sources as PMS stars (Wolk, 1996; Walter et al., 1997). This association is young, nearby and affected by low reddening, making it an ideal target to analyse the PMS population even down to brown dwarfs (e.g. Béjar et al., 2001; Barrado y Navascués et al., 2003; Kenyon et al., 2003) and isolated planetary mass objects (Zapatero Osorio et al., 2000). Furthermore, at an age of 3$``$5 Myr (e.g. Oliveira et al., 2002; Zapatero Osorio et al., 2002; Jayawardhana et al., 2003), the $`\sigma `$ Orionis cluster is at a crucial stage in terms of disk evolution and it is therefore a key case to better constrain disk dissipation timescales. Recently, a possible proto-planetary disk, apparently on the process of being dissipated, has been discovered very close to $`\sigma `$ Ori (van Loon & Oliveira, 2003).

The $`K_\mathrm{s}`$-excess disk frequency is 5$`12`$ % for the low-mass and brown dwarf members of the $`\sigma `$ Orionis cluster (Oliveira et al., 2002; Barrado y Navascués et al., 2003). On the other hand, the presence of strong H$`\alpha `$ emission suggests accretion disk frequencies as high as 30% (Zapatero Osorio et al., 2002). However, the most reliable method to determine the disk frequency in a low-mass population is by measuring the $`(KL)`$ colours and deriving colour excesses (e.g. Wood et al., 2002). Jayawardhana et al. (2003) have obtained L-band observations of 6 $`\sigma `$ Ori cluster members, finding two with a $`(KL)`$ excess. The significance of this result is obviously limited by the size of the sample. We have performed L-band (3.8 $`\mu `$m) observations of a representative sample of 28 cluster members, using the newly installed imager UIST at the United Kingdom Infrared Telescope (UKIRT). Young stars are well known for their variability across the spectrum including infrared (IR) wavelengths (Carpenter et al., 2001, 2002), therefore we have obtained nearly simultaneous K-band observations for all our targets. In this paper, we describe the results of this survey, and discuss our derived disk frequency within the framework of disk destruction timescales by comparing with similar surveys in other young clusters.

## 2 Cluster members and properties

### 2.1 Sample of cluster members

We have an on-going program to observe in the L-band $`\sigma `$ Ori cluster members identified at optical wavelengths. We describe here the observations of 28 of the brightest cluster members: their positions, $`I_\mathrm{c}`$ magnitudes, 2MASS (Two Micron All Sky Survey) $`J,H`$ and $`K_\mathrm{s}`$ magnitudes, the new K- and L-band magnitudes and identifications are listed in Table 1. Some sources were first identified as ROSAT X-ray sources and photometric cluster candidates by Wolk (1996, W96) while other objects are photometric candidates identified by Béjar et al. (2001, B01). Zapatero Osorio et al. (2002, ZO02) have spectroscopically confirmed cluster membership for both these sets of objects. The remaining objects are spectroscopic cluster members found by Kenyon et al. (2003, K03). The I-band magnitudes are from either Béjar et al. (2001) or Kenyon et al. (2003). The brighter objects in the sample are mostly from the X-ray selected sample (Wolk, 1996) while the fainter objects were photometrically and spectroscopically selected. In Sect. 5.1 we discuss the effects of possible selection biases on our results.

Searching for circumstellar disks around these 24 objects is the main goal of these observations. To this sample we have added 4 objects that are known IRAS sources in the region (no reference entry in Table 1). They have been confirmed by Oliveira & van Loon (2003) (see also Oliveira et al. 2003) as mid-infrared sources with spectral energy distributions (SEDs) consistent with them being young stars with dusty circumstellar disks; thus, based on their location, youth and infrared excesses they are also likely to be members of the $`\sigma `$ Ori cluster.

### 2.2 Cluster age

The age determination for clusters at very young ages is always an uncertain affair. In OB associations, several methods can be used: the properties of the massive O-stars, low-mass stellar isochrone fitting, and lithium abundance evolution in the association.

The age of the multiple system $`\sigma `$ Ori is estimated to be 1.7$``$7 Myr based on stellar properties and its membership of the Orion OB1b group (Brown, de Geus & de Zeeuw, 1994, and references therein). Using isochrones from several authors, Béjar, Zapatero Osorio & Rebolo (1999, and references therein) determined a low-mass star cluster age between 1$``$5 Myr; using a different sample of low-mass cluster members, Oliveira et al. (2002) determined an isochronal median age of 4.2$`{}_{1.5}{}^{}{}_{}{}^{+2.7}`$ Myr, where the quoted errors are due to the uncertainties of the Hipparcos distance to $`\sigma `$ Ori. Zapatero Osorio et al. (2002) compared the lithium abundance in low-mass cluster members with theoretical predictions from several authors. They found no evidence of appreciable lithium destruction and they inferred an upper limit to the cluster age of 8 Myr.

In Fig. 1 we plot the $`I/(IJ)`$ colour-magnitude diagram for the present sample. We use the Baraffe et al. (1998) evolutionary models with the mixing length parameter set to 1.0 pressure scale height for lower mass objects and to 1.9 for $`M>0.62`$ M. We have opted to use these models because they incorporate model atmospheres and predict PMS magnitudes without the use of empirical effective temperature-colour relations. We have computed the age and mass of each object; two objects (from the 24 with I-band magnitudes) fall outside the Baraffe grid of models — they appear to be very young (age $`<`$ 1 Myr). The objects have masses in the range 1.0$``$0.13 M. The sample has a median age of 3.6 Myr (at the Hipparcos distance of 352 pc), consistent with previous age determinations (see above). The conservatively large uncertainty in the Hipparcos distance (352$`{}_{85}{}^{}{}_{}{}^{+166}`$ pc) produces median age uncertainties as computed by Oliveira et al. (2002).

The most obvious thing from this colour-magnitude diagram is the large apparent age spread ($`<`$ 1 to 20 Myr). Hartmann (2001) has analysed many possible sources of uncertainty in the determined stellar ages in star forming regions. He concluded that observational errors probably account for a large fraction of the observed age spread in such regions. Therefore, it is still not firmly established whether observed age spreads are real effects.

The main sources of uncertainty in the case of the present data set are photometric variability and unresolved binaries. Unresolved binaries would tend to make objects appear brighter and redder, making them look younger than isolated cluster siblings. This is one reason why we consider the true age of the cluster likely to be older than the median age of the cluster (3.6 Myr, see above). Another effect, seldom considered, is described by Comerón et al. (2003): they propose that underluminous cluster members in the Lupus 3 dark cloud provide evidence for accretion-modified evolution. These objects, if not a priori excluded from a cluster sample, would appear older than the rest of the cluster members.

The effect that we consider to be the major cause of the large age spread we observe in the $`\sigma `$ Ori cluster is variability. We make use of 2MASS J-band magnitudes and I-band magnitudes from the literature. PMS stars are known to be variable in $`J,H`$ and $`K`$ (Carpenter et al., 2001; Carpenter, 2001) and we found evidence of K-band variability in this sample (Sect. 3.2). They found a dispersion among variable stars of about 0.1 mag in $`J`$; if the effect in the I-band is of similar amplitude and the variability in the two bands is independent, an uncertainty in the $`(IJ)`$ colour of more than 0.1 mag would not be surprising. If variability is related to the presence of a circumstellar disk, then the uncertainty could be even larger (Carpenter et al., 2001; Carpenter, 2001). Indeed, from the 6 objects we found to be variable in the K-band (Sect. 3.2), two objects have no I-band measurements; one object falls outside the model grid i.e. appears too young; of the remaining 3 objects, 2 have computed ages of 16 Myr and 20 Myr. Although limited by small number statistics, 3 out of 4 K-band variable objects have “anomalous” ages, compatible with the spread being caused by variability.

Therefore, we believe that most of the spread we observe in Fig. 1 is not a real age spread, and that the likely age of the cluster is larger than 3.6 Myr with an upper limit of 8 Myr imposed by Zapatero Osorio et al. (2002). We will take into account the age uncertainty in our subsequent analysis. This short summary of the problems in determining the age of young clusters highlights the importance of using consistent age determinations, for instance when trying to infer an overall disk destruction timescale (Sect. 5.3).

### 2.3 Reddening

The average reddening towards the O-star $`\sigma `$ Ori is quite low: $`E(BV)=0.05`$ (e.g. Brown et al., 1994). However, a superficial analysis (see below) of traditional colour-colour diagrams seems to hint at some amount of reddening towards cluster members. Therefore, we have decided to investigate the question of reddening for the target sample.

Fig. 2 shows the $`IJH`$ colour-colour diagram for the target sample. This diagram has been used with success to compute individual reddening to young stars in clusters (Lucas & Roche, 2000; Thompson, 2003). The solid lines are the empirical loci for main-sequence and giant stars (from Bessell & Brett (1988) converted to 2MASS magnitudes using the transformations from Carpenter et al. (2002)) and the dashed lines are the reddening bands (Rieke & Lebofsky, 1985). Traditionally, one attempts to de-redden individual objects with respect to main-sequence colours. But PMS stars have lower gravities than field dwarfs and thus populate the region between the main-sequence and giant loci. The dotted lines are isochrones from Baraffe et al. (1998) for 3 Myr and 5 Myr; the target colours in the diagram are quite well represented by the isochrones. According to I. Baraffe (private communication) we should be cautious when using H-band magnitudes produced by their models because of possible shortcomings in the water line-list used (the H-band includes several strong water absorption bands). Nevertheless, very few objects show any evidence of significant reddening. The few objects that seem reddened actually are found to exhibit a $`(KL^{})`$ excess that indicates the presence of a circumstellar disk (see Sect. 4). We could use the reddening measured towards $`\sigma `$ Ori, but it is so small that its effect is negligible — especially at IR wavelengths. Therefore, the results described in this paper were obtained without de-reddening the target colours and magnitudes.

## 3 Infrared Observations

### 3.1 K and L-band imaging

In order to compute meaningful $`(KL^{})`$ excesses, it is essential that we are able to detect the stellar photospheres both in the K- and L-bands. We used the 2MASS $`K_\mathrm{s}`$ magnitudes of the targets to estimate exposure times in the K-band and use the Baraffe et al. (1998) evolutionary models (see previous section) to estimate L’-band magnitudes. Every target was detected in both bands. The observations in the two filters were consecutive, in order to minimise the effect of variability in the measured $`(KL^{})`$ colours.

The observations were performed at the UKIRT on the 17th and 19th January 2003, with the newly installed imager-spectrometer UIST. We observed with the broad-band K and L filters on the Mauna Kea Observatory Near-Infrared (MKO-NIR) system (Tokunaga, Simons & Vacca, 2002). For both K- and L-band observations, we used the 0.06$`\mathrm{}`$ pixel scale and windowed the detector to 512 $`\times `$ 512 pixels to reduce overheads, resulting in a field of view of 31$`\mathrm{}\times 31\mathrm{}`$. Typically this results in no more than one cluster member per field, so that exposure times are tailored for each target. Adverse weather conditions meant that we only could observe for about 30% of the observing time allocated to this project; from our original target list we observed the 28 brightest (in $`K_\mathrm{s}`$) targets. On the first night conditions were stable and dry, with typical seeing of the order of 0.7 arcsecond (measured as the full-width at half-maximum of the K-band images); in the second night conditions were less stable, with seeing of about 1 arcsecond.

In the L-band, on-target individual exposure times ranged from 1.5 min to 15 min; the brightest targets were observed on a 4-point dither pattern and the fainter targets on an 8-point dither pattern. At each dither position, 27 (sky-limited) exposures of 0.81 s were co-added. An 8-point dither pattern is complete in approximately 3 min (on-target) and the cycle was repeated until the requested exposure time was reached for each target. The frames were reduced and combined using ORAC-DR (the UKIRT data reduction and high level instrument control software), with dedicated recipes for the reduction of thermal imaging data. Sky background was removed by subtracting a median-average of the closest-in-time exposures at the different dither positions; flat-field frames were created by combining normalised object frames using the median at each pixel. Aperture photometry was performed on the mosaiced frames using the Starlink package GAIA (Graphical Astronomy and Image Analysis Tool), with apertures of 1.2$`\mathrm{}`$ and 1.8$`\mathrm{}`$ for the first and second nights respectively. The statistical signal-to-noise ratio of the L-band photometry ranges from 50 to 9. To calibrate the observations we used the following standards (see Leggett et al., 2003): HD 40335 ($`L^{}=6.45`$ mag) and SAO 112626 ($`L^{}=8.56`$ mag), respectively for the 4-point and 8-point dither patterns. Standards were observed every one to two hours.

In the K-band the targets are rather bright (Table 1) so the exposure time was set to 10 s. These observations were performed in a 5-point dither pattern. The K-band images were also reduced and analysed within ORAC-DR and GAIA. These observations were calibrated using the following standards: FS 11 ($`K=11.252`$ mag) and FS 121 ($`K=11.302`$ mag) from the list of UKIRT Faint Standards. The signal-to-noise ratio of the K-band photometry exceeds 30 and is more precise than the corresponding $`L^{}`$ measurement in all cases.

### 3.2 K-band variability

For the targets in Table 1 we have two independent measurements of K-band magnitude: $`K_\mathrm{s}`$ from the 2MASS All-Sky Catalog of Point Sources (Cutri et al., 2003) — data taken in October 1998 — and the $`K`$ from the present observations. To convert the 2MASS photometry to the MKO system we first converted from 2MASS to the old UKIRT system using the transformation in Carpenter (2001) and then from this system to the new MKO system using the transformation provided in Hawarden et al. (2001).

Fig. 3 shows the histogram of the difference between the two sets of measurements in the MKO system. We can see that there is a systematic effect in the $`K`$ magnitude of the order of 0.04 mag, that is acceptable within the measurement uncertainties — it is not possible to say which of the steps of the transformation (2MASS to UKIRT or UKIRT to MKO) introduces this effect. Small differences can be easily explained by uncertainties in the calibrations and photometric system transformations. Nevertheless some objects show differences that clearly hint at K-band variability, namely targets 5, 6 and 8 and more marginally targets 3, 22 and 26. PMS stars are known to be variable at near-infrared wavelengths. Carpenter et al. (2001, 2002) found that $``$ 30% of objects in Orion A and Chamaeleon I molecular clouds are variable. For more than 70% of variable PMS stars, variability is characterised by small-amplitude ($``$0.1 mag), essentially colourless, fluctuations that can be attributed to rotational modulation by cool spots. However, when the observed fluctuations are of larger amplitude and occur both in magnitudes and colours, then a combination of effects related with the presence of circumstellar disks is a likely cause: hot spots, variable circumstellar extinction and varying disk properties (accretion rate and/or inner disk temperature). Objects 5, 6 and 8 could fall in this latter category (they also have large $`(KL^{})`$ excesses) and objects 3, 22 and 26 in the former but our data really does not allow us to say much more.

## 4 Analysis of the $`(KL^{})`$ colours of the cluster members

### 4.1 Colour-colour diagrams

The most immediate method to search for circumstellar disks is to look for stars with $`(HK)`$ or $`(KL)`$ excesses in the $`JHK`$ or $`JHKL`$ colour-colour diagrams, that would indicate the presence of circumstellar dust material. Near-IR colours trace the warmer dust and the magnitude of a near-infrared excess depends critically on the position and temperature of the inner disk boundary. Therefore, disk frequencies determined from K-band excesses tend to be lower limits to the true frequency. At longer wavelengths, IR excesses grow rapidly and disk properties interfere less with the ability to detect disks. In particular for the L-band, disks remain optically thick over a wide range of disk masses, resulting in measurable $`(KL)`$ excess down to very low disk masses (Wood et al., 2002). Therefore, a $`(KL)`$ excess is a very reliable disk diagnostic (Haisch et al., 2001a; Wood et al., 2002).

Fig. 4 shows several colour-colour diagrams for the 28 cluster members. As a clarification, throughout this paper $`(JH)`$ and $`(HK_\mathrm{s})`$ colours are always from 2MASS, while $`(KL^{})`$ is always from our UKIRT observations. The IRAS sources are known to possess disks so they will be excluded from the discussion on disk frequencies (so as not to introduce any bias), but they are shown here to illustrate the position in the diagrams of classical T Tauri stars. The diagrams show the position of the objects with respect to the main-sequence and giant star loci and their reddening bands (see Fig. 2), converted to the 2MASS and MKO systems using transformations from Bessell & Brett (1988), Hawarden et al. (2001) and Carpenter (2001). When analysing such diagrams, stars are considered to exhibit an excess when they appear to the right of the reddening vector for late-type stars. From the $`JHK_\mathrm{s}`$ diagram (top left) it can be seen that (with the exception of the IRAS sources) only one object exhibits a considerable K-band excess. This scenario changes dramatically when looking at the $`JHKL^{}`$ diagram (top right): maybe as many as 13/24 objects (or 54%) seem to have an L-band excess. Simply counting stars with excess in such diagrams may still underestimate the true disk frequency: the reddening bands are rather broad, meaning that early-type stars with an excess might not be distinguishable from later-type stars with small reddening and no excess.

The two diagrams at the bottom of Fig. 4 largely remove this ambiguity. In the $`HK_\mathrm{s}KL^{}`$ and $`IJKL^{}`$ diagrams, the reddening bands are narrow and almost parallel to the empirical loci, so that stars with an excess in the L-band are clearly above these bands. The Baraffe et al. (1998) isochrones for 3 Myr and 5 Myr predict which region of the diagrams would be populated by diskless PMS stars.

The $`IJKL^{}`$ diagram (bottom right) is particularly clear, with two populations of stars: one with colours in between the giant and main-sequence loci and one of objects with colours more than 0.1 mag above the empirical main-sequence locus. A legitimate question in the light of our concern about variability, is whether we can make use of $`(IJ)`$ colours that mix magnitudes obtained at different epochs ($`J`$ magnitudes are from 2MASS and $`I`$ magnitudes are from sources in the literature). The strength of this diagram, however, is that variability in $`(IJ)`$ would displace objects horizontally across the diagram (the same way as reddening) while an excess moves the objects almost vertically, so these effects act orthogonally. We can see from this diagram that, with respect to empirical main-sequence colours, 11/24 stars show no excess, 13/24 stars (or 54%) show excesses $`0.1`$ mag and 8/24 stars (or 33%) show excesses $`0.2`$ mag.

### 4.2 $`(KL^{})`$ excesses

PMS objects are young and have lower surface gravities than the field dwarfs used to calibrate the empirical relations used in the previous section (Luhman, 1999). We investigated if theoretical models (see Sect. 2.2) could be used to predict K- and L-band magnitudes for our sample. If, in the $`IJKL^{}`$ diagram (Fig. 4), we compare the colours of objects with no apparent excess with the theoretical isochrones of Baraffe et al. (1998), the model $`(KL^{})`$ colours seem to be too blue (even bluer than corresponding giant colours), especially for earlier spectral types ($`(IJ)<1.7`$ mag). In terms of surface gravity, we would expect PMS stars to have colours between main-sequence and giant stars, therefore we opt to use empirical spectral type calibrations to compute $`(KL^{})`$ excesses. Furthermore, a disk frequency based on the comparison with model colours would not be readily comparable to other surveys that tend to use the main-sequence calibration. If model colours were used, the $`IJKL^{}`$ diagrams implies that a larger fraction ($``$ 80%) of objects would have an $`(KL^{})`$ excess.

Spectral types are known for 18 objects in the sample (Zapatero Osorio et al., 2002, targets labelled ZO02 in Tables 1 & 2). Spectral types are also known for 2 of the IRAS sources. For the 6 objects identified by Kenyon et al. (2003) spectral types are not known. Their $`(IJ)`$ and $`(HK_\mathrm{s})`$ colours may be used to estimate spectral types, as none of these objects show considerable excess in the K-band. We use tabulated relations between spectral type and colours (again from Bessell & Brett 1988) to define linear relations between spectral type and each colour (these relations are only valid for spectral types later than M2). These are then interpolated and averaged to obtain the spectral type for the sources. The results are in column 2 in Table 2 (indicated by $``$). We apply the same procedure to the stars with known spectral types in order to test its validity. For the 10 sources with known spectral type later than M2, we obtain the following: for 6/10 stars we compute the correct spectral type within 0.5 subclass, for 3/10 within 1 subclass and for 1/10 within 1.5 subclass. So we are confident that our computed spectral types are mostly correct to within 1 subclass.

Even though these may not be the most accurate spectral type determinations, they are easily adequate for the purposes of computing $`(KL^{})`$ excesses. Between spectral types M2 and M6, one subclass is equivalent to at most a variation of 0.055 mag in photospheric $`(KL^{})`$ colour and 0.05 mag in $`(HK_\mathrm{s})`$. We adopt these values as estimates of the errors in the excesses introduced by this procedure; spectral types from the literature are uncertain by 0.5 subclass. For each target we have now spectral types and we have computed $`(HK_\mathrm{s})`$ and $`(KL^{})`$ excesses that are listed in columns 4 & 7 in Table 2. For the 2 IRAS sources without known spectral types, we compute lower limits for their excesses by assigning them spectral type M6 (even though their brightness clearly indicates earlier spectral type) and photospheric colours of respectively $`(KL^{})=0.46`$ mag and $`(HK_\mathrm{s})=0.40`$ mag.

In Fig. 5 we plot the computed excesses against K-band magnitude for both $`(KL^{})`$ and $`(HK_\mathrm{s})`$. The error bars are obtained by adding in quadrature the photometric uncertainty in the colours with the uncertainty introduced by the spectral type determination. As expected all 4 IRAS objects have very large excesses in the L-band, consistent with their mid-IR colours (Oliveira & van Loon, 2003). We consider that an object has a circumstellar disk if it exhibits an excess detected at a 2-$`\sigma `$ level. Thus 11/24 (or 46%) objects show a $`(KL^{})`$ excess. If we arbitrarily adjust for the systematic effect discussed on Sect. 3.2, then the computed excesses would be larger by about 0.04 mag and 14/24 objects would have a significant excess.

## 5 Discussion

### 5.1 How representative is this sample?

The sample we discuss here is not complete, but we believe it is representative. As we discuss below, there is no serious bias either for or against the detection of disks in this sample (apart from the 4 IRAS sources).

#### 5.1.1 X-ray selection bias

The equivalent width (EW) of the H$`\alpha `$ emission line has traditionally been used to classify T Tauri stars (TTS): classical T Tauri stars (CTTS) are defined as having large EW\[H$`\alpha `$\] as evidence of circumstellar accretion and thus of the presence of a circumstellar disk; weak-lined T Tauri stars (WTTS) have chromospheric EW\[H$`\alpha `$\] thus showing no signatures of accretion. CTTS are found to be underluminous in X-rays when compared with WTTS in the same star-forming regions (e.g. Neuhaüser et al., 1995; Flaccomio et al., 2003). This would imply that surveys for circumstellar disks in X-ray selected samples may be biased towards objects with no accretion disk signatures, depending on the sensitivity of the X-ray surveys.

But many PMS stars classified as WTTS based on their H$`\alpha `$ emission were actually found to have circumstellar disks. Haisch et al. (2001b) found that a large fraction of WTTS in IC 348 have considerable $`(KL)`$ excesses, i.e. WTTS are not necessarily naked TTS. Furthermore, Preibish & Zinnecker (2002) provided statistically significant evidence that accreting stars in IC 348 have lower X-ray luminosity when compared with non-accreting objects, but found no evidence for a difference between stars with or without $`(KL)`$ excess. This suggests that X-ray selection will not impair the detection of circumstellar disks in the L-band, even though it might be biased against strongly accreting objects. In fact, as we discuss in the next section, about half of the X-ray selected cluster members have large $`(KL^{})`$ excesses, similar to the overall disk fraction.

#### 5.1.2 Spectroscopic selection bias

Another potential source of bias has to do with the spectroscopic identifications of cluster members. The presence of the Li i 6708Å feature in stellar spectra is a sure sign of youth at least in the mass range we are considering here. However, photospheric spectral lines in CTTS can be heavily veiled at optical wavelengths (e.g. Gullbring et al., 1998). In particular, the hot continuum attributed to disk accretion can fill-in the Li i 6708 Å line (e.g. Magazzù, Rebolo & Pavlenko, 1992). The sample was selected on the basis of the detection of the lithium feature, therefore our estimated disk frequencies could be regarded as a lower limit. Kenyon et al. (2003) found a few objects in their sample of cluster candidates that show no lithium in their spectra but have radial velocities and Na i doublet strengths consistent with cluster membership. However, their measured H$`\alpha `$ line widths are not consistent with accreting PMS objects (White & Basri, 2003). No evidence has been found for heavily-veiled strongly-accreting cluster members, therefore we consider that the spectroscopic selection of cluster members is not a significant source of bias.

### 5.2 L-band excess versus other disk indicators

In this section we compare $`(KL^{})`$ excesses with other disk indicators. Zapatero Osorio et al. (2002) describe spectra of 18 of the stars in our sample (Table 1); in particular they have measured equivalent widths for the H$`\alpha `$ emission and several forbidden emission lines. A classification criterion for PMS stars based on EW\[H$`\alpha `$\] has to take into account the different levels of chromospheric H$`\alpha `$ emission observed for different spectral types (Martín, 1998). White & Basri (2003) have proposed the following (empirical and spectral-type dependent) classification: a PMS star is classified as a CTTS if EW\[H$`\alpha `$\]$`>`$3 Å for K0$``$K5 stars, EW\[H$`\alpha `$\]$`>`$10 Å for K7$``$M2.5 stars, EW\[H$`\alpha `$\]$`>`$20 Å for M3$``$M5.5 stars and EW\[H$`\alpha `$\]$`>`$40 Å for M6$``$M7.5 stars. This classification was devised using high-resolution spectra, so when applying it to intermediate- and low-resolution spectra (as is the case here) we have to take into account that EW\[H$`\alpha `$\] is probably overestimated due to blending with a nearby TiO feature. The presence of forbidden lines (e.g. \[O i\], \[N ii\] and \[S ii\]) in the spectrum of a PMS star is evidence of jets and outflows that are also related with accretion processes and circumstellar disks (e.g. Edwards et al., 1987).

In Fig. 6 (left), we plot EW\[H$`\alpha `$\] and spectral types for the objects from Zapatero Osorio et al. (2002); on the right we plot EW\[H$`\alpha `$\] against $`(KL^{})`$ excesses. We can see that 4 objects show EW\[H$`\alpha `$\] in excess of what is expected for their spectral types; these 4 objects have $`(KL^{})`$ excesses of the order of 0.17$``$0.98 mag, consistent with the presence of circumstellar disks (objects 3, 14, 17 and 19 in Tables 1 & 2). Of these 4 objects, 3 exhibit spectra with forbidden line emission (objects 3, 14 and 19). Of the 10 objects that were X-ray selected, only one object has a large EW\[H$`\alpha `$\]; half of these objects have a $`(KL^{})`$ excess. The IRAS sources have, as expected, large $`(KL^{})`$ excesses, as they also have large mid-IR excesses (Oliveira & van Loon, 2003).

There are also 3 objects that have $`(KL^{})`$ excesses larger than 0.4 mag, but do not show large EW\[H$`\alpha `$\]. This again illustrates that $`(KL^{})`$ is the most efficient and robust disk indicator. H$`\alpha `$ emission is known to be extremely variable (Guenther & Emmerson, 1997) and objects that have $`(KL^{})`$ excesses indicative of circumstellar disks could alternate from episodes of high accretion activity (CTTS state) to episodes of undetected accretion activity (WTTS). More interestingly, if we interpret this as an evolutionary sequence (i.e. the disk-star system evolves from an accreting system, through a weakly or non-accreting system, to a naked PMS star), this could be an indication that circumstellar disks can survive for some time after accretion onto the stellar surface has stopped. The $`\sigma `$ Ori cluster is at a crucial stage of disk destruction (see next section) and we would indeed expect to find a mixture of objects that are still accreting and others where accretion has strongly diminished or ceased.

### 5.3 Comparison with other young clusters

The most complete analysis of clusters’ disk frequencies is described in Haisch et al. (2001a, and references therein). They compile the results of several cluster surveys in the K- and L-bands. These authors have consistently used the same method to determine which objects possess a $`(KL)`$ excess: they place cluster members in a $`JHKL`$ colour-magnitude diagram — only those objects with $`K`$ magnitudes brighter than the completeness limit of the $`L`$ survey are considered, to assure that stellar photospheres are detected in both bands; they count the stars that lie to the right of the reddening vector that passes through the position in this diagram of an M5 main-sequence star. We can apply exactly the same procedure in Fig. 4 (top right): 13/24 stars (54 $`\pm `$ 15)%<sup>1</sup><sup>1</sup>1Statistical errors are conservatively estimated as $`\sqrt{N_{\mathrm{disk}}}/N`$. are to the right of the reddening band. This result depends on the adopted reddening law — the reddening law defines the position of the boundary between objects with and without excesses — but the different determinations are within the quoted statistical errors (Haisch et al., 2001a). Depending on the spectral type distribution of the sample, this method might underestimate the true disk frequency particularly for earlier spectral types (Lada et al., 2000).

According to this technique, the $`\sigma `$ Ori cluster has a (54 $`\pm `$ 15)% $`JHKL^{}`$-excess disk frequency, at an age of 3$``$8 Myr. How does this compare with other clusters? Such comparison can be done by placing the $`\sigma `$ Ori cluster in Fig. 1 from Haisch et al. (2001a) that plots the fraction of $`JHKL`$ excess objects against stellar age for several young clusters and associations. In terms of its age, the $`\sigma `$ Ori cluster can be compared with NGC 2264 (age $``$ 3.2 Myr, disk frequency 52% $`\pm `$ 4%) and NGC 2362 (age $``$ 5 Myr, disk frequency 12% $`\pm `$ 10%). Crucially, the measurements for these two clusters — particularly NGC 2362 — allowed the authors to estimate that the overall disk lifetime is about 6 Myr, with 50% of the disks dispersed by about 3 Myr. In this scenario, the $`\sigma `$ Ori cluster can play a key role to better constrain disk destruction timescales. In Haisch et al. (2001a) analysis, the ages of the clusters were determined using different PMS models and these authors estimate that this introduces an overall systematic uncertainty in the ages of the order of 1.2 Myr. The age of $`\sigma `$ Ori is also not well established (see discussion Sect. 2.2): 3$``$5 Myr is favoured by several authors (Béjar et al. 1999; Oliveira et al. 2002; Jayawardhana et al. 2003; this work) with an upper limit of 8 Myr (Zapatero Osorio et al., 2002). If the age of the cluster is $``$ 4 Myr then our estimate of disk frequency is consistent with the above mentioned timescales. However, if the cluster is older than 4 Myr, then our result suggests a slower disk destruction, i.e. a longer overall disk lifetime. The NGC 2264 and NGC 2362 surveys are for objects more massive than $``$ 0.85 M (i.e. earlier spectral types) while our sample populates the range 0.13$``$1.0 M. As mentioned above, this technique might be at fault for earlier spectral types. Furthermore, if disk-destruction timescales are mass-dependent then one has to be careful when comparing these results.

Lyo et al. (2003) analysed $`(KL)`$ excesses in the $``$ 9 Myr, sparsely-populated $`\eta `$ Chamaeleontis cluster. They found that of the 12 late-type stars in the central part of the cluster 7 objects have $`(KL)`$ excess. Their mass range is comparable to the mass range of our sample, but even so the disk frequency in $`\eta `$ Cha is remarkably high for its age. Environmental effects might explain such result: disk destruction timescales might be controlled also by the stellar density in the star-forming region and by the photoevaporation power of nearby O-stars. If this is the case, then $`\eta `$ Cha cannot be straight forwardly compared with the other clusters.

### 5.4 Comparison with disk model predictions

Wood et al. (2002) investigates the observational signatures of circumstellar disks in a simple evolutionary scenario by which disk mass (of small particles) decreases with time (for a stellar effective temperature of 4000 K, spectral type K7-M0). Except for high disk masses, the main contribution to the spectral energy distribution (SED) is from reprocessing of starlight (or disk irradiation), as inferred from observations (Hartmann et al., 1998). For the $`(KL)`$ excess (or $`\mathrm{\Delta }(KL)`$) their circumstellar disk model predicts that: i) large excesses ($`\mathrm{\Delta }(KL)`$ 0.7 mag) can only be achieved with the contribution of accretion luminosity from massive disks; ii) for passive disks, $`\mathrm{\Delta }(KL)`$ is insensitive to disk mass over several orders of magnitude; iii) for disk masses lower than 10<sup>-7</sup> M, $`\mathrm{\Delta }(KL)`$ decreases rapidly.

In our sample, 9 of the objects have spectral types K7$``$M0.5 and 3 of these objects have large excesses ($``$ 0.7 mag) that can hint at a contribution from accretion luminosity to the SEDs. Only one of these objects shows signatures of accretion in its spectrum (Sect. 5.2) but variability likely plays an important role. The objects with later spectral types tend to have modest excesses ($``$ 0.4 mag), consistent with disk irradiation.

## 6 Summary

In this paper we present new K- and L-band photometry for a representative sample of $`\sigma `$ Orionis cluster members. Different methods and disk-excess criteria will provide somewhat different disk frequencies. We computed a disk frequency for this cluster of (54$`\pm `$15)%, if measured directly from the $`JHKL^{}`$ colour-colour diagram, or (46$`\pm `$14)%, if excesses are computed with respect to predicted photospheric colours (according to the objects spectral types, 2-$`\sigma `$ excess detections). Therefore, when trying to determine a disk destruction timescale, one has to be consistent in the methods used.

An overall disk lifetime of 6 Myr has been proposed, but this relies heavily on the age and disk frequency derived from NGC 2362 (age $``$ 5 Myr and disk-frequency $``$ 12%, Haisch et al. 2001a). The disk frequency for the $`\sigma `$ Ori cluster could increase support for this timescale, but we find that the crucial factor is not so much the derived disk frequency, but the age determination of the cluster. As is the case for most young clusters (Hartmann, 2001), the age of the $`\sigma `$ Ori cluster is uncertain, with a likely age in the range 3$``$5 Myr (Oliveira et al., 2002; Zapatero Osorio et al., 2002; Jayawardhana et al., 2003). If the cluster age is about 3 Myr then it lends support to the above mentioned disk lifetimes. However, if the cluster is older than 4 Myr then it points towards slower disk destruction timescales.

If disk dissipation timescales are mass dependent then the disk frequencies derived for the different clusters are not so readily comparable. Another factor that is probably relevant is the environment of the star forming region, in particular the rates of close stellar encounters and the amount of ionizing radiation produced by massive (O-type) cluster members.

## Acknowledgments

We thank the staff of the United Kingdom Infrared Telescope (UKIRT) for their support during the observing run. The UKIRT is operated by the Joint Astronomy Centre on behalf of the U.K. Particle Physics and Astronomy Research Council (PPARC). This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. JMO acknowledges the financial support of PPARC.


text:
# Formation of soliton trains in Bose-Einstein condensates as a nonlinear Fresnel diffraction of matter waves

## Abstract

The problem of generation of atomic soliton trains in elongated Bose-Einstein condensates is considered in framework of Whitham theory of modulations of nonlinear waves. Complete analytical solution is presented for the case when the initial density distribution has sharp enough boundaries. In this case the process of soliton train formation can be viewed as a nonlinear Fresnel diffraction of matter waves. Theoretical predictions are compared with results of numerical simulations of one- and three-dimensional Gross-Pitaevskii equation and with experimental data on formation of Bose-Einstein bright solitons in cigar-shaped traps.

Discovery of Bose-Einstein condensate (BEC)  has created new active field of research of quantum macroscopical behavior of matter. Among most spectacular evidences of such macroscopic behavior one can mention formation of interference fringes between two condensates  and creation of dark  and bright  solitons. The interference phenomenon is usually considered in framework of a linear wave theory, whereas solitons are treated as a nonlinear wave effect. At the same time, basically, these two phenomena have much in common. For example, formation of bright soliton trains in nonlinear wave systems is often explained as a result of modulational instability, where selection of the most unstable mode is a result of interplay of interference and nonlinear effects (see, e.g. ). Such interconnection of interference and soliton phenomena is demonstrated most spectacularly in formation of solitons in vicinity of a sharp edge of density distribution. In this case, at linear stage of evolution the linear diffraction provides an initial modulation of the wave and further combined action of interference and nonlinear effects leads to formation of soliton trains. Without nonlinear effects, such kind of time evolution of a sharp wave front would be a temporal counterpart of usual spatial Fresnel diffraction and therefore soliton train formation at the sharp front of nonlinear wave can be called a nonlinear Fresnel diffraction.

Similar formation of oscillatory structures at sharp wave front or after wave breaking in modulationally stable systems described by the Korteweg-de Vries equation is well known as a “dissipationless shock wave” (see, e.g. ). Its theoretical description is given  in framework of Whitham theory of nonlinear wave modulations , where the oscillatory structure is presented as a modulated nonlinear periodic wave which parameters change little in one wavelength and one period. Then slow evolution of the parameters of the wave is governed by Whitham equations obtained by averaging of initial nonlinear wave equations over fast oscillations of the wave. Application of this method to modulationally unstable systems has been given for important particular case of soliton train formation at the sharp front of a long step-like initial pulse . Here we shall consider by this method formation of solitons in BEC with negative scattering length (attractive interaction of atoms).

We suppose that condensate is confined in a very elongated cigar-shaped trap whose axial frequency $`\omega _z`$ is much less than the radial frequency $`\omega _{}`$. In the first approximation we can neglect the axial trap potential and suppose that condensate is contained in a cylindrical trap ($`\omega _z=0`$) and its initial density distribution has a rectangular form. Evolution of BEC is governed by three-dimensional (3D) Gross-Pitaevskii (GP) equation

$$i\mathrm{}\psi _t=\frac{\mathrm{}^2}{2m_a}\mathrm{\Delta }\psi +\frac{1}{2}m_a\omega _{}^2(x^2+y^2)\psi +g|\psi |^2\psi ,$$
(1)
for the condensate wave function $`\psi `$, where we use standard notation $`g=4\pi \mathrm{}^2a_s/m_a`$ for the effective nonlinear coupling constant, $`a_s<0`$ is the $`s`$-wave scattering length, and $`\psi `$ is normalized on the number of particles in BEC, $`|\psi |^2𝑑𝐫=N.`$ For analytical treatment of nonlinear Fresnel diffraction it is important to determine conditions when the 3D equation (1) can be reduced to its one-dimensional (1D) approximation (see, e.g. )

$$i\mathrm{}\mathrm{\Psi }_t=\frac{\mathrm{}^2}{2m_a}\mathrm{\Psi }_{zz}+g_{1D}|\mathrm{\Psi }|^2\mathrm{\Psi },|\mathrm{\Psi }|^2𝑑z=N,$$
(2)
where

$$g_{1D}=\frac{g}{2\pi a_{}^2}=\frac{2\mathrm{}^2a_s}{m_aa_{}^2},a_{}=\sqrt{\frac{\mathrm{}}{m_a\omega _{}}},$$
(3)
that is the transversal degrees of freedom are frozen. It is well known (see, e.g. ) that a homogeneous distribution with linear density $`n_0=|\mathrm{\Psi }|^2=\mathrm{const}`$ described by (2) with negative $`g_{1D}`$ ($`a_s<0`$) is unstable with respect to self-modulation with increment of instability equal in our present notation to

$$\mathrm{\Gamma }=\frac{\mathrm{}K}{2m_aa_{}}\sqrt{8|a_s|n_0(a_{}K)^2},$$
(4)
where $`K`$ is a wavenumber of small periodic modulation. The most unstable mode has the wavenumber

$$K_{max}=2\sqrt{|a_s|n_0}/a_{}$$
(5)
and the corresponding increment is equal to

$$\mathrm{\Gamma }_{max}=4|a_s|n_0\omega _{}.$$
(6)
This means that after time $`1/(|a_s|n_0\omega _{})`$ the homogeneous condensate splits into separate solitons (diffraction fringes) and each soliton (diffraction fringe) contains about $`N_sn_0/K_{max}`$ atoms. If in 3D GP equation (1) the nonlinear energy $`gN_sK_{max}/a_{}^2gn_0/a_{}^2`$ in each solitons is much less than the kinetic energy in the transverse direction, $`\mathrm{}^2/m_aa_{}^2`$, then the transverse motion is reduced to the ground state oscillations and the 3D condensate wave function can be factorized into $`\psi =\varphi _0(x,y)\mathrm{\Psi }(z,t)`$, where $`\varphi _0=(\sqrt{\pi }a_{})^1\mathrm{exp}[(x^2+y^2)/(2a_{}^2)]`$ is the ground state wave function of transverse motion, and $`\mathrm{\Psi }(z,t)`$ obeys to the effective 1D nonlinear Schrödinger (NLS) equation (2). Thus, the condition of applicability of 1D equation (2) for description of solitons formation is

$$n_0|a_s|1,$$
(7)
which means that the instability wavelength $`1/K_{max}`$ is much greater than the transverse radius $`a_{}`$ of BEC. If (7) is not satisfied, then the transverse motion has to be taken into account which may lead to collapse of BEC inside each separate soliton. Therefore we shall confine ourselves to the BEC described by the 1D NLS equation under supposition that the initial distribution satisfies the condition (7).

To simplify formulae in the analytic theory, we transform (2) to dimensionless variables $`\tau =2(|a_s|n_0)^2\omega _{}t,`$ $`\zeta =2|a_s|n_0z/a_{},`$ $`\mathrm{\Psi }=\sqrt{2|a_s|}n_0u,`$ so that (2) takes the form

$$iu_\tau +u_{\zeta \zeta }+2|u|^2u=0,$$
(8)
and $`u`$ is normalized to the effective length $`L`$ of the condensate $`|u|^2𝑑\zeta =L/a_{}`$ measured in units of $`a_{}`$. We are interested in the process of formation of solitons (nonlinear Fresnel diffraction fringes) at the sharp boundary of initially rectangular distribution. Since this process takes place symmetrically at both sides of the rectangular distribution, we can confine ourselves to the study of only one boundary. This limitation remains correct until the nonlinear waves propagating inside the condensate collide in its center. If the initial distribution is long enough, this time is much greater than the time of solitons formation. Thus, we consider the initial distribution in the form

$$u(\zeta ,0)=\{\begin{array}{c}\gamma \mathrm{exp}(2i\alpha \zeta ),\mathrm{for}\zeta <0\hfill \\ 0,\mathrm{for}\zeta >0,\hfill \end{array}$$
(9)
where $`\gamma `$ is the height of initial step-like distribution and $`\alpha `$ characterizes the initial homogeneous phase. The problem of this kind has already been considered in some other problems of nonlinear physics  and we shall present here only the main results.

Due to dispersion effects described by the second term in Eq. (8), the sharp front transforms into slightly modulated wave which describes usual Fresnel diffraction of atoms. In our case the diffraction pattern evolves with time rather than is “projected” on the observation plane. The linear stage of evolution is followed by the nonlinear one in which combined action of dispersion and nonlinear terms yields the pattern which can be represented as a modulated nonlinear periodic wave or, in other words, a soliton train. We suppose that this soliton train contains large enough number of solitons, so that their parameters change little in one wavelength and one period. Then, in framework of Whitham theory, the density of BEC can be approximated by a modulated periodic solution of Eq. (8) (see )

$$n=|u(\zeta ,\tau )|^2=(\gamma +\delta )^24\gamma \delta \text{sn}^2(\sqrt{(\alpha \beta )^2+(\gamma +\delta )^2}\theta ,m),$$
(10)
where $`\text{sn}(x,m)`$ is the Jacobi elliptic function,

$`\theta =\zeta V\tau ,V=2(\alpha +\beta ),`$ (11)
$`m=4\gamma \delta /[(\alpha \beta )^2+(\gamma +\delta )^2],`$ (12)
the parameters $`\alpha `$ and $`\gamma `$ are determined by the initial condition (9), and $`\beta `$ and $`\delta `$ are slow functions of $`\zeta `$ and $`\tau `$. Their evolution is governed by the Whitham equation

$$\frac{(\beta +i\delta )}{\tau }+v(\beta ,\delta )\frac{(\beta +i\delta )}{\zeta }=0,$$
(13)
where Whitham velocity $`v(\beta ,\delta )`$ is given by the expression

$$v(\beta ,\delta )=2(\alpha +\beta )\frac{4\delta [\gamma \delta +i(\beta \alpha )]\mathrm{K}}{(\beta \alpha )(\mathrm{K}\mathrm{E})+i[(\delta \gamma )\mathrm{K}+(\delta +\gamma )\mathrm{E}]},$$
(14)
$`\mathrm{K}=\mathrm{K}(m)`$ and $`\mathrm{E}=\mathrm{E}(m)`$ being the complete elliptic integrals of the first and second kind, respectively. Since our initial condition (9) does not contain any parameters with dimension of length, the parameters $`\beta `$ and $`\delta `$ can only depend on the self-similar variable $`\xi =\zeta /\tau `$. Then Eq. (13) has the solution

$$\zeta /\tau =\xi =v(\beta ,\delta )$$
(15)
with $`v(\beta ,\delta )`$ given by (14). Separation of real and imaginary parts yields the formulae

$`\zeta /\tau =4\beta 2(\gamma ^2\delta ^2)/(\beta \alpha ),`$ (16)
$`{\displaystyle \frac{(\alpha \beta )^2+(\gamma \delta )^2}{(\alpha \beta )^2+\gamma ^2\delta ^2}}={\displaystyle \frac{\mathrm{E}(m)}{\mathrm{K}(m)}},`$ (17)
which together with Eq. (12) determine implicitly dependence of $`\beta `$ and $`\delta `$ on $`\xi =\zeta /\tau `$. It is convenient to express this dependence in parametric form

$`\beta (m)=\alpha \gamma \sqrt{4A(m)(1+mA(m))^2},`$ (18)
$`\delta (m)=\gamma mA(m),`$ (19)
where

$$A(m)=\frac{(2m)\mathrm{E}(m)2(1m)\mathrm{K}(m)}{m^2\mathrm{E}(m)}.$$
(20)
Substitution of these expressions into (10),(11) yields the density $`n`$ as a function of $`m`$. Since the space coordinate $`\zeta `$ defined by Eq. (16) is also a function of $`m`$ at given moment $`\tau `$, we arrive at presentation of dependence of $`n`$ on $`\zeta `$ in parametric form. The limit $`m0`$ corresponds to a vanishing modulation, and this edge point moves inside the condensate according to the law

$$\zeta _{}=(4\alpha +4\sqrt{2}\gamma )\tau .$$
(21)
The other edge with $`m1`$ moves according to the law

$$\zeta _+=4\alpha \tau ,$$
(22)
and corresponds to the bright solitons (or fringes of nonlinear diffraction pattern) at the moment $`\tau `$. The whole region $`\zeta _{}<\zeta <\zeta _+`$ describes the oscillatory pattern arising due to nonlinear Fresnel diffraction of the BEC with initially sharp boundary at $`\zeta =0`$.

We have performed numerical simulation of 1D and 3D GP equations with the aim to compare approximate Whitham theory with numerical results. The 1D density distributions calculated numerically from (8) and analytically are shown in Fig. 1. We see excellent agreement between the theoretical and numerical predictions of the height of the first soliton generated from initially step-like pulse, but its position given by analytical formula is slightly shifted with respect to numerical result. This is well-known feature of asymptotic Whitham approach  which accuracy in prediction of location of the oscillatory pattern cannot be much better than one wavelength. Thus, we see that the above theory reproduces the numerical results quite well for period of time $`\tau 2`$. For much greater time values some other unstable modes different from one-phase periodic solution (10) can also give considerable contribution into wave pattern. Nevertheless, the qualitative picture of soliton pattern remains the same.

For 3D numerical simulation, the GP equation (1) was transformed to dimensionless form by means of substitutions $`x=a_{}x^{},`$ $`y=a_{}y^{},`$ $`z=a_{}z^{},`$ $`t=2t^{}/\omega _{},`$ $`\psi ^{}=(N^{1/2}/a_{}^{3/2})\psi ,`$ so that it takes the form

$$i\psi _t=\mathrm{\Delta }\psi +r^2\psi (8\pi N|a_s|/a_{})|\psi |^2\psi ,$$
(23)
where primes are omitted for convenience of the notation and $`|\psi |^22\pi r𝑑r𝑑z=1`$, $`r^2=x^2+y^2`$. Evolution of the density distribution $`\rho (z)=_0^{\mathrm{}}|\psi (r,z)|^22\pi r𝑑r`$ along the axial direction is shown in Fig. 2 for the values of the parameters corresponding to the experiment  ($`a_s=3a_0,`$ $`\omega _{}=2\pi 625\mathrm{Hz},`$ $`L=300a_{}`$) except for the number of atoms which was chosen to be $`N=510^3`$ in order to satisfy the condition (7), so that $`|a_s|n_0=1.710^3`$. We see that diffraction (soliton) pattern arises after the dimensionless time $`t400`$ which corresponds after appropriate scaling transformation to $`\tau 2`$ in Fig. 1. The width of solitons in Fig. 2 also agrees with the width predicted by 1D analytical theory and numerics. The spatial distribution of the condensate density $`|\psi (r,z)|^2`$ is illustrated by Fig. 3. The 3D nonlinear interference pattern is clearly seen. For greater values of the condensate density, when 1D theory does not apply, numerical simulation demonstrates similar evolution of the diffraction pattern up to the moment when collapse starts in each separate soliton. Thus, formation of solitons in the experiment  with large initial number of atoms $`N10^5`$ goes through collapses with loss of atoms until the remaining atoms can form stable separate soliton-like condensates. The present theory emphasizes the importance of the initial stage of evolution with formation of the nonlinear Fresnel diffraction pattern.

Formation of soliton trains in BEC confined in a cigar-shaped trap has also been studied numerically in . The results of 1D simulation in  agree qualitatively with our results. In numerics of  strong losses were introduces to prevent fast collapse of BEC with large number of atoms. Nevertheless, formation of soliton trains was also observed.

The above theory is correct for evolution time much less than period of oscillations $`2\pi /\omega _z`$ in the axial trap. When the axial trap is taken into account, solitons acquire velocities in axial direction even if initial phase is equal to zero. The number of solitons produced ultimately from some finite initial BEC distribution can be found by means of quasiclassical method applied to an auxiliary spectral problem associated with the NLS equation (8) in framework of the inverse scattering transform method . If the initial wave function is represented in the form $`u_0(\zeta )=\sqrt{n_0(\zeta )}\mathrm{exp}(i\varphi _0(\zeta ))`$, then the total number of solitons is equal approximately to

$$N_s=\frac{1}{\pi }\sqrt{n_0(\zeta )+\frac{v_0^2(\zeta )}{4}}𝑑\zeta \frac{1}{2},$$
(24)
where $`v_0(\zeta )=\varphi _0(\zeta )/\zeta `$ is the initial velocity distribution of BEC. If there is no initial phase imprinted in BEC, then the total number of solitons is given by the formula

$$N_s=(\sqrt{2|a_s|}/\pi a_{})|\mathrm{\Psi }|𝑑z,$$
(25)
which is written in dimensional units and we have neglected a “one-half” term in (24).

In experiment, the initial stage is usually obtained by sudden change of the sign of the scattering length from positive to negative one, so that initial density distribution has, for large enough number of atoms, the Thomas-Fermi form

$$|\mathrm{\Psi }|^2=(3N/4Z)\left(1z^2/Z^2\right),$$
(26)
where $`Z`$ is the Thomas-Fermi half-length of the condensate. Then substitution of (26) into (25) gives

$$N_s=\sqrt{3N|a_s|L}/(4a_{}),$$
(27)
where $`L=2Z`$ is the total length of the condensate. Up to constant factor, this estimate coincides with one obtained in  by division of $`L`$ by the instability wavelength $`1/K_{max}`$. Note that this estimate includes also very small solitons which cannot be observed in real experiments, so that it must be considered as an upper limit of the number solitons which can be produced from a given initial distribution. The same property of this kind of estimate for number of dark solitons has been observed in comparison of analytical theory with numerical simulations in . The axial potential influences mainly on velocities of solitons, so the above estimate can be applied to the condensate in a cigar-shape trap under condition that inequality (7) is fulfilled.

In conclusion, we have studied theoretically and numerically the process of formation of soliton trains near the sharp edges of the density distribution of BEC. The arising oscillatory regions can be considered as nonlinear Fresnel diffraction fringes of matter waves.

This work was supported by FAPESP (Brazil) and CNPq (Brazil). A.M.K. thanks also RFBR (grant 01–01–00696) for partial support. A.G. would like to thank Randy Hulet and his group for useful discussions.