text:
# A New Method for ISOCAM Data Reduction – I. Application to the European Large Area ISO Survey Southern Field: Method and Results

## 1 Introduction

The Infrared Space Observatory (ISO; Kessler et al. 1996) was the successor to the Infrared Astronomical Satellite (IRAS). ISO, besides carrying out detailed studies of individual objects and small regions, has provided an opportunity to perform survey work at sensitivities of several orders of magnitude better than its precursor. Thus, a significant fraction of the mission time was spent on field surveys. The largest survey conducted with ISO is the European Large Area ISO Survey (ELAIS), which provides a link between the IRAS survey and the deeper ISO surveys. ELAIS is a collaboration between 20 European institutes which involves a deep, wide-angle survey at high Galactic latitudes, at wavelengths of 6.7 $`\mu `$m, 15 $`\mu `$m, 90 $`\mu `$m and 175 $`\mu `$m with ISO (see Oliver et al. 1997 and Oliver et al. 2000 for a detailed description of the Survey). In particular, the 15 $`\mu `$m survey was carried out with the ISO-CAM camera (Cesarsky et al. 1996) over a total area of $``$ 13 deg2, divided into 4 main fields and several smaller areas. One of the main fields, S1, and one of the smaller areas, S2, are located in the southern hemisphere. S1 is centered at $`\alpha `$(2000) = 00h 34m 44.4s, $`\delta `$(2000) = -43 28 12′′ and covers an area of $`2^{}\times 2^{}`$, while S2 is centered at $`\alpha `$(2000) = 05h 02m 24.5s, $`\delta `$(2000) = -30 36 00′′ and covers an area of $`21^{}\times 21^{}`$. The whole S1 and S2 areas have been surveyed in the radio (at 1.4 GHz, Gruppioni et al. 1999; Gruppioni, Ciliegi, Oliver et al. 2000 in preparation) in several optical bands and in the near-infrared (La Franca et al. 2000 in prep.; Heraudeau et al. 2000, in prep.).

Since ELAIS is the largest survey performed by ISO and covers just the gap in flux density that exists between the IRAS Survey and the ISOCAM Deep and Ultra-Deep Surveys (Elbaz et al. 1999), it was extremely important to obtain the best and most reliable possible results from these data through an accurate data reduction.

To this purpose, we have developed a new ISOCAM data reduction technique (the LARI technique) especially designed for the detection of faint sources. This method, designed by C. Lari and based on the assumption of the existence of two different time scales in ISOCAM transients, has been tested on ISOCAM-HDF data, providing excellent results in agreement with those obtained with the PRETI technique (Starck et al. 1999).

Before attempting to reduce the entire ELAIS survey, we decided to apply the LARI technique to one single field, in order to test the capabilities of our method and to adapt some of its tasks for this specific set of data. In particular, we have applied the LARI technique to the 15 $`\mu `$m data in the southern ELAIS field S1, where most of the available multi-wavelength follow-up observations are available. Here we present the results of the LARI method in S1, as well as the complete 15 $`\mu `$m source catalogue obtained with this technique.

The paper is structured as follows: in section 2 we present the survey strategy and parameters; in section 3 we give a detailed description of the new data reduction technique that we have developed and used; in section 4 we describe the reduction and analysis of our data; in section 5 we present the results of tests made on simulated data; in sections 6, 7 and 8 we discuss the source photometry the calibration accuracy and the astrometric corrections respectively, while in section 9 we describe our source catalogue and in section 10 we present our conclusions.

The source counts obtained from these data will be presented and discussed in a companion paper (Gruppioni, Lari, Pozzi et al. 2000, Paper II).

## 2 The ELAIS survey observation strategy

The S1 field, as well as the other ELAIS survey areas, was selected for its high Ecliptic latitude ($`|\beta |>40^{}`$, to reduce the impact of Zodiacal dust emission), for its low cirrus emission ($`I_{100\mu m}<1.5`$ MJy/sr) and for the absence of any bright IRAS 12 $`\mu `$m sources ($`S_{12\mu m}>0.6`$ Jy). In figure 1 the location of the S1 survey field is shown, overlaid on a Cirrus map (the COBE normalized IRAS maps of Schlegel et al. 1998). IRAS sources with 12 $`\mu `$m fluxes brighter than 0.6 Jy are also plotted.

The ELAIS ISOCAM survey was conducted in raster mode with the LW2 (6.7 $`\mu `$m) and LW3 (15 $`\mu `$m) filters. The ISOCAM detector was stepped across the sky in a grid pattern, with about half detector width steps in one direction and the whole detector width steps in the other. In this way, the reliability was improved as each sky position was observed twice in successive pointings and the overheads were reduced because each raster covered a relatively large area (40$`{}_{}{}^{}\times 40^{}`$). At each raster pointing (i.e. grid position of the raster) the 32$`\times `$32 ISOCAM detector was read out several times. Table 1 describes the observation parameters for the LW3 filter.

## 3 Lari Technique: Generalities

As already described in detail by Starck et al (1998), ISOCAM data obtained with the long wavelength detector (LW) are affected by several problems. The two main effects, which become more important the deeper we push for source detection, are produced by cosmic ray impacts (‘glitches’) and transient behaviour (slow response of the detector to flux variations).

Usually, ‘glitches’ can be divided into three categories: common, faders and dippers, according to their behaviour, decay time and influence on the pixel responsivity. Slow decreases of the signal following cosmic ray impacts are called faders, while depletions in the detector gain, followed by a reduction of the pixel sensitivity very slowly recovering afterwards (see figure 2, $`top`$ panel), are called dippers. These two effects are believed to be associated with proton or $`\alpha `$ particle impacts on the detector, while cosmic ray electrons produce common ‘glitches’. Common glitches only last one readout and their decay time is relatively fast (lasting only a few readouts), while faders and dippers have much longer lasting impact on the pixel sensitivities. So, the number of frames affected by the latter is much higher than in the case of common glitches, the sensitivity of pixels taking from tens to hundreds of seconds to recover completely. However, common glitches are much more frequent than faders and dippers and, if not correctly removed, they may look like sources on the maps and produce false detections. For this reason, the data cleaning is an extremely delicate process, which requires great care in order to produce highly reliable final maps and source lists.

The LARI method was mainly developed to overcome the main problems affecting ISOCAM LW data and to give better quality maps and as complete and reliable source catalogues as possible. Analogously to the PRETI method (Aussel et al. 1999), our algorithm corrects the cube of ISOCAM data for cosmic rays and transient effects before reconstructing the images and carrying out source detection.

The model on which the LARI method is based (described in detail in Appendix A), rests on the assumption that the incoming flux of charged particles generates transient behaviour producing two different time scale effects: a fast (breve) and a slow (lunga) one. The latter component accounts for the slow response of the detector and is essential in recovering the transient effects of the dippers. Each of the two time scales is associated with an independent reservoir of charge, which decays with this characteristic time-scale towards the contacts (i.e. a multi-component model for semiconductors). These two reservoirs of charge are fed by both incident infrared photons and cosmic rays. The latter are also able to trigger a fast charge release towards the contacts (‘glitches’). When a cosmic ray particle hits the detector, the quickly-varying charge reservoir breve is on average increased, while the slowly decaying charge reservoir is quickly forced to release part of its charge content. Thus, while the breve component is fed by a large fraction of the incident photons (around 40-45 %), the lunga one is fed only by a few percent of them. The remaining fraction is very quickly forced towards the contacts (prompt component). Due to differences between the two time scales of about a factor of 20 when the process reaches stabilisation, the lunga component collects a higher total amount of charge with respect to the breve one.

The value of both time constants depends on the signal level (which is fixed by observations) such that the lower is the signal, the larger is the time constant. Our model simply assumes the time constant to be inversely proportional to the amount of accumulated charge.

To first order, faders are described in this model as discontinuities mainly in the breve charge reservoir, caused by the cosmic ray impact. Similarly, to first order the dippers are discontinuities mainly in the lunga charge reservoir. The maximum depth of the dippers is determined by the fraction of the flux feeding the lunga reservoir. The overwhelmingly large majority of pixels are well-fit by a lunga fraction of $`0.1`$, implying that dippers cannot exceed one-tenth of the sky background level. Very occasionally however some dippers exceed this threshold. To account for these, an additional zero-point dark current ‘offset’ can be set, so the maximum depth is not larger than one-tenth of the revised total background. Incidentally, the presence of dippers in the dark current records (that have zero background) shows that this ‘offset’ is a general property of the detector, almost certainly fed by the thermal noise.

## 4 Application to ISOCAM LW ELAIS data

The application of our model to the ISOCAM LW data obtained in the ELAIS fields required some particular adaptation of the algorithms and the construction of some ‘ad hoc’ procedures necessary to overcome the specific problems generated by the chosen observational strategy.

The main cause of problems in the ELAIS data is the very short integration time. In fact, the time spent by the detector on each readout in these observations is only 2 seconds (see Table 1) and the total time spent on each raster pointing is 10$`\times `$ the integration time. The short observing time over each raster position (10$`\times `$2 sec) not only affects the signal-to-noise ratio, but it has two major negative effects on the data:

* since it is shorter than the fast time scale, it makes strong glitches hide real sources;
* only a fraction (60%) of the total incoming flux is recorded during each exposure, thus causing large photometric errors.

In our model, glitches are treated as discontinuities in the charge ($`Q`$) reservoir, with constant parameters $`a`$ and $`e`$ (see Appendix A). However, immediately after the maximum of a glitch, the detector is considered to behave normally under a constant (over the raster pointing) flux $`I`$. This is not completely correct for very strong glitches, which may cause the signal immediately following the maximum to be higher than predicted and this is mostly true for short integration time observations (i.e. 2 seconds like ELAIS).

Moreover, the relation between the increment of the breve component and the decrement of the lunga one is not constant. In fact, cosmic rays producing a higher increment in the breve reservoir than in the lunga one look like faders, while in the opposite case we have dippers. In the ELAIS data there are dippers without an initial glitch spike: generally we find the glitch feature in a contiguous pixel but, very rarely, it is completely absent.

Another problem arises from the fact that the Point Spread Function (PSF) is spatially under-sampled in all ISOCAM LW3 observations with pixel-field-of-view (PFOV) lens greater than 1.5 arcsec (in our case PFOV = 6 arcsec). Thus, any position determination method applied to individual point sources gives biased results for under-sampled data (the worse sampled the data, the more the resulting position is centered on a pixel). This bias can be corrected to some extent and the source position can be improved up to a fraction of the pixel size by taking into account the a priori raster pattern. In any case, the PSF is not unique for all sources, but it depends strongly on the source location within a pixel. This PSF, corresponding to the position of a source in the raster map and with an average FWHM of $``$10\`\`, will be referred as “effective PSF” throughout the paper. Moreover, in the ELAIS data each pixel in the final raster map comes from the combination and projection on the sky of different overlapping single images. For this reason, a source in the raster map is produced by the combination of different source images, where the source has different flux distributions, depending on its location within the pixel of the single images. This is a serious problem which affects source detection, source photometry and the completeness of the catalogue. In this work, we have carefully analysed and tried to quantify this combination effect by performing simulations (see section 5).

### 4.1 ELAIS Data Reduction with the Lari Method

All the codes developed for data reduction with the LARI method are written in the Interactive Data Language (IDL) and the whole data reduction and analysis has been performed with IDL software. The main process of ISOCAM LW data reduction with the LARI method consists of several basic steps. First, the raw data are converted into a raster structure containing all the information about the observation (pointings, instrument configuration, etc.) as well as the single images (one for each pointing, which are combined together to give the final raster image). The images are then converted to ADU/gain/s and the dark current subtracted. These first two steps are performed using the CAM Interactive Analysis (CIA) package. Next, the data are corrected for short time cosmic rays and the affected readouts are masked before copying the “de-glitched” data to a new structure (called “liscio”). This structure contains also the initial guess for the parameters defining the lunga and breve reservoirs of charge and information about the main ‘glitches’ and ‘dippers’ (derived by the deglitching process). The task which performs the first guess for parameters, also evaluates the background and the minimum ‘offset’ to be added to the data in order to have the ‘dipper’ depth in the range allowed by the model (one tenth of the sky background, see section 3). This is done for each pixel. The code not only finds the stabilization background level, that is, the zero level for data fully recovered from transients, but models the ‘glitches’, the sources and the background with all the transients over the whole pixel history. Moreover, our code is able to predict the trend we would have on each raster position if only the stabilization background flux was hitting the detector (starting from the previously accumulated charges, i.e. the ‘local background’). The excess with respect to this ‘local background’ represents the flux excess not recovered from transients. The maps of this excess, after flat fielding, are called ‘unreconstructed’ maps and (in case of a good enough fit) represent the effective flux collected by the detector during the raster exposure. In this paper, fluxes obtained from ‘unreconstructed’ maps will be named fs, while fluxes measured on ‘reconstructed’ maps (i.e. reconstructed from transient effects) will be named fsr. These two fluxes that we can measure for a source are shown in figure 2($`bottom`$ panel): the dot-dashed line represents the ‘unreconstructed’ data, while the dashed line represents the data ‘reconstructed’ for transients.

With our code, we created a model data-set for the deglitched data, reproducing not only the source signal, but also all the transient effects affecting the data. In figure 2 an example of pixel history is shown, together with the background and data models obtained with our algorithm.

In outline, the fitting algorithm starts with the brightest glitches in the raster, assumes discontinuities at these positions, and tries to find a fit to the time-lines that satisfies the solid-state physics of the detector. If no acceptable fit is found, the next fainter glitch is considered as a potential discontinuity, and so on. Because of the reduced number of useful readouts in the ELAIS raster data, in the fit we use fixed default values for all the pixels (the physical parameters scaled only for the background level), leaving as free parameters only the charge values at the beginning of the observations and at the top of glitches.

By successive iterations, the parameters and the background for each pixel are adjusted to fit the data better, until the rms of the difference between model and real data is smaller than a given amount (e.g. 0.2 ADU/gain/s ). Note also that the effects generated by the presence of glitches in the nearby pixels are considered by the fitting algorithm. The code recognizes sources above a given flux level, which decreases as the reduction improves the fit. In the pixels around relatively strong sources ($`>`$ 1.3 ADU/gain/s) we force the fit to find sources, leaving the fit level free. Once a satisfying fit is obtained for all the pixels over the whole pixel history, the flat field is computed from the stabilization level of the background. In the raster structure we set the flat fielded smoothed differences of: a) data readouts minus local fitted background (‘unreconstructed data’); b) fully recovered intensities minus stabilization level (‘reconstructed data’). Glitches and bad data are masked and this mask is stored in the raster structure to be used later in the map creation.

Then the reduced images per raster pointing are computed and corrected for flat-field distortions. Finally, the images are combined together to create the final raster maps (one for each raster position), where we then look for source detection.

A final reduction stage is performed after source extraction, simulating the data we would have from these detections and correcting the pixel fit, forcing the algorithm to recognize the source whenever this had not happened correctly (i.e. the source had been recognized only in one of the two overlapping single pointing images and not in the other).

We will now enter into more detail on the map creation and source detection processes.

### 4.2 Map creation

Once the images, corresponding to each raster position, have been created by averaging together all the time-scans relative to that pointing, they are converted from ADU/gain/s to mJy/pixel using the ISOCAM User’s Manual calibration factor (e.g. dividing by 1.96) and flat fielded. Also the number of un-masked times-scans (raster.npix) are scaled with flat-field coefficients, on the simple assumption of a constant noise in the data prior to flat-fielding.

After that, they are projected onto a sky map (raster image) using a simple TAN projection. The algorithm used is part of the CIA package (projette.pro). It computes the values of pixels on the sky map by averaging the pixel values in the single images, with a weight equal to the number of useful time-scans, to give the pointing image (raster.npix). For each raster map, two corresponding maps are also constructed, with the same size and sky orientation. One is the map containing, for each pixel, the number of frames coadded together (excluding the masked pixels) to obtain that pixel value in the raster (“NPIX” map). The other is the “RMS” map, where the rms of each pixel has been computed by scaling the measured mean rms of the central part of the map according to the inverse square root of “NPIX”.

In projecting each raster pointing onto the sky, the algorithm takes into account the field distortions of ISOCAM, as measured by Aussel et al. (1999). These distortions are a chromatic effect which cause the pixel size to be non-uniform on the sky (border pixels are larger in area than central pixels). This effect must be considered also when computing the flat-field. For a more detailed description of the ISOCAM field distortions see Okumura (2000).

In figure 3 the grey-scale image of the central raster map (S1$`\mathrm{\_}`$5) is shown. Note that this image, $`S1\mathrm{\_}5`$, is the combination of three single observations, this field having been observed with a redundancy of 3 compared to the other rasters.

The projection algorithm strongly affects the appearance of point sources on the map, having the general effect of smoothing the PSF over several pixels. As we will show in section 5.2, the peak flux values of sources with the same total flux can differ significantly, by a factor of up to 2, depending on the source position within the raster pixel.

### 4.3 Source detection

Before performing the proper source detection (on the final maps) to produce our definitive source catalogue, we have identified candidate sources inside the pixel histories. This was also very useful to check the confidence level of our fit to the data. Since during data reduction we have created a model for the background, we identified sources in the history of pixels from their flux excess above the background over the single time-scans. We inspected by eye every excess greater than 0.7 ADU/gain/s, correcting the few cases corresponding not to real sources, but to algorithm failures (by re-setting the parameters and starting again with model fitting until convergence is achieved). We have found that our method is very conservative, in the sense that cases where a spurious source is created by the algorithm constitute a very small fraction of the total number of correct detections, while the fraction of good sources missed by the fitting algorithm is rather significant (since sources are normally seen on several pixels, a lost detection on the pixel history does not necessarily mean the source is lost on the map).

Only faint sources remaining undetected on the map (because of these failures) contribute to real incompleteness, while for most of the brighter sources the failures result in a decrease of their total flux. In the final stage of our reduction (see below), we re-project the sources detected on the raster map into the pixel time-line, allowing a better fit of the data for all the sources that will appear above the interactive check threshold. This job significantly reduces the flux defect for the detected sources (but is of course unable to recover sources that fell below the detection threshold, i.e. to correct for incompleteness).

Concrete determination of the fraction of detected sources versus real sources, which leads to the estimate of completeness and reliability of our source catalogue, have been performed using simulations. This will be discussed in section 5.

After this preliminary check, we have searched for detections in the single calibrated images by selecting and visually inspecting all pixels with flux $`>0.4`$ mJy / pixel. In this case too, we have performed again the reduction for those pixels where the algorithm had failed to fit the data, thus producing a false detection.

These two checks on candidate sources, which required corrections and further cleaning for some pixels, provided very reliable (not complete!) source lists and images. After that, we could be confident that almost no spurious sources were present in our data-set. Therefore, we could proceed to the proper source detection. We must now point out that all the checks performed on the single pointing images and pixel histories do not guarantee that all these (and only these) sources will then be detected on the final raster map. In fact, as already discussed, the raster pixels are produced by combining together different single pointing images, where the same source could be located in different positions inside the pixels, thus being for example above the 0.4 mJy/pixel check threshold in one image and below this threshold in the other. For this reason, the list of sources obtained in the single pointing images for our preliminary checks are not always coincident with the final list derived on the raster map, where each source is determined by different effects. Moreover, in the raster map creation there are also distortion effects. Thus, the completeness and reliability of our final source lists can be tested only through simulations.

The source detection is done on the final raster maps, but is done using the signal-to-noise ratio. First, we have selected all pixels above a low flux threshold (0.1 mJy / pixel) using the IDL Astronomy Users Library (accessible via the World Wide Web home page http://idlastro.gsfc.nasa.gov/homepage.html) task called find. This algorithm finds positive brightness perturbations (i.e. stars) in an image, returning centroids and shape parameters (roundness and sharpness). Then, we have extracted from the selected list only those objects having a signal-to-noise ratio $``$ 5.

As discussed earlier, the LARI method is able to ‘reconstruct’ the source flux from transient effects. However, as we will clearly show with simulations (see section 5.2), the algorithm does not ‘reconstruct’ faint fluxes, corresponding to sources that it is not able to recognize. Therefore, faint sources have ‘reconstructed’ fluxes similar to their ‘unreconstructed’ ones, while for bright and correctly ‘reconstructed’ sources the ‘unreconstructed’ flux is, on average, about 1.7 times smaller than the corresponding ‘reconstructed’ flux. For this reason, to have a homogeneous flux determination for all our sources (both bright and faint), we have chosen to run the detection algorithm on the ‘unreconstructed’ maps. The correction for transients effects has then been performed individually in a second step, by using for each source its “effective PSF” when deriving its total flux (“auto-simulation” procedure, see section 5.1).

In order to achieve better position determination, we have run the detection algorithm on higher resolution maps, obtained by re-binning the original raster maps with a pixel size of 2 arcsec. The positions and fluxes given in output by find are determined by a convolution with a Gaussian PSF of given full width at half maximum (FWHM). We have chosen a FWHM = 9.8 arcsec, which is slightly smaller than the average FWHM of the ELAIS LW3 “effective PSF”.
The fluxes given by find are peak fluxes (mJy / pixel), which, coming from a Gaussian convolution, do not always correspond to the real source peaks. Therefore, for source peak flux (fs) we have given the maximum pixel value found within a box of 4 $`\times `$ 4 arcsec around the maximum given by find. To obtain the total fluxes (in mJy) we have used (and compared) two different methods: direct aperture photometry on the maps and “auto-simulations”, as discussed in detail in section 2.

## 5 Simulations

Because the raster maps on which we have performed the source detection are derived from the combination of several single pointing images, the only way to evaluate the effects produced on sources in the combined maps is through simulations.

With simulations, we can study the completeness and reliability of our detections at different flux levels and estimate the source positional accuracy and the internal calibration of the source photometry.

We added randomly distributed point sources to each of the three overlapping central raster maps (S1$`\mathrm{\_}`$5, S1$`\mathrm{\_}`$5$`\mathrm{\_}`$B and S1$`\mathrm{\_}`$5$`\mathrm{\_}`$C) at five different total fluxes (200 at 0.7, 150 at 1, 200 at 1.4, 150 at 2 and 150 at 3 mJy). It must be pointed out that our simulations have not been performed over the entire flux range covered by our survey, but only at the faint end. The reason is that we choose to sample with a high statistical significance the flux range more affected by incompleteness effects due either to mapping undersampling or data reduction method failures.

In a similar way we get simulations on the combined mosaiced map for which we have 50 sources at the same five flux levels.

To perform the simulations, we have created a high resolution map ($`1^{\mathrm{`}\mathrm{`}}`$) with simulated source, taking the ISOCAM PSF into account. The PSF has been successfully modeled by Okumura (1998) for stars. The PSF varies with the wavelegth, and since the ISOCAM filters are large, the shape of the PSF depends on the assumed spectrum of the point source. For our purpose, we have recomputed a model, following the prescriptions of Okumura (1998) but using a spectrum of the form $`f_\nu =constant`$, that is a closer match to the expected galaxy spectrum than the Rayleigh-Jeans form used for stellar spectra. The resulting PSF is larger than the one computed on stars.

Inverting the flat field and converting in ADU/gain/s we obtained the flux excess corresponding to the simulated sources. The value of this flux excess was then added to the real pixel histories (containing glitches and noise) by using Lari model, to obtain the simulated data cube. This simulated flux, included in the “liscio structure” (see section 4.1), have been reduced exactly in the same way as we did for all the original data structures, doing the same checks and repairs. In the produced maps we extracted the simulated sources with the same procedures used for the real rasters, measuring the resulting positions and peak fluxes. The peak fluxes measured after the data reduction will be referred as fs and fsr (as well as for real sources) respectively for ‘unreconstructed’ and ‘reconstructed’ maps. The corresponding theoretical peak fluxes associated to the excess flux maps, not reduced and containing neither glitches or noise, will be named f0 and f0r. These two sets of parameters then allow to evaluate separately the effects produced by the ELAIS observational strategy and the ISOCAM instrument only (f0 and f0r; see section 5.1) from the effects produced by the LARI reduction method applied to ELAIS data (fs and fsr; see section 5.2).

### 5.1 Theoretical transient behaviour of the detector

By simulations of the theoretical transient behaviour of the detector, we mean simulations of the effects due to the finite spatial resolution (6 ′′) and to the finite integration time (10$`\times `$2 sec) of our observations. We need to consider the spatial resolution of our observations, since the PSF is comparable in size to a pixel, causing the observed point source to depend strongly on its position on scales smaller than the pixel size.

Regarding the finite integration time, we need to simulate the fact that the CAM detector does not reach immediately the level corresponding to a given input flux, but needs a certain time to stabilize (see section 3). This stabilization effect, which means that only a fraction of the incident flux is detected, is not constant, but depends on the length of the time spent by the detector on target (not always the same) and on the amount of masking in a pixel (due to ‘glitches’ and uncertainties on the time spent on positions).

With the positions measured on simulated maps we can simulate how sources would appear in the ELAIS rasters if no noise and ‘glitches’ were present. To do this, we created two maps for each raster. The first is obtained by projecting back the simulated sources, injected on measured positions, onto the single pointing images and then computing the resulting raster map. The peak fluxes measured on it are the ‘reconstructed’ peak fluxes: $`f0r`$. For the second map we went back to the pixel history, predicting the behaviour due to the finite integration time transient (only theoretical, without any noise) and then producing a raster map. The peak fluxes measured on it are the ‘unreconstructed’ peak fluxes: $`f0`$.

Figure 4 shows the distribution of the peak fluxes (‘unreconstructed’ peak flux: $`f0`$; ‘reconstructed’ peak flux: $`f0r`$) normalized to 1 mJy (divided by the corresponding input flux, i.e. the effective PSFs), for simulated sources. The effect of mapping is the main responsible for the large spread of values in the figure, because the distribution of the $`f0r/f0`$ ratios has a rms of only 0.09 (0.05 for the repeated field S1$`\mathrm{\_}5`$) around a mean value of 1.66. The position uncertainty is only a minor contributor to the observed spread, since it causes only $``$14% of the $`f0`$ and $`f0r`$ distribution dispersion.

The ratios of sources detected over the total number of injected sources, due to PSF under-sampling and finite integration time are reported in the second column of Table 2 for each input flux.

The simulation of the mapping and mapping $`+`$ transient effects provides the estimate of the individual PSF for each source and gives a technique to be used to derive total fluxes for all the sources: for each source we will have two individual PSFs, one from $`f0r`$ and one from $`f0`$.

As we will show in section 2, there is a tight correlation between the peak flux of a source and its theoretical peak flux due to mapping and transients only (“effective PSF”). These “effective PSFs” can be used for aperture flux determination (using small radii) and we will show that the total fluxes derived from the observed peak fluxes correspond very well to the aperture photometry ones.

### 5.2 Real transient behaviour of the detector

Since the data reduction method can cause incompleteness in the final source list, we must take into account the effects produced by the fit when deriving the corrections to be applied to our catalogue. In fact, our data reduction method is based on a fitting algorithm and, depending on how well this is able to model the background, the ‘glitches’ and the sources, our catalogue will be more or less complete. For this reason, with simulations we have also tested the effects of Lari model on the final data products.

As stated above, to test the data reduction method, we have followed the same procedure for the simulated data that we used for real data. These simulated data cubes contain both real sources and simulated ones. They have also the same rms noise, all the ‘glitches’, ‘faders’, ‘dippers’ and background transients as the original data. The confusion will be slightly increased, but this effect is not critical for ELAIS data.

By comparing the output fluxes (per pixel) obtained for the simulated sources affected only by the mapping effects ($`f0`$) with the output fluxes of the reduced simulated sources ($`fs`$), we find a correlation (although not a 1 to 1 correlation, since the reduced fluxes are always slightly lower than the unreduced ones, see figure 5(top)). A similar correlation is observed for the corresponding ‘reconstructed’ peak fluxes ($`f0r`$ and $`fsr`$, reconstructed from the transients), although for faint sources our algorithm is not able to reconstruct correctly the fluxes (see figure 5($`bottom`$)).

The dispersion of values in figure 5(top) is caused by two kind of errors:

* an error proportional to the flux caused by the reduction method limits or by the mapping and finite integration time effects. In either cases, this error affects the peak fluxes.
* an additive error caused by the presence of noise and confusion. This error is more effective at low fluxes than at high ones.

At low fluxes, the combination of these errors may cause the total loss of a source (i.e. incompleteness).

Figure 6 shows the distribution of the ratio between $`fs`$ and $`f0`$, for the 424 simulated sources detected above $`5\sigma `$, compared to the same ratio for the 147 detections at 3 mJy only. As already noted in figure 5($`top`$), the ratio is not centered on the value of 1, thus causing a general underestimation of the total fluxes (derived from $`fs`$). There are two reasons for this effect: one is the fact that $`f0`$ is computed on the measured positions of $`fs`$ (about 11% higher than on real positions); the other is the under-evaluation of the wings of faint sources by the LARI method. The former effect is caused by the fact that projection effects cause $`f0`$ to be enhanced at favorable positions (i.e. center of a pixel on the single images), affecting also the centroid position even in the absence of noise. This overestimates the peaks of simulated sources which do not fall on the center of a pixel, thus leading to a bias in the ratio between the real $`f0`$ peak flux and the measured $`f0`$ simulated peak flux. The values of $`f0`$ computed at the positions measured for $`fs`$ are on average $``$11% higher than the ones computed at real positions.
The total distribution is larger and with a longer tail than the 3 mJy one. This is caused by the noise contribution to the errors, which is more significant at low than at high fluxes, as shown in figure 7. In this plot, the predicted distributions of the ratio $`fs/f0`$ in presence of a noise of 26 $`\mu `$Jy are shown for three different values of $`f0`$, corresponding to different mean values of $`f0`$ for different input fluxes (i.e. $`<f0>=605\mu `$Jy for 3 mJy input flux). As we can notice, the presence of noise broadens the flux distributions and this effect becomes stronger towards fainter fluxes. The predicted and unbiased distribution of $`fs/f0`$ for the brighter sources (3 mJy) peaks at 0.78 $`\pm `$ 0.03, value subsequently assumed to correct our measured fluxes (see section 2).

If we assume that the 3 mJy distribution reflects all the multiplicative error components due to the data reduction and we correct for the noise effect, we can predict the distribution of the detection rate for all the total and peak fluxes considered and derived in the simulation. Figure 8 shows the ratio of found-to-predicted detections as a function of $`f0`$, with the predictions obtained by considering two sources of error: (1) the multiplicative errors only (due to the reduction method and mapping, see section 2; solid line); (2) the multiplicative errors plus the additive noise contribution, assuming a typical noise level on maps of 26 $`\mu `$Jy (dashed line). It is visible a decrease of detection rate/predicted rate below $``$300 $`\mu `$Jy/pixel, which corresponds roughly to 1.5 mJy in total input flux. This deficiency is the nominal incompleteness of our data reduction method.

The detection rates at different fluxes derived with our simulations are reported in Table 2, where the predictions if only mapping smearing would be present, the predictions considering the data reduction but without taking into account the incompleteness of our method (see figure 8), the predictions considering also the incompleteness curve and the found values are given as detection fractions (i.e. detected sources / input sources) for the five different input fluxes. In figure 9 the detection rate curves relative to the values given in Table 2 are shown as function of the input flux.

As shown both in Table 2 and in figure 9, almost all the injected 3 mJy sources are detected at $``$ 5$`\sigma `$. Thus, the detection rate is 98.7% above 3 mJy, and it remains above 85.9% at 2 mJy. However, the detection rate drops quickly at fainter fluxes, in fact it reaches 52.3% at 1.4 mJy and 24.8% and 4.4% at 1 and 0.7 mJy respectively. Both sampling and reduction method are responsible for the large source undetectability at faint fluxes, although the contribution due to LARI method seems to become more important around 1.4 - 1 mJy, then it keeps almost constant, while the PSF sampling effect significantly decreases the detection rate for fluxes fainter than 1 mJy.

It must be pointed out that these incompleteness factors cannot be directly translated to the real catalogue, which has not a monochromatic flux distribution as the simulations. These factors can however be used to obtain the completeness of the catalogue and the source counts corrections, assuming a model for the $`logNlogS`$ (see Gruppioni, Lari, Pozzi et al. 2001). When applying simulations to real sky maps we must also remember that there are other sources of error not included in the simulations, as flat field corrections and distortion tables. The latter causes a higher smearing of the images and a larger uncertainty on centroid positions.

## 6 Flux determination

The simulation procedure described in section 5.1 has been used also to estimate the effective PSF on a source, real or simulated, and its total flux, given its position and peak flux. This procedure, performed on sources to determine their total flux from their peak flux and position is called “auto-simulation”.

The linear relation existing between $`fs`$ and $`f0`$ (see figure 5($`top`$)) allows the definition of a flux estimate for both simulated and real sources:

$$f_{total}=fs\times (flux/f0).$$
(1)
While for simulations $`flux`$ is the injected total flux, for real sources we need to adopt a rough estimate for $`flux`$ to derive $`f_{total}`$. $`Flux`$ is the total injected flux used to compute $`f0`$ (for simulations is 0.7, 1, 1.4, 2 and 3 mJy). Since transient corrections depend (slowly) on $`flux`$, for real sources we started with a rough estimate for $`flux`$ and then we iterated equation 1 to obtain a good estimate of $`f_{total}`$ also for strong sources.
The starting rough estimate of $`flux`$ is obtained by:

$$flux=fs/<fs/flux>_{sim}$$
(2)
where $`<fs/flux>_{sim}`$ = 0.132 was the average value taken from simulations.
Given this input total flux, we can derive $`f0`$ for real sources exactly as we did for simulated sources (see section 5.1). Then, by using formula 1, corrected for the systematic bias of the $`fs/f0`$ distribution (i.e. divided by 0.78; see section 5.2), we obtain the value of the total flux, $`f_{total}`$, for our sources. Given the relation 1, figure 5 could also be seen as the representation of the $`f_{total}/flux`$ (i.e. measured flux / true flux) distribution.

In figure 10 the total fluxes obtained with this procedure are plotted as function of the reduced peak fluxes obtained for all the simulated sources. As in figure 5, the open circles represent the sources detected above 5 $`\sigma `$, while the dots are the sources below the 5 $`\sigma `$ threshold.

The total fluxes obtained with the auto-simulations for the simulated sources are then compared with the total fluxes obtained with aperture photometry of the same simulated sources.

Concerning the aperture photometry fluxes, we found that the better determination was achieved with an aperture radius of 8\`\`, after correcting for the missing flux outside the aperture (40% in the PSF wings at distance $`>8^{\mathrm{`}\mathrm{`}}`$). We have found a good agreement between the two flux determinations (see figure 11 for a comparison between the two results). As total flux estimates for our real data sources, we have then decided to adopt the fluxes obtained from the auto-simulations.

In addition to the systematic bias affecting the $`f0`$ measured values, there is also a flux-dependent bias at low signal to noise levels, which derives from the fact that only sources with a high $`fs/f0`$ value and with positive noise fluctuations can be detected. However, only constant bias corrections are applied to our catalogue data.

### 6.1 Flux errors

As already mentioned in section 5.2, there are two main sources of uncertainty on our source fluxes: a multiplicative error due to mapping and data reduction method and an additive error due to the presence of noise in the map (neglecting the uncertainties due to flat-fielding corrections and field distortions, the latter always depressing the peak fluxes). As mentioned in section 5.1, the spread in the $`f0`$ distribution caused by position errors is about 14%. This spread is not only an important cause of the total flux bias, but it contributes significantly to the width of the $`fs/f0`$ distribution ($``$0.18) at high fluxes. The extra contribution from data reduction is about 11%. Since simulations show that this spread is rather insensitive of fluxes, we assumed that the multiplicative error is constant.

Because our total fluxes are obtained from the ratio between peak fluxes and auto-simulated peak fluxes, the combination of the two errors leads to a flux-dependent distribution like the one shown in figure 8.

Being the width of the $`fs/f0`$ distribution equal to 0.18 at high signal to noise levels, the distribution convolved with noise will have a width:

$$w^2=(0.18)^2+(rms/f0)^2$$
(3)
We used this relation to obtain the relative flux errors for sources.

## 7 Test of the photometry

The photometric accuracy of our reduction of the S\_1 area can be tested using the stars of the field. Aussel & Alexander (in prep.; see also Alexander & Aussel, 2000) have performed a detailed study of the mid-infrared emission of stars, from large sample of sources drawn from the IRAS Faint Source Calalog with excellent counterparts in the Tycho-2 catalog (Hog et al., 2000). They show that the B-V color of stars is extremely well correlated with the B- color, where  is a magnitude scale constructed from the IRAS flux, following the prescriptions of Omont et al. (1999). This relation allows to predict accurately the 12 $`\mu `$m IRAS flux of a star, provided that its B-V is known, and is lower than 1.3. Stellar atmosphere models (Lejeune et al., 1998) show that for the spectral types hotter than K3 that the color criteria select, the ratio of the 15 $`\mu `$m flux to the 12 $`\mu `$m flux of stars is constant. We have therefore used the relation calibrated on IRAS data by Aussel & Alexander to predict the fluxes of stars in the field, and we compare them to the product of our analysis.

The area surveyed in S\_1 contains 170 stars in the Tycho-2 catalogue, 145 of which with $`BV<1.3`$. In our analysis we detect 63 of them, 48 with $`BV<1.3`$. We plot on figures 12 and 13 respectively the measured fluxes versus the predicted fluxes and the histogram of the ratio of the measured flux to the predicted flux. In figure 12 the dashed line shows the one-to-one relation, followed by our data over more than two order of magnitude in flux. In figure 13, the dotted line shows the ratio of measured over predicted IRAS 12 $`\mu `$m fluxes, for a sample of 3950 stars from the study of Aussel & Alexander. The distribution is a skewed log-normal, dominated by the error of the IRAS FSC photometric error of the order of 10% on average. It is skewed toward observed fluxes higher than predicted fluxes, because some stars present an excess of IR emission due to the presence of a disk or shells. Dashed line is the result of the covolution of the former distribution with the $`fs/f0`$ strong sources distribution to simulate the spread of values we would expect in our analysis, neglecting noise. The mean value of this distribution is 1.047 while the 48 stars in S\_1 with $`BV<1.3`$ have a mean value of .955 leading to a relative flux scale of $`1.096\pm .044`$.

The solid line is the ratio of the measured LW3 flux and predicted fluxes for the 48 stars detected in S\_1.

The shape of the distribution is the same as the dashed one, apart the small scale factor, with the same skewness and we are confident our fluxes are correct, over a large range of fluxes since these stars go from 0.85 mJy to 135 mJy in LW3.

## 8 Positional accuracy

The positional errors in RA and DEC for our sources can be considered as the combination of three different sources of uncertainty: the finite spatial sampling ($`\sigma _s`$), the reduction method ($`\sigma _r`$), and the uncertainties in the pointing accuracy ($`\sigma _p`$). The latter is due to errors in the ISOCAM lens position (the wheel jitter) and results in an offset of about 1.2 pixels from the optical axis, that translates to $``$7 arcsec with a pixel size of 6′′. Moreover, ISO absolute pointing accuracy is about 3 arcsec.

The effect of the finite spatial sampling ($`\sigma _s`$) has been evaluated from the “theoretical” simulation (see 5.1), considering the distribution of the differences between the positions of the injected sources (RA,DEC) and the positions of the (same) sources detected in the projected map (RA0, DEC0).
The sum of this effect plus that produced by the method of reduction ($`\sigma _{s+r}`$) has been evaluated from the “real” simulation (see 5.2), considering the distribution of the differences between the positions of the injected sources (RA,DEC) and the positions of the sources detected in the projected map after the reduction (RAS, DECS). The widths of these distributions are 0.63 (RA) and 0.91 arcsec (DEC) for sampling only and 1.17 (RA) and 1.27 arcsec (DEC) for sampling and reduction effects. In figures 14 we plot the distributions of the differences in RA and DEC between the injected and detected positions.

By using our simulations at different input fluxes, we have also checked the dependence of the positional errors on source signal-to-noise ratio, as shown in figure 15. While the positional accuracy due to sampling only is almost constant with signal-to-noise ($`\sigma _s`$0.9′′ for DEC and $`\sigma _s`$0.65′′ for RA), as expected being a pure geometrical factor, the positional accuracy after the reduction is strongly dependent on signal-to-noise, increasing by about 50% from $`S/N10`$ to $`S/N57`$ (i.e. $`\sigma `$1.0′′ for RA at $`S/N10`$; $`\sigma `$1.5′′ at $`5S/N7`$). These dependences can be approximated by exponential laws of the form :

$$\sigma _{s+r}(RA)=1.0+17.17\times e^{(0.57\times S/N)}$$
(4)
$$\sigma _{s+r}(DEC)=1.06+1.21\times e^{(0.16\times S/N)}$$
(5)
These laws, plotted in figure 15 as solid lines and found with a non-linear least squares fit, have been used to estimate the positional errors due to the mapping and reduction method as a function, for each source, of its signal-to-noise.

The errors introduced by uncertainties in the ISOCAM pointing can be estimated by performing optical identifications for the sources found in each raster and computing an offset with respect to the optical astrometric reference system. As optical reference list we have used the PMM USNO-A2.0 Catalogue (Monet 1998).

With between 28 and 36 ISOCAM/USNO coincidences per raster found, we derived the median offsets for each of the 11 frames using the following procedure. First, we have cross-correlated the ISOCAM and USNO lists using a maximum distance of 60 arcsec to obtain the best value for the maximum distance for reliable identification. This distance resulted to be 12 arcsec for all the rasters. Then, for each raster we have obtained the median offset values in RA and DEC ($`med_{ra1}`$ and $`med_{dec1}`$) using all the ISO-USNO sources with a maximum distance of 12 arcsec. We applied the offset values to the ISO positions and we have cross-correlated again the ISO and the USNO catalogues.

We have selected again all the sources within a maximum distance of 12 arcsec and we have used these sources to calculate new median offset values ($`med_{ra2}`$ and $`med_{dec2}`$). The total offsets for each raster have then been obtained as $`RA\mathrm{\_}offset=med_{ra1}+med_{ra2}`$, $`DEC\mathrm{\_}offset=med_{dec1}+med_{dec2}`$. The values of the offsets and their relative errors (computed as the standard errors on median: $`\sigma _{med}=1.2533\frac{\sigma }{\sqrt{N}}`$ (Akin & Colton 1970), where $`\sigma `$ and $`N`$ are respectively the standard deviation and the number of sources considered in each raster) have been reported in Table 3. Each source position has then been corrected for the offset found for the corresponding raster. The error ($`\sigma _p`$) introduced on source positions by the presence of the systematic offset is given by the error on the offset determination (see column 4 and 6 in Table 3). This error has been added to the positional uncertainty due to mapping and reduction method ($`\sigma _{s+r}`$) to obtain the total position error for each source:

$$\sigma _{RA}^2=\sigma _{s+r}^2(RA)+\sigma _p^2(RA)$$
(6)
$$\sigma _{DEC}^2=\sigma _{s+r}^2(DEC)+\sigma _p^2(DEC)$$
(7)
In figure 16 is shown an example of our ISOCAM 15 $`\mu `$m contour levels superimposed to DSS optical images, after correcting for the systematic offsets computed above. The plot can give an idea of the astrometric accuracy of our catalogue and images. In fact, as clearly visible from the figure, the positions of our sources after offsetting appear as accurate as we estimated from simulations (see above), thus allowing reliable identifications within a few arcseconds or less.

## 9 Repeated central region, S1$`\mathrm{\_}`$5

The central field of the S1 area, as mentioned in section 4.2, has been observed three times in order to reach a deeper flux limit with respect to the rest of the area and to allow reliability checks on sources.

The reduction of the three observations was carried out in the standard way (see section 4.1) until the stage of map creation. After the creation of the three single raster maps, some particular routines have been applied for combining them and for performing simulations in the combined map.
To obtain a unique combined map from the three single observations, first we have projected all of them on the sky with the same orientation (north-south). The three single rasters have then been corrected for the relative astrometric offsets (see section 8) and then coadded. For the coaddition, each raster map has been weighted, with a weight proportional pixel-by-pixel to its relative “NPIX” map. The combined “NPIX” map was just the pixel-by-pixel sum of each “NPIX” map. The “RMS” distribution over the mosaic map has been obtained with the standard procedure (see 4.1). The average rms value in the central part of the combined S1$`\mathrm{\_}`$5 map is about 0.016 mJy.

Once obtained the combined map, source extraction have been performed in the same way as for the single observation raster maps (see section 4.3). In S1$`\mathrm{\_}`$5 we have detected 93 sources.

To derive the total fluxes from the detected peak fluxes, we have performed the “auto-simulations” (see section 2) for the 93 sources, by injecting point sources into each of the three single fields and then combining the resulting images with the same weight used to coadd the real maps. The total-to-peak ratio found for the simulated sources have then been applied to the peak flux obtained for the real sources to get their total flux, exactly in the same way as we did for the single rasters in the rest of the S1 area.

Total fluxes in the combined map range between 0.57 and $``$100 mJy.

### 9.1 Simulations in S1$`\mathrm{\_}`$5

To perform simulations of the repeated raster we have used and appropriately combined the simulations performed separately on the three single rasters, S1$`\mathrm{\_}`$5, S1$`\mathrm{\_}`$5$`\mathrm{\_}`$B, S1$`\mathrm{\_}`$5$`\mathrm{\_}`$C (section 5). The positions of the sources injected in each of the three individual fields have been chosen in order to simulate the effects caused by the application, in the coaddition phase, of a relative astrometric offset among the rasters. 50 sources for each of the above 5 total fluxes (0.7, 1, 1.4, 2 and 3 mJy) were injected. The “auto-simulations” were performed on the positions found on the combined map.

In figure 17 the $`fs`$ and $`f0`$ peak fluxes of the detected sources are plotted superimposed to the ones found for the simulated sources in the main S1 area. Apart from the deeper detection level, the general trend is the same with a somewhat smaller dispersion. For the 3 mJy input sources the $`fs/f0`$ observed distribution peaks at 0.82 $`\pm `$ 0.03 and its width is 0.11, while the corresponding values for single fields simulation are 0.78 $`\pm `$ 0.03 and 0.18.

Following the same procedure as before we can predict completeness and detection rates also for sources in the combined map.

The results of the simulations have been reported in table 4 and shown in figure 18.

The combination of the three maps does not only reduce the errors in flux determination and increase the fraction of detected sources at faint fluxes, but provides more precise positions in the sky.

Figure 19 shows the distributions of the differences in RA and DEC between the injected and detected positions, while figure 20 shows the dependence of position errors on the signal-to-noise ratio. The widths found for the distributions in RA and DEC for sampling and reduction effect, are 0.92 and 0.85 respectively, smaller than those found for the main survey, 1.17 and 1.27. Considering the dependence of the position errors on the signal-to-noise ratio, the laws found are of the form:

$$\sigma _{s+r}(RA)=0.74+1.59\times e^{(0.2\times S/N)}$$
(8)
$$\sigma _{s+r}(DEC)=0.54+0.86\times e^{(0.07\times S/N)}$$
(9)
These law are less steep than those found in the all survey and, as we can see from the figure 20, the values of the positional errors near the limit of the survey (signal-to-noise = 5) are, for both coordinates, less than 1.5′′.

The combination of the three maps, changing the repetition factor from 2 to 6 for each single pointing image, not only reduces the errors, but also the effects due to mapping.

## 10 The Source Catalogue in S1

The final catalogue obtained with our method contains a total of 462 sources detected at 15 $`\mu `$m (LW3) in the ELAIS region S1. All the sources detected over the whole $`2^{}\times 2^{}`$ area have signal-to-noise ratio greater than 5. The catalogue reports the source name, the offset-corrected position (right ascension and declination at equinox J2000), the positional accuracy on the final images, the source peak flux (in mJy/pixel), the detection level (signal-to-noise ratio), the total flux and its error (in mJy), the raster name and eventually a note, indicating whether the source is identificated with a star, its flux have been obtained through aperture photometry, etc. In case of extended of very bright sources, the total flux reported in the table is computed by aperture photometry instead of by “auto-simulations”, the latter method providing a correct measurement mainly for unresolved sources. Note that a few sources (belonging to the border part of a raster, overlapping with an adjacent raster) might appear in two different rasters. In this case, the repeated sources have been reported twice in the catalogue and the corresponding additional raster number is quoted in the notes.

As described in section 4.3, before extracting sources on the final maps, we have identified candidate sources inside the pixel histories (flux excesses above the background over the single time-scans greater than 0.7 ADU/gain/s) and on the single calibrated images (all pixels with flux $`>0.4`$ mJy / pixel), providing 100% reliable lists of sources above these two thresholds. Although there is not a perfect correspondence between these two flux thresholds and total fluxes in the final raster maps (due to flat-fielding, distortions, etc.), we have found that by splitting the S1 catalogue in two (above and below 1 mJy) we can provide a highly reliable catalogue, where all the sources have been checked before extraction and a less reliable but deeper catalogue, where most sources could not be checked with the above criteria. The majority of sources fainter than 1 mJy, in fact, might have flux excesses in the single pixel histories below 0.7 ADU/gain/s and peak fluxes on the single images fainter than 0.4 mJy / pixel, the limits chosen for visual inspection, below which is almost impossible to distinguish a flux excess on the pixel history from local background fluctuations. This does not apply to S1$`\mathrm{\_}5`$, because sources have been extracted on the combination of three single observations, that have been separately checked before the coaddition.

In Tables 5 and 6, the first page (corresponding to the first raster) of the catalogues in S1, respectively above and below 1 mJy, are shown as examples. The full ELAIS S1 + S1$`\mathrm{\_}5`$ catalogues at 15 $`\mu `$m obtained with Lari method will be available from http://boas5.bo.astro.it/$``$elais/catalogues/.

## 11 Conclusions

A new data reduction technique (the Lari method) has been successfully applied to the 15 $`\mu `$m ISOCAM observations of one of the four main ELAIS fields (S1). This technique, based on the existence of two different time-scales in ISOCAM transients, was particularly efficient in overcoming the main problems affecting the ISOCAM LW data and in detecting faint sources. Its application to the southern ELAIS field has produced a catalogue of 462 sources, detected above the 5 $`\sigma `$ threshold over an area of about 4 square degrees. The integrated fluxes of these sources cover the range 0.5 - 100 mJy, filling the existing gap between the Deep ISOCAM Surveys and the Faint IRAS Survey. The completeness and photometry accuracy of our catalogue have been tested through accurate simulations performed at different flux levels. The results of these simulations showed that our catalogue is highly reliable and $`>98.5`$% complete at 3 mJy. The completeness, due either to the mapping effects or to the data reduction method, then decreases at fainter fluxes. The positional accuracy, estimated with simulations, resulted to be about 1 arcsec in both right ascension and declination for signal-to-noise ratios $`>7`$, while it increases to $`2`$ and 1.6 at signal-to-noise ratios close to the survey threshold (5), respectively for right ascension and declination. The photometric accuracy of our data reduction has also been tested using the stars of the field, comparing the measured fluxes with the ones predicted by the relation calibrated on IRAS data by Aussel & Alexander (2000). Our fluxes resulted well consistent with the predicted ones over a large range of fluxes, since these stars go from 0.85 mJy to 135 mJy in LW3.

In a forthcoming paper (Gruppioni, Lari, Pozzi et al. 2001) we will present the source counts obtained from this survey in the crucial uncovered flux range 0.45 – 100 mJy, dividing the Deep/Ultra-Deep ISOCAM Surveys from the fainter IRAS Surveys.

## Acknowledgments

This work was supported by the EC TMR Network programme FMRX–CT96–0068. CG acknowledges partial support by the Italian Space Agency (ASI) under the contract ARS–98–119 and by the Italian Ministry for University and Research (MURST) under grants COFIN98 and COFIN99. We would like to thank Gianni Zamorani for helpful suggestions and for a careful reading of the manuscript and David Elbaz for constructive refereeing, which improved the quality of the paper.

## Appendix A Lari Model Description

The process is governed by two differential equations, one for each charge reservoir, of the form:

$$\frac{dQ}{dt}=eIaQ^2$$
(10)
where $`I`$ is the incident flux of photons, $`e`$ is the efficiency of the process feeding the component and $`a`$ is a time constant which depends on the detector pixel-size. Note that $`e`$ and $`a`$ assume different values for the two components: $`a(`$breve) $`>a(`$lunga) and $`e(`$breve) $`>e(`$lunga).

We have not attempted to model $`I`$ and $`e`$ for glitches. In principle, glitches could be described by the physics of ionising particles. However their effect strongly depends on the nature and energy of the incident cosmic particle. For example, high energy incident cosmic rays could produce saturation on the detector, thus causing the parameter $`e`$ to depend on the values of $`I`$ and $`Q`$ . For simplicity, we have neglected this effect in our model, considering both $`e`$ and $`a`$ as constants.

Our model is completely conservative (no decay of the accumulated charges is considered, except toward the contacts) and homogeneous (the charge reservoir involves all the detector parts that do not contribute to polarising the contacts). In fact, at stabilisation we have $`aQ^2=eI`$ and $`S=I`$ ($`S`$ being the signal), while generally $`I=S+\mathrm{\Delta }Q/\mathrm{\Delta }t`$. The quantity $`aQ^2`$ in the accumulated charge equation is exactly the same amount of charge which that component feeds the contacts with.

Considering charges as fluxes in ADUs ($`q=Q/\mathrm{\Delta }t`$), we have:

$$\mathrm{\Delta }t\frac{dq_1}{dt}=e_1Ia_1\mathrm{\Delta }t^2q_1^2$$
(11)
for the breve component. Integrating over an observation integration time $`\mathrm{\Delta }t`$ (with $`I=`$ constant):

$$\mathrm{\Delta }q_1=e_1Ia_1\mathrm{\Delta }t^2<q_1^2>$$
(12)
and

$$\mathrm{\Delta }t\frac{dq_2}{dt}=e_2Ia_2\mathrm{\Delta }t^2q_2^2$$
(13)
for the lunga component, which integrated over an integation time becomes

$$\mathrm{\Delta }q_2=e_2Ia_2\mathrm{\Delta }t^2<q_2^2>.$$
(14)
We then have

$$S(obs.signal)=e_0I+a_1\mathrm{\Delta }t^2<q_1^2>+a_2\mathrm{\Delta }t^2<q_2^2>$$
(15)
where $`e_0+e_1+e_2=1`$

The two differential equations are of Riccati type, which, for $`I=constant`$ have a general analytical solution:

$$q(t)=A\times \frac{((A+q(0))e^{bt}(Aq(0)))}{((A+q(0))e^{bt}+(Aq(0))}$$
(16)
where $`q`$ can be either $`q_1`$ or $`q_2`$. $`A`$ represents the asymptotic value for $`q`$ ($`A=\frac{\sqrt{eI/a}}{\mathrm{\Delta }t}`$), while $`b`$ ($`=2\sqrt{eIa}`$) is the inverse of the time–scale at stabilisation. Note that the time–scale in this model is inversely proportional to the square root of the (constant) incident flux (at stabilisation only!), while the observed flux $`S`$ tends to $`I`$ as $`t`$ tends to infinity.

In our code we make use also of an approximate equation with finite difference values to obtain valid solutions also for the general case of variable $`I`$:

$$(q(t+\mathrm{\Delta }t)q(t))=eI(t)a\mathrm{\Delta }t^2q(t+\mathrm{\Delta }t)q(t)$$
(17)
where $`q`$ is either $`q_1`$ or $`q_2`$, $`e`$ is either $`e_1`$ or $`e_2`$, $`a`$ is either $`e_1`$ or $`e_2`$, respectively for the breve and for the lunga components. $`q(t)`$ and $`q(t+\mathrm{\Delta }t)`$ are the charges that are accumulated respectively at the beginning and at the end of and integration. $`I(t)`$ is the average intensity over the whole integration.

The error we commit by use of this second-order approximation instead of the exact equation is less than 1%, so this can be considered a good approximation. The same kind of approximation is used in other parts of our code when calculating derivatives.

In equation 15 the observed flux is the flux subtracted by the dark current. In dark observations ‘glitch’ transients show as ‘faders’ and ‘dippers’, the latter having time-scales larger than the integration time, but not infinite (as the model requires). Thus, both lunga and breve charge productions are fed also by the thermal component of the dark current. This effect might be important when the photon flux is small. It is very difficult to estimate this extra source of transient signal and, since we could not find any documentation on ISOCAM thermal dark current measurements, we tried to estimate it from the data. We have thus associated the thermal dark current to the minimum amount of extra signal that is needed in order to keep the parameter $`e_2`$ below the value of the lunga fraction ($`0.1`$, implying that dippers cannot exceed one-tenth of the sky background counts).

By giving the name ‘offset’ to the thermal feeding, equations 12, 14 and 15 become:

$$\mathrm{\Delta }q_1=e_1(I+offset)a_1\mathrm{\Delta }t^2<q_1^2>$$
(18)
$$\mathrm{\Delta }q_2=e_1(I+offset)a_2\mathrm{\Delta }t^2<q_2^2>$$
(19)
$`S(obs.signal)+offset=`$ (20)
$`e_0(I+offset)+a_1\mathrm{\Delta }t^2<q_1^2>+a_2\mathrm{\Delta }t^2<q_2^2>`$
In practice, these two equations can be solved by successive iterations, provided that $`q_{1/2}(0)`$, $`e_{1/2}`$ and $`a_{1/2}`$ are known. Estimates for these parameters, characterising each pixel, can be obtained by minimising the $`\chi ^2`$ estimator, under the condition of having a constant incident flux at each raster position. In this model, all the past history of each pixel is contained in the intial charge, as long as there is no other source of electrons (eg from the surrounding pixels or from longer time relaxation processes inside the pixel).


format: text

text:
# Trapping oscillations, discrete particle effects and kinetic theory of collisionless plasma

## I Introduction

It is tempting to expect that kinetic equations and their numerical simulation provide a fair description of the time evolution of systems with long-range or ‘global’ interactions. A typical, fundamental example is offered by wave-particle interactions, which play a central role in plasmas. In this Letter we test this opinion explicitly.

Collisionless plasma dynamics is dominated by collective processes. Langmuir waves and their familiar Landau damping and growth  are a good example of these processes, with many applications, e.g. plasma heating in fusion devices and laser-plasma interactions. For simplicity we focus on the one-dimensional electrostatic case, traditionally described by the (kinetic) coupled set of Vlasov-Poisson equations . The current debate on the long-time evolution of this system shows that further insight in this fundamental process is still needed .

The driving process (induced by the binary Coulomb interaction between particles) is the interaction of the electrostatic waves in the plasma with the particles at nearly resonant velocities, which one analyses canonically by partitioning the plasma in bulk and tail particles. Langmuir modes are the collective oscillations of bulk particles, with slowly varying complex amplitudes in an envelope representation ; their interaction with tail particles is described by a self-consistent set of hamiltonian equations . These equations already provided an efficient basis  for investigating the cold beam plasma instability and exploring the nonlinear regime of the bump-on-tail instability . Analytically, they were used to give an intuitive and rigorous derivation of spontaneous emission and Landau damping of Langmuir waves . Besides, as it eliminates the rapid plasma oscillation scale $`\omega _\mathrm{p}^1`$, this self-consistent model offers a genuine tool to investigate long-time dynamics.

As we follow the motion of each particle, we can also address the influence of the finite number of particles in the long run. This question is discarded in the kinetic Vlasov-Poisson description, for which the finite-$`N`$ correction is the Balescu-Lenard equation  formally derived from the accumulation of weak binary collisions, with small change of particle momenta. It implies a diffusion of momenta, driving the plasma towards equilibrium. However, when wave-particle coupling is dominant, the Balescu-Lenard equation is not a straightforward approach to finite-$`N`$ effects on the wave evolution.

Here we investigate direct finite-$`N`$ effects on the self-consistent wave-particle dynamics. It is proved  that the kinetic limit $`N\mathrm{}`$ commutes with evolution over arbitrary times. As one might argue that finite $`N`$ be analogous to numerical discretisation in solving kinetic equations, we also integrate the kinetic system with a ‘noise-free’ semi-lagrangian solver . In this Letter we compare finite grid effects of the kinetic solver and granular aspects of the $`N`$-particles system, whose evolution is computed with a symplectic scheme .

We discuss the case of one wave interacting with the particles. Though a broad spectrum of unstable waves is generally excited when tail particles form a warm beam, the single-wave situation can be realized experimentally  and allows to leave aside the difficult problem of mode coupling mediated by resonant particles .

## II Self-consistent wave-particle model and kinetic model

Consider a one-dimensional electrostatic potential perturbation $`\mathrm{\Phi }(z,\tau )=[\varphi _k(\tau )\mathrm{exp}i(kz\omega _k\tau )+\mathrm{c}.\mathrm{c}.]`$ (where c.c. denotes complex conjugate), with complex envelope $`\varphi _k`$, in a plasma of length $`L`$ with periodic boundary conditions (and neutralizing background). Wavenumber $`k`$ and frequency $`\omega _k`$ satisfy a dispersion relation $`ϵ(k,\omega _k)=0`$. The density of $`N`$ (quasi-)resonant electrons is $`\sigma (z,\tau )=(nL/N)_{l=1}^N\delta (zz_l(\tau ))`$, where $`n`$ is the electron number density and $`z_l`$ is the position at time $`\tau `$ of electron labeled $`l`$ (with charge $`e`$ and mass $`m`$). Non-resonant electrons contribute only through the dielectric function $`ϵ`$, so that $`\varphi _k`$ and the $`z_l`$’s obey coupled equations

$`d\varphi _k/d\tau `$ $`=`$ $`{\displaystyle \frac{ine}{ϵ_0k^2N(ϵ/\omega _k)}}{\displaystyle \underset{l=1}{\overset{N}{}}}\mathrm{exp}[ikz_l+i\omega _k\tau ]`$ (1)
$`d^2z_l/d\tau ^2`$ $`=`$ $`(iek/m)\varphi _k\mathrm{exp}[ikz_li\omega _k\tau ]+\mathrm{c}.\mathrm{c}.`$ (2)
where $`ϵ_0`$ is the vacuum dielectric constant. With $`\alpha ^3=ne^2/[mϵ_0(ϵ/\omega _k)]`$ , $`t=\alpha \tau `$, $`\dot{}=d/dt`$, $`x_l=kz_l\omega _k\tau `$ and $`V=(ek^2\varphi _k)/(\alpha ^2m)`$, this system defines the self-consistent dynamics (with $`N+1`$ degrees of freedom)

$`\dot{V}`$ $`=`$ $`iN^1{\displaystyle \underset{l=1}{\overset{N}{}}}\mathrm{exp}(ix_l)`$ (3)
$`\ddot{x}_l`$ $`=`$ $`iV\mathrm{exp}(ix_l)iV^{}\mathrm{exp}(ix_l)`$ (4)
for the coupled evolution of electrons and wave in dimensionless form. This system derives from hamiltonian $`H(𝐱,𝐩,\zeta ,\zeta ^{})=_{l=1}^N(p_l^2/2N^{1/2}\zeta e^{ix_l}N^{1/2}\zeta ^{}e^{ix_l})`$, where a star means a complex conjugate and $`\zeta =N^{1/2}V`$. An efficient symplectic integration scheme is used to study this hamiltonian numerically .

The system (3)-(4) is invariant under two continuous groups of symmetries. Invariance under time translations implies the conservation of the energy $`𝖧=H`$. The phase $`\theta `$ of $`\zeta =|\zeta |e^{i\theta }`$ plays the role of a position for the wave, and system (3)-(4) is also invariant under translations $`\theta ^{}=\theta +a`$, $`x_l^{}=x_l+a`$. This translation invariance leads to the conservation of momentum $`𝖯=_lp_l+|\zeta |^2`$, where the contribution from the wave is analogous to the Poynting vector of electromagnetic waves (which is quadratic in the electromagnetic fields) . Conservation of these invariants constrains the evolution of our system, and we checked that the numerical integration preserves them.

In the kinetic limit $`N\mathrm{}`$, electrons are distributed with a density $`f(x,p,t)`$, and system (3)-(4) yields the Vlasov-wave system

$$\dot{V}=ie^{ix}f(x,p,t)𝑑x𝑑p$$
(5)
$$_tf+p_xf+(iVe^{ix}iV^{}e^{ix})_pf=0$$
(6)
For initial data approaching a smooth function $`f`$ as $`N\mathrm{}`$, the solutions of (3)-(4) converge to those of the Vlasov-wave system over any finite time interval . This kinetic model is integrated numerically by a semi-lagrangian solver, covering $`(x,p)`$ space with a rectangular mesh : the function $`f`$ (interpolated by cubic splines) is transported along the characteristic lines of the kinetic equation, i.e. along trajectories of the original particles .

Let us first study linear instabilities. One solution of (3)-(4) corresponds to vanishing field $`V_0=0`$, with particles evenly distributed on a finite set of beams with given velocities. Small perturbations of this solution have $`\delta V=\delta V_0e^{\gamma t}`$, with rate $`\gamma `$ solving

$$\gamma =\gamma _\mathrm{r}+i\gamma _\mathrm{i}=iN^1\underset{l=1}{\overset{N}{}}(\gamma +ip_l)^2.$$
(7)
For a monokinetic beam with velocity $`U`$, (7) reads $`\gamma (\gamma +iU)^2=i`$ ; the most unstable solution occurs for $`U=0`$ (with $`\gamma =(\sqrt{3}+i)/2`$). For a warm beam with smooth initial distribution $`f(p)`$ (normalized to $`f𝑑p=1`$), the continuous limit of (7) yields $`\gamma =i(\gamma +ip)^2f(p)𝑑p`$. For a sufficiently broad distribution ($`|f^{}(0)|1`$), we obtain $`|\gamma _\mathrm{r}|\gamma _\mathrm{r}=\gamma _\mathrm{r}\pi f^{}(\gamma _\mathrm{i})`$, where $`f^{}=df/dp`$, and $`\gamma _\mathrm{i}\pi \gamma _\mathrm{r}f^{\prime \prime }(0)`$ for $`|f^{\prime \prime }(0)|\pi ^1`$. Except for the trivial solution $`\gamma _\mathrm{r}=0`$, other solutions can only exist for a positive slope $`f^{}(0)`$. Then the perturbation is unstable as the evolution of $`\delta V`$ is controlled by the eigenvalue $`\gamma `$ with positive real part, i.e. with growth rate $`\gamma _\mathrm{r}\gamma _\mathrm{L}=\pi f^{}(0)>0`$. Negative slope leads to the linear Landau damping paradox : the observed decay rate $`\gamma _\mathrm{L}=\pi f^{}(0)<0`$ is not associated to genuine eigenvalues, but to phase mixing of eigenmodes , as a direct consequence of the hamiltonian nature of the dynamics.

Now, this linear analysis generally fails to give the large time behavior. This is obvious for the unstable case as non-linear effects are no longer negligible when the wave intensity grows so that the trapping frequency $`\omega _\mathrm{b}(t)=\sqrt{2|V(t)|}`$ becomes of the order of the linear rate $`\gamma _\mathrm{r}`$.

We used the monokinetic case as a testbed . Finite-$`N`$ simulations show that the unstable solution grows as predicted and saturates to a limit-cycle-like behavior where the trapping frequency $`\omega _\mathrm{b}(t)`$ oscillates between $`1.2\gamma _\mathrm{r}`$ and $`2\gamma _\mathrm{r}`$. In this regime, some of the initially monokinetic particles have been scattered rather uniformly over the chaotic domain, in and around the pulsating resonance, while others form a trapped bunch inside this resonance (away from the separatrix) . This dynamics is quite well described by effective hamiltonians with few degrees of freedom .

In this Letter, we discuss the large time behavior of the warm beam case, with $`f^{}(p_0)0`$ at the wave nominal velocity $`p_0=0`$. Fig. 1 displays three distribution functions (in dimensionless form) with similar velocity width : (i) a function (CD) giving the same decay rate for all phase velocities, (ii) a function (CG) giving a constant growth rate for all phase velocities , (iii) a truncated Lorentzian (TL) with positive slope $`f^{}(0)>0`$.

## III Damping case

For the damping case, the linear description introduces time secularities which ultimately may break linear theory down : the ultimate evolution is intrinsically nonlinear, not only if the initial field amplitude is large, as in O’Neil’s seminal picture , but also if one considers the evolution over time scales of the order of the trapping time (which is large for small initial wave amplitude). The question of the plasma wave long-time fate is thus far from trivial . Though some simulations  infer that nonlinear waves eventually approach a Bernstein-Greene-Kruskal steady state  instead of Landau vanishing field, the answer should rather strongly depend on initial conditions. Our $`N`$-particle, 1-wave system is the simplest model to test these ideas.

A thermodynamical analysis  predicts that, for a warm beam and small enough initial wave amplitude, $`\omega _\mathrm{b}N^{1/2}`$ at equilibrium in the limit $`N\mathrm{}`$. Fig. 2 shows the evolution of a small amplitude wave launched in the beam. The $`N`$-particle system (line N) and the kinetic system (line V) initially damp the wave exponentially as predicted by perturbation theory , for a time of the order of $`|\gamma _\mathrm{L}|^1`$.

After that phase-mixing time, trapping induces nonlinear evolution and both systems evolve differently. For the $`N`$-particle system, the wave grows to a thermal level that scales as $`N^{1/2}`$, corresponding to a balance between damping and spontaneous emission . For the kinetic system, initial Landau damping is followed by slowly damped trapping oscillations around a mean value which also decays to zero, at a rate decreasing for refined mesh size. Fig. 2 reveals that finite-$`N`$ and kinetic behaviors can considerably diverge as spontaneous emission is taken into account. The time $`\tau _N`$ after which the finite-$`N`$ effects force this divergence is found to diverge as $`N\mathrm{}`$.

## IV Unstable case

Now consider an unstable warm beam ($`f^{}(0)>0`$). Line N1 (resp. N2) of Fig. 3 displays $`\mathrm{ln}(\omega _\mathrm{b}(t)/\gamma _\mathrm{r})`$ versus time for (3)-(4) with a CG distribution with $`N=128000`$ (resp. 512000) and $`\gamma _\mathrm{r}=0.08`$. Line V1 (resp. V2) shows $`\mathrm{ln}(\omega _\mathrm{b}(t)/\gamma _\mathrm{r})`$ versus $`\gamma _\mathrm{r}t`$ for the kinetic system and the same initial distribution with a $`32\times 128`$ (resp. $`256\times 1024`$) grid in $`(x,p)`$ space. All four lines exhibit the same initial exponential growth of linear theory with less than 1% error on the growth rate. Saturation occurs for $`\omega _\mathrm{b}/\gamma _\mathrm{r}3.1`$ . Lines N1 and V1 do not superpose beyond the first trapping oscillation after saturation. Note that, in our system, oscillating saturation does not excite sideband Langmuir waves as our hamiltonian incorporates only a single wave, not a spectrum.

After the first trapping oscillation, kinetic simulations exhibit a second growth at a rate controlled by mesh size. Line V2 suggests that a kinetic approach would predict a level close to the trapping saturation level on a time scale awarded by reasonable integration time. This level is fairly below the equilibrium $`V_{\mathrm{th}}`$ predicted by a gibbsian approach  ; such pathological relaxation properties in the $`N\mathrm{}`$ limit seem common to mean-field long-range models . Both kinetic simulations also exhibit a strong damping of trapping oscillations, which disappear after a few oscillations, whereas finite-$`N`$ simulations show persistent trapping oscillations.

One could expect that finite-$`N`$ effects would mainly damp these oscillations, so that the wave amplitude reaches a plateau. Actually, we observe persistent oscillations for all $`N`$, and the wave amplitude slowly grows further, whereas the velocity distribution function flattens over wider intervals of velocity.

This spreading of particles is due to separatrix crossings, i.e. successive trapping and detrapping by the wave . Indeed, when the wave amplitude grows (during its pulsation), it captures particles with nearby velocity, i.e. with a relative velocity $`\mathrm{\Delta }v_{\mathrm{in}}\pm \sqrt{8|V|}`$ ; the trapped particles start bouncing in the wave potential well. When the wave amplitude decreases, particles are released, but if they experienced only half a bouncing period, they are released with a relative velocity (with respect to the wave) opposite to their initial one, i.e. $`\mathrm{\Delta }v_{\mathrm{out}}\mathrm{\Delta }v_{\mathrm{in}}`$. Now notice that a particle which has just been trapped would oscillate at a longer period than the nominal bouncing period (namely the one deep in the potential). Moreover, if the recently trapped particle had just adiabatic motion in the well, it would have to recross the separatrix when the resonance would enclose the same area as at its trapping . Thus one expects the particle to be unable to complete a full bounce, and the fraction of particles for which $`\mathrm{\Delta }v_{\mathrm{out}}\mathrm{\Delta }v_{\mathrm{in}}`$ is significant.

During this particle spreading process in $`(x,p)`$ space, the wave pulsation is maintained by the bunch of particles which were initially trapped, and are deep enough in the potential well to remain trapped over a whole bouncing period. These particles form a macroparticle, as is best seen in the case of a cold beam . Note that, over long times, the macroparticle must slowly spread in the wave resonance, following two processes. One acts if the trapped particle motion is regular : the trapped motions are anisochronous, i.e. have different periods (only the harmonic oscillator has isochronous oscillations). The other one works if the motion is chaotic : nearby trajectories diverge due to chaos. Both processes contribute to the smoothing of the particle distribution for long times, but over much longer times than those over which we follow the system evolution and observe the wave modulation.

This second growth after the first trapping saturation depends on the shape of the initial distribution function. In Fig. 3(b), line N2 is the same as in Fig. 3(a), computed over a longer duration, and line N3 corresponds to $`N=64000`$ with the TL distribution of Fig. 1. Although N3 corresponds to 8 times fewer particles than N2, the final level reached at the end of the simulation is lower. In the second growth, particles are transported further in velocity, so that the plateau in $`f(p)`$ broadens with time. As the wave grows, it can trap particles with initial velocity further away from its phase velocity. Since the TL distribution reaches its maximum at $`v1.06`$ and decreases significantly beyond this velocity (while CG is still growing for larger $`v`$), fewer particles (with TL than with CG) can give momentum to the wave when being trapped ($`𝖯`$ is conserved) ; hence the second growth is slower for the TL distribution.

We followed the evolution of the wave amplitude for N3 up to $`\gamma _\mathrm{r}t=1750`$ : starting from the first trapping saturation level ($`0.4V_{\mathrm{th}}`$), fluctuations persist with a growth rate that slowly decreases as we reach $`0.78V_{\mathrm{th}}`$ at the end of the computation. Line N4 of Fig. 3 corresponds to the TL distribution with 2048000 particles and shows persistent oscillations with approximately the same amplitude as for $`N=64000`$.

## V Conclusion

These observations clearly indicate that the kinetic models are an idealization and do not contain all the intricate behavior of a discrete particles system. Now, we must also admit that the kinetic simulation schemes do not exactly reproduce the analytic implications of the kinetic equation. It is then legitimate to ask whether the numerical implementation of the kinetic equations reproduce the difference between the finite-$`N`$ dynamics and the kinetic theory.

A basic property of the collisionless kinetic equation is that it transports the distribution function $`f(x,p)`$ along the particle trajectories (or characteristic lines in $`(x,p)`$ space). As long as the kinetic calculation of $`f`$ is accurate, one expects the kinetic scheme to follow closely the $`N`$-particle dynamics too. However, the kinetic scheme is bound to depart from the analytic predictions of the kinetic equation, because the (chaotic or anisochronous) separation of particle trajectories implies that constant-$`f`$ contours eventually evolve into complex, interleaved shapes. This filamentation is smoothed by numerical partial differential equation integrators, while $`N`$-body dynamics follows the particles more realistically, sustaining the trapping oscillations. Hence both types of dynamics will depart from each other when filamentation reaches scales below the semi-lagrangian kinetic code grid mesh.

The onset of filamentation is easily evidenced in kinetic simulations. Indeed, whereas the kinetic equation analytically preserves the 2-entropy $`(1f)f𝑑x𝑑p`$, numerical schemes increase entropy significantly when constant-$`f`$ contours form filaments in $`(x,p)`$-space . As this is also the time at which trapping oscillations are found to damp in our simulations, it appears that vlasovian simulations must be considered with caution from that time on – and it turns out that it is also the time from which the second growth starts.

In summary, discussing the basic propagation of a single electrostatic wave in a warm plasma, we presented finite-$`N`$ effects which do not merely result from numerical errors and elude a kinetic simulation approach. Their understanding depends crucially on the dynamics in phase space. The sensitive dependence of microscopic evolution to the fine structure of the initial particle distribution in phase space  implies that the interplay between limits $`t\mathrm{}`$ and $`N\mathrm{}`$ requires some caution. Somewhat paradoxically, refining the grid for the Vlasov simulations does not solve this problem.

The driving process in the system evolution is separatrix crossing, which requires a geometric approach to the system dynamics. Further work in this direction  will also shed new light on the foundations of common approximations, such as replacing original dynamics (1)-(2) by coupled stochastic equations, in which particles undergo noisy transport.

## VI Acknowledgments

The authors thank D.F. Escande for fruitful discussions, and J.R. Cary and I. Doxas for computational assistance. MCF and DG were supported by the French Ministère de la Recherche. Computer use at Institut Méditerranéen de Technologie and IDRIS was granted by Région Provence-Alpes-Côte d’Azur and CNRS. This work is part of the european network Stability and universality in classical mechanics and CNRS GdR Systèmes de particules chargées (SParCh).


text:
# STATUS OF MAGIC AND EXPECTED CAPABILITIES

## 1 Introduction

Imaging Air Cherenkov Telescopes (IACTs) have established in recent years the existence of a few galactic and extragalactic sources above 300 GeV $`^\mathrm{?}`$ (the so-called Very High Energy gamma ray range). New technical developments currently allow to reduce the threshold of this kind of telescopes to energies around 30 GeV, thus overlapping those typical of $`\gamma `$-ray satellites like EGRET or the forthcoming GLAST and AGILE. Here we describe a new generation 17 m diameter IACT dubbed MAGIC $`^{\mathrm{?},\mathrm{?}}`$ (Major Atmospheric Gamma Imaging Telescope) with a threshold of 30 GeV in its first development phase (phase I) and around 10 GeV in phase II.

## 2 The general design

MAGIC is a second generation IACT with a number of improvements and innovative elements to reduce the energy threshold and the different backgrounds.

### 2.1 The reflector

Thanks to its 17 m diameter disk, MAGIC will collect three times more light than a typical first generation 10 m diameter IACT. The shape of the reflector is parabolic, hence isochronous: this improves the light background reduction. A carbon fiber mount, weighing less than 10 tons, is instrumental in rapid repositioning the telescope in the search for GRBs. This material makes the disk also hard to deform (sagging $`<`$ 3mm). In addition, an active mirror control system $`^\mathrm{?}`$ provides maximum shape stability. The reflector is tessellated, with 50 cm $`\times `$ 50 cm all-alluminium, diamond-turned, quartz-coated, mirror elements $`^\mathrm{?}`$ of 85% reflectivity in the 300-650 nm wavelength region.

### 2.2 The camera

New compact (2.5 cm, 0.1 diameter) photomultipliers $`^\mathrm{?}`$ with approximately 20% quantum efficiency (QE) in the 300-500 nm range and minimal time dispersion shall be used in phase I. In order to minimize the impact of DC like night sky background (NSB) and moonlight the dynodes were reduced to six and operated with a total gain of only 15000. The camera of 4 diameter is composed of 600 PMTs with coarser sampling in the outer region.

In phase II the camera will be provided with Hybrid Photo Detectors of $`<`$QE$`>`$ $``$ 40% in an wider wavelength range of 330-650 nm. Extending the light detection range to red wavelengths is particularly useful for high zenith angle observations where most of the short wavelength photons are absorbed. Finally, in its phase III, MAGIC aims at equipping the camera with Avalanche Photo Diodes, provided with an extremely high $`<`$QE$`>`$ $``$ 80% (330-680 nm).

### 2.3 The readout

Cherenkov light pulses are typically very short (in the order of 1-2 ns). However as the pulse shape contains important physical information and can be used for NSB and hadron reduction, it is crucial to register it faithfully. MAGIC endeavors to do so by applying a number of improvements. To begin with the signal is transported over optical fiber $`^\mathrm{?}`$. Apart from negligibly dispersing the pulse, optical transport reduces cable weight and enables optical decoupling and noise immunity.

In addition, the pulse is digitized using very fast (300 MHz, and 1 GHz in a second phase) 8 bit FlashADCs $`^\mathrm{?}`$. The maximum sustained event rate will be 1 kHz (with zero deadtime). Dual ranging of the analog signals at the FADC inputs extends the range of the FADCs to about 70-80 dB.

### 2.4 Sensitivity and performance of MAGIC phase I

The telescope will have a threshold of about 30 GeV (peak of the differential flux) for phase I. The gamma-ray sensitivity ranges from $``$10-10 at around 10 GeV (set by the cosmic electron background $`^\mathrm{?}`$) down to 8 $``$ 10-12 cm-2 s-1 at 1 TeV (determined by the cosmic hadron background). The collection area flattens at around 105 m2 above 100 GeV. Close to the threshold the energy resolution will be around 30% and improve to around 10% at 1 TeV $`^\mathrm{?}`$.

In the presence of moonlight, the threshold has to be increased to 60-100 GeV (half-moon, $`>`$ 30 away). Operation under moonlight significantly increases the potential of the telescope to detect GRBs and follow up other transient sources.

MAGIC is under construction in the Roque de los Muchachos Observatory at La Palma (Canary Islands, Spain). First light is expected for summer 2001 and first physics results may be available on 2002.

## 3 Prospects for Astrophysics and Fundamental Physics

Here are briefly described some of the astrophysics, cosmology and fundamental physics questions MAGIC is expected to address.

### 3.1 Active Galactic Nuclei

AGNs exhibit their highest variability in X-rays and $`\gamma `$-rays. Multiwavelength campaigns including IACTs have already proved central to understand the physics of the source, specially because GeV $`\gamma `$-rays probably come from very close to the supermassive black hole.

### 3.2 The cosmological $`\gamma `$-ray horizon

High energy $`\gamma `$-rays are absorbed in the infrared background ($`\gamma _{VHE}+\gamma _{IR}e^{}+e^+`$). This means that there is a maximum observable redshift, zhor(E), which is usually referred to as $`\gamma `$-ray horizon. Only AGNs below this redshift are observable. Due to its extremely low energy threshold, MAGIC will see the bulk of the cosmological AGNs. A large number of detections shall allow to characterize the AGN population as a function of distance and determine if $`\gamma `$-rays are indeed only externally absorbed or there is some internal absorption. If a large enough number of AGNs are detected, a measure of the infrared background may also be attempted.

### 3.3 The diffuse extragalactic $`\gamma `$-ray background

EGRET has measured the $`\gamma `$-ray background up to 100 GeV. The general opinion is that it is due to AGNs. But it has been proved $`^\mathrm{?}`$ that AGNs only contribute 25% of the background above 100 MeV. The remaining background may come from FSRQs and BL Lacs $`^\mathrm{?}`$. In this case MAGIC may detect more than 100 more sources than EGRET thanks to its higher sensitivity. The background could alternatively arise from topological defects (which may in turn be also the origin of the cosmic rays above 1019 eV).

### 3.4 Pulsars

Three pulsars are among the strongest EGRET sources. Observations at energies above 10 GeV may help clarify the mechanisms producing the $`\gamma `$-rays. The outer gap or polar cap models differ in the maximum energy attainable for the pulse emission. In addition, and owing to its improved angular resolution, MAGIC can clarify if many of the $``$ 100 EGRET unidentified sources are associated to pulsars. A special pulsar trigger has been developed to attain a low energy threshold around 10 GeV based on pulsar timing analysis.

### 3.5 The origin of cosmic rays: shell-type SNRs

It is widely believed that CRs are produced in SNR blast shocks. X-rays and $`\gamma `$-rays are essential probes, since, unlike CRs, they point to the source. There is good evidence of electron acceleration. Synchrotron X-rays and inverse Compton $`\gamma _{VHE}`$’s produced by electrons have been observed in SN-1006. However no firm evidence has been found for hadron acceleration. Observation of $`\gamma `$-rays coming from dense clouds serving as targets for cosmic rays accelerated in the shell (through $`p,Fe+N\pi ^{}s`$ and $`\pi ^0\gamma \gamma `$) would provide an unambiguous proof of the presence of hadrons in the shell.

### 3.6 Gamma Ray Bursts

The mechanism producing GRBs have been not resolved yet. MAGIC counts on its huge collection area to detect fast transients with high statistics. The telescope has actually been optimized for GRB searches by ways of its fast repositioning capability. Its potential for GRB detection is illustrated by the very high rates which would have been expected for some typical GRBs. For example, GRB-910503, with a maximum measured energy of 10 GeV and lasting for 84 s, would have produced a rate of 6$``$104 Hz in MAGIC, whilst GRB-930131, reaching at least 1.2 GeV and lasting for 100 s, would have produced around 2$``$104 Hz. Such strong signals allow us to study the details of the light curve and further characterize the source structure.

### 3.7 Invariance of the speed of light

Several quantum gravity models predict an energy dispersion of the speed of light $`^\mathrm{?}`$. This results in a arrival time difference of high energy photons as respect to low energy ones which depends on particle energy and quantum gravity energy scale. By detecting sub-second time differences between e.g. keV and 30-300 GeV photons for z$``$0.1-1, energy scales in the order of 1018-1020 GeV are well measurable. A similar analysis may also be applied to pulsar pulsations.

### 3.8 Cold dark matter

One of the most plausible candidates for dark matter is the neutralino ($`\chi `$). The particle physics lower limit for mχ is around 30-50 GeV and further evidence suggests that its mass should be below 1 TeV. An annihilation line ($`\chi \chi \gamma \gamma `$) may thus be expected from the center of our Galaxy $`^\mathrm{?}`$ in the most sensitive range of the telescope.

## Acknowledgments

I am deeply grateful to my colleagues from the MAGIC collaboration for providing extensive information. In particular I would like to thank N. Magnussen, E. Lorenz and R. Mirzoyan for their kind help. The support of the German BMBF and Spanish CICYT is also acknowledged.

## References