format: text

text:
# CHANDRA Observations of “The Antennae” Galaxies ( NGC4038/39): II. Detection and Analysis of Galaxian X-ray sources

## 1 Introduction

The first imaging X-ray observations of nearby galaxies showed that a large fraction of their X-ray emission is due to a population of discrete X-ray sources: Supernova remnants (SNRs) and X-ray binaries (XRBs) (e.g. Fabbiano 1989). However, the properties of the X-ray source population and its dependence on galaxian parameters (such as star-formation) could not be investigated before the deployment of Chandra because of the poor sensitivity and spatial resolution of the pre-Chandra X-ray observatories (e.g. Fabbiano 1995).

At a distance of 29 Mpc ($`\mathrm{H}_\mathrm{O}=50`$ km s-1 Mpc-1) the Antennae (NGC 4038/39), the prototypical merging galaxies (Toomre & Toomre 1972), are a prime target for the study of the X-ray source population in star-forming galaxies. The Antennae have been observed in the X-rays with almost every imaging X-ray telescope. The first Einstein observations detected a few discrete sources as well as diffuse, possibly soft, emission associated with this system (Fabbiano et al. 1982, Fabbiano & Trinchieri 1983). Subsequent observations with ROSAT HRI (with a resolution of $`5^{\prime \prime }`$) resolved 12 discrete sources or emission regions and yielded for the first time indications that some of them might be Ultra Luminous X-ray sources (ULXs), emitting well in excess of the Eddington limit for an accreting neutron star (Fabbiano et al. 1997). One of these sources was found to be variable between two observations obtained one year apart. Observations with the ROSAT PSPC and ASCA (Read et al. 1995, Sansom et al. 1996) did not provide significant information on the population of the discrete sources, due to their poor spatial resolution. However, the ASCA observations showed that there is a significant hard X-ray component (which could derive from a population of X-ray binaries and/or from hot thermal gas), as well as a softer gaseous emission component.

NGC4038/39 were observed with Chandra ACIS-S, and the overall results of the observations were reported by Fabbiano, Zezas & Murray (2001) (Paper I). As reported there, a large population of extremely luminous ($`\mathrm{L}_\mathrm{X}10^{38}10^{40}`$ erg s-1) point-like sources is detected, along with a diffuse component. The two components contribute equally to the total X-ray luminosity of the system. Here, we describe the detailed analysis of the X-ray properties of the discrete sources. These results will be elaborated upon and discussed in the companion paper (Paper III; Zezas, Fabbiano, Rots & Murray 2002), where we also present a comparison with multiwavelength properties of the emission regions. The derivation of the X-ray luminosity function (XLF) of the X-ray sources of the Antennae Galaxies and a comparison with the XLFs of other galaxies is presented in Paper IV (Zezas & Fabbiano 2002).

Throughout this paper we use a distance of 29 Mpc ($`\mathrm{H}_0=50`$ km s-1 Mpc-1). For $`\mathrm{H}_0=75`$ km s-1 Mpc-1, the distance becomes 19.3 Mpc and all the cited luminosities will be a factor of 2.25 lower. At a distance of 29 Mpc, 1 arcsec corresponds to a physical distance of $`140`$ pc (at 19 Mpc it corresponds to 92 pc). The quoted errors in the spectral fits are at the 90% confidence level for one interesting parameter, unless otherwise stated.

## 2 Observations and Data Analysis

This study is based on observations (OBSID 315) obtained with the ACIS-S camera (Garmire et al. 1997) on board the Chandra X-ray Observatory (Weisskopf et al. 2000). The Antennae were observed for a total of 72 ks on December 1, 1999. The first results from these observations have been reported in Paper I. The initial processing of the raw data was performed by the Chandra X-ray Center (CXC), with the reprocessing pipeline software (version R4CU5UPD6.5), using the updated calibration files as of September 2000. In the following analysis we mainly used the CIAO package 111http://cxc.harvard.edu/ciao. The geral data analysis techniques are presented in the Chandra Data Analysis Threads 222http://cxc.harvard.edu/ciao/documents\_threads.html. For our analysis we used the Level-2 event files which include only events of ASCA grades 0,2,3,4,6 and status 0. We searched the data for periods of enhanced background radiation by extracting a light curve from a large source-free region of the ACIS-S3 chip. The background as measured from the event2 file (prior to any processing and including events of all energies) was found to be fairly constant, at a level of $`0.456\pm 0.004`$ counts/sec in an area of 7.06 arcmin2. This is roughly half of the quiescence background level reported in the Chandra Proposers’ Observatory Guide. We screened the raw data to exclude events with energies outside the (0.3 - 10) keV band. The 0.3 keV cut-off is chosen because of the on-board event rejection, whereas the 10.0 keV cutoff is set because the telescope’s High Resolution Mirror Assembly (HRMA; Van Speybroeck et al. 1997) effective area is effectively zero above 10 keV. The ACIS-S3 chip which was in the focal point of the HRMA has two “hot columns” at the borders of each node, which we excluded (see Paper I).

In order to check the absolute astrometry of these observations, we searched the U.S. Naval Observatory (USNO) star catalogue for stars within 4 arcmin from the galaxy. We found three stars with X-ray counterparts. The offsets between the optical and the X-ray sources were not systematic and less than $`0.5\mathrm{}`$. Therefore we conclude that the absolute astrometry of these observations is good to within $`0.5^{\prime \prime }`$, in accordance with estimates from the calibration team (Aldcroft et al. 2001)333http://cxc.harvard.edu/cal/ASPECT/celmon/.

### 2.1  Source Detection

After the initial processing we created images in the (0.3-2.0) keV, (2.0-4.0) keV, and (4.0-10.0) keV bands. The lower energy band allows direct comparison with the ROSAT results. We decided to divide the (2.0-10.0) keV band in two intervals at 4.0 keV, in order to obtain more information on the spectral shape of the X-ray emission. The 4.0 keV boundary was chosen because at higher energies the effective area of the ACIS-S3 chip begins to drop and we wanted to retain a relatively good signal to noise ratio (S/N) for the objects detected in the hard band. We smoothed each image using the adaptive smoothing algorithm implemented in the csmooth tool of the CIAO v2.0 suite. This algorithm adjusts the size of the smoothing kernel (in our case a two-dimensional Gaussian) in order to maintain a uniform S/N ratio over the smoothed image. For these images the S/N ratio was set to be between 3 and 5. The smoothing scales range between 0.5 arcsec (for strong point sources) and 256 arcsec (for very low surface brightness diffuse emission). The majority of the image is smoothed with a Gaussian with a FWHM of $`15`$ pixels. The images in the three energy bands, together with a full band unsmoothed image, are presented in Fig. 1.

We searched the unsmoothed images in the three energy bands for discrete sources, using both the sliding cell and wavelet algorithms implemented in the CIAO celldetect and wavdetect tools, respectively. We ran celldetect with a cell corresponding to 90% of the encircled energy of a point source which varied with off-axis angle. The cell size ranges between 6 pixels for an on-axis source to 12 pixels for a source at 2 arcmin off-axis (which is the maximum in the case of the Antennae). The background was estimated from an annulus around the source defined to cover the same area as the source cell. Because the detection cell is not adjustable to the size of the source, the effectiveness of this algorithm is limited to unconfused regions and point-like sources. It also strongly depends on the local background, which may well include other nearby sources. The 3$`\sigma `$ best fit ellipses to the spatial distribution of the events for each source are plotted as red ellipses in Fig. 1 on the raw (0.3-10.0) keV band image.

As the X-ray emission of the Antennae is very patchy, with crowded regions of sources embedded in diffuse emission, we also used the wavedect source detection tool which performs better in these “difficult” conditions (Freeman et al. 2001a). The advantage of this tool is that it parameterizes the shape of a point-like source with the “Mexican Hat” (MH) function, consisting of a narrow core and relatively extended wings. This function can be spatially rescaled to match the extent of each source, allowing one to obtain information on the spatial properties of the detected sources and to detect and accurately measure fluxes for sources at various spatial scales. In our analysis, we searched for sources on scales of 1, 2, 4, 8 and 16 pixels. These scales correspond to the radius R of the core of the MH function (its Full Width at Zero Intensity is $`2\sqrt{2}\times R)`$. The determination of the background was done from a background map produced by counts detected in the wings of the MH function (see the Detect guide for more details; Dobrzycki et al. 2000).

In order to establish whether a source is point-like or extended and to estimate its total emission (corrected for the Point Spread Function; PSF), wavedetect compares the scale-length which matches best the shape of the source with the size of the PSF at its position. Below 4.0 keV we used models of the PSF at 1.49 keV as most of the photons are detected between 1.0-2.0 keV, where the ACIS-S effective area peaks. Moreover, at this energy the HRMA PSF is most accurately calibrated. For higher energies we used the calibration at 4.5 keV since above this energy the effective area of the detector falls rapidly and even for the hardest sources most of the photons are detected below 6.0 keV. The limiting chance detection probability was set to $`10^6`$, which corresponds to $`0.5`$ spurious sources in the searched area (Dobrzycki et al. 2000). Fig. 1a shows the 3$`\sigma `$ ellipse fit to each source found by wavdetect as white ellipses.

The merged source list including results from both celldetect and wavdetect results in the (0.3-10.0) keV band is presented in Table 1. For completeness we also include in this table sources detected only in the soft or medium band. Column (1) gives a number indentifying the sources internally in this paper (sources 1-49), sources 1c, 2c and 3c at the bottom of the table are sources detected only by  celldetect, Column (2) gives the CXO source name, Columns (3) and (4) give the RA and Dec (J2000), Column (5) gives the background subtracted source counts together with the $`1\sigma `$ statistical error following the Gehrels approximation for low count statistics (Gehrels 1986), Columns (6) and (7) list the background contribution to the total source counts in the detection region of each source and the S/N ratio of the detection, Column (8) gives the extent of each source in terms of the fraction of the PSF size at its position on the detector (defined as the geometrical mean of the minor and major axis of the $`1\sigma `$ ellipse divided by the 39% encircled energy radius of a point source in the soft (0.3-2.0 keV band)). Columns (9) to (14) give similar information as Columns (3)-(9) but for the sources detected with celldetect. Columns (15) and (16) give the logarithm of the observed and absorption-corrected (0.1-10.0) keV band luminosity respectively for each source in erg s-1 (assuming a 5 keV bremsstrahlung model and galactic line of sight absorption $`\mathrm{N}_\mathrm{H}=3.4\times 10^{20}\mathrm{cm}^2`$; Stark et al. 1992, and using the AO2 release of the ACIS-S3 effective area). For sources detected with both algorithms these luminosities are based on the count rate measured with wavdetect, since it estimates more accurately the source intensity and the local background. For sources detected only in the soft band we estimate the total band luminosity assuming the same spectrum as for sources detected in the full band. However, for sources detected only the medium band we also use a 5 keV thermal bremsstrahlung model but with a column density of $`2\times 10^{21}`$ cm-2 based on spectral fitting results (section 2.4). Table 1 shows that for the sources detected with both celldetect and wavdetect the derived parameters agree fairly well: the measured net source counts are typically consistent within 30%. The greatest differences are for resolved sources which are slighlty extended or confused and therefore not well modeled with the detection cell used by celldetect.

Fig. 1a shows the position of each source detected with wavdetect or celldetect (marked by the 3$`\sigma `$ ellipse) overlaid on the full-band raw data. The sources are numbered following the convention of Table 1. As it was expected, wavdetect is more efficient than celldetect, in detecting and separating close sources. This figure also shows that the three sources detected with celldetect alone and listed at the bottom of Table 1 are local enhancements of the diffuse emission which are rejected by the source selection algorithm of wavdetect, and therefore we do not include them in our following analysis.

In Table 2 we present the results from the wavdetect run for each energy band separately. Column (1) has the source identification number used in this paper, columns (2) to (5) give the net source counts (with 1$`\sigma `$ statistical error), the background counts in the source area, the S/N of the source (in $`\sigma `$) and its extent (defined in as in Table 1). Columns (6)-(9) and (10)-(13) give the same information for the medium (2.0-4.0 keV) and hard (4.0-10.0 keV) bands. In both Tables 1 and 2 we included only sources detected at the 3$`\sigma `$ level above the local background and contained within the optical outline of the Antennae galaxies. Figs. 1b,c,d show the position of each source (marked by their $`3\sigma `$ best fit ellipse in the full band), following the numbering convention of Table 1, overlaid on adaptively smoothed soft, medium and hard band images.

We detect a total of 49 sources, possibly associated with the Antennae. The limiting luminosity is $`10^{38}`$ erg s-1 ($`\mathrm{H}_0=50`$ km s-1 Mpc-1; $`5\times 10^{37}`$ erg s-1, for $`\mathrm{H}_0=75`$ km s-1 Mpc-1). 31 sources (38 for $`\mathrm{H}_0=75`$ km s-1 Mpc-1) have luminosities below $`10^{39}`$ erg s-1 while 18 (11) have higher luminosities, so, as reported in Paper I, they are Ultraluminous X-ray Sources (ULXs), emitting well in excess of the Eddington limit of a neutron star accretion binary ($`3\times 10^{38}\mathrm{ergs}\mathrm{s}^1`$). The slight difference between these results and those presented in Paper I are due to the use of more appropriate effective areas and unabsorbed fluxes for the determination of the luminosities. To estimate the number of sources not physically associated with the Antennae we performed a simple calculation using both the logN$``$logS relation of Giacconi et al. (2000) and source counts of background sources from this observation. We found that only 2.6 serendipitous sources are expected in the area of the Antennae, down to a luminosity level of $`1\times 10^{38}`$ erg s-1 (Paper I). This number reduces to 0.6 sources for sources brighter than $`10^{39}`$ erg s-1. We therefore conclude that most of the detected sources, and especially the higher luminosity ones, do belong to NGC4038/39.

### 2.2 Spatial Analysis

Determining if a source is point-like or extended depends on the S/N of the source, because this limits the accuracy of our measurements of the spatial distribution of source counts. We used 100 counts as a threshold for our investigation of source extent. A more detailed analysis of the extent of fainter sources will be presented in Paper V (Fabbiano et al. 2002, in prep). Due to the very high spatial resolution of Chandra ,sources with fewer counts are detected over typically a $`2\times 2`$ or $`3\times 3`$ pixel cell, which is too small to allow any spatial profile fitting. Spatial fitting is also difficult because the ACIS undersamples the HRMA PSF (the 50% encircled energy radius of HRMA at the focal aim-point is $`0.3^{\prime \prime }`$ while the ACIS pixel size is $`0.49^{\prime \prime }`$). For this reason we consider all the sources with fewer than 100 counts and small source ellipses in the (0.1-10.0) keV band as point-like, without doing any profile fitting. Instead we obtained radial profiles for all sources detected with more than 100 counts. We also derived radial profile for four sources detected with fewer than 100 counts, but having $`3\sigma `$ best fit ellipses with major axis larger than 8 pixels (sources 5,6,7 and 10), in order to investigate if they have a point-like conponent.

For these sources we created models of the PSF, appropriate to the position of each source on the detector, using the CIAO mkpsf tool which generates a PSF at a given energy and off-axis angle, based on pre-flight calibration. As the best calibration of the PSF of ACIS is at the energy of 1.459 keV, where also the sensitivity of ACIS-S3 is optimized, we extracted our model PSFs at this energy and we compared them with the respective source profiles in the 0.3-2.0 keV band. Although these model PSFs do not take into account the telescope motion during the observation, the residuals after the application of the aspect solution are usually less than $`0.75^{\prime \prime }`$ ($`1\sigma `$ level) 444http://asc.harvard.edu/mta/ASPECT/celmon/index.html, and therefore we can directly compare these models with the actual observations. After calculating the model PSFs, we rescaled them to the total number of photons detected in each source. The PSF becomes asymmetric at large off-axis angles ($`\theta 1.5^{}`$). This mainly affects the sources associated with NGC 4039 (Southern galaxy) since the optical axis was roughly on the position of the nucleus of NGC4038 (Northern galaxy). For this reason we used elliptical annuli instead of circular annuli for the profile extraction. The PSF consists of a narrow core and more extended wings, which rotate depending on the position of the source relative to the aimpoint, thus making the determination of the position angle of the major axis non-trivial. For our analysis we determined the angle of the major axis directly from the model PSF.

We used the CIAO application Sherpa (Freeman et al. 2001b) to fit the radial profiles of the sources using as models the radial profiles of the simulated PSFs (extracted using exactly the same regions as for the source profiles). In the cases which show residual wings, due to incomplete subtraction of the source background we added a constant to the PSF. For each source, identified in Column (1), Table 3 presents the results of these fits: Column (2) gives the $`\chi ^2`$ together with the number of degrees of freedom for each fit with the model PSF, and Column (3) gives the same information for fits with the PSF plus a constant background. Figure 2 shows the results of the spectral fits with the PSF plus background model for 3 examples of point-like sources and for all the extended sources. The points display the source profile, whereas the solid line represents the model PSF. The bottom panel of each plot shows the fit residuals in $`\sigma `$. The fits with the PSF alone are in most cases very poor and the inclusion of a constant improves them in the 99% confidence level.

In order to obtain a measurement of the physical scale of each source we also fitted both the source and the model PSFs with a Gaussian. The results of these fits are also presented in Table 3, where Columns (4) and (6) give the FWHM of the source and the model PSF profile in ACIS-S pixels (with $`1\sigma `$ errors for one interesting parameter), and Columns (5) and (7) give the $`\chi ^2`$ and number of d.o.f for the source and the model PSF respectively. For some sources the resulting reduced $`\chi ^2`$ is vey high. This is due to residual local diffuse emission and also to the wings of the PSF which are not well modeled by a single Gaussian. For this reason we fitted the profiles with a Gaussian and a constant. These results are presented in Columns (8)-(11). In Column (12) we give the deconvolved FWHM of the Gaussian profile of the extended sources in arcseconds and parsecs (assuming a distance of 29 Mpc), point like sources are marked as “p”. In general the latter models gave a statistically significant improvement in the fit based on an F-test with 1 additional parameter. We find that all but 6 sources for which we measured the profiles are consistent with the model PSFs and therefore we consider them as point-like. The FWHM of the Gaussian used to fit the profiles of sources 5, 6, 7, 10, 24 and 29 are larger than those of the PSF, above the $`3\sigma `$ confidence level (source 10 has a weak point-like core but very extended wings and we consider it as extended). This clearly suggests that these sources have an extended component. Their extent ranges between $`3^{\prime \prime }`$ and $`15^{\prime \prime }`$ corresponding to a physical scale of $`400`$ pc to $`2.1`$ kpc, as measured by the FWHM of the best fit Gaussian. Sources 25 and 29 correspond to the Northern and Southern nuclear regions respectively, as determined by their radio positions (Neff & Ulvestand, 2001).

### 2.3 Timing Analysis

Because of the low count rates of the sources in the Antennae, it is very hard to search for short term variability. We extracted light-curves for the 5 brightest sources (net count rate $`>0.001\mathrm{counts}/\mathrm{sec})`$ binned with bin sizes of 500 s, 1000 s and 5000 s and we compared them with the standard deviation of the mean count rate for each light-curve. In all cases the data points were consistent with the respective mean values at the 99% confidence level. We also split the observation in five exposures of 15 ks. Again no variability was detected.

We also searched for variability in all the sources by comparing the cumulative distribution of the photon arrival times for each one with the photon arrival times of events in a large source free background region around the galaxy. In order to minimize the background we only used events in the (0.3-7.0) keV band, because the particle background dominates above 7 keV. The comparison was performed using the Kolmogorov-Smirnov test (KS test). This test showed that for the vast majority of the sources the two distributions are consistent with being drawn form the same parent population, i.e., the sources are not variable at the 99% significance level. Only for sources 14 and 44 we detect variability at a significance level of 99.2% and 99.7%, respectively. Source 14 is a relatively bright source in the southern part of the galaxy. A light-curve of this source in bins of 15ks suggests a decline of its flux, although because of the large error bars a constant still gives an acceptable fit. Source 44 is the third brightest source in the Antennae. A light curve in bins of 15ks suggests that there is a flare in the second bin, but again because of the low S/N a constant gives a good fit. Fig. 3 shows the cumulative distributions of the two sources (solid lines) compared with the distributions of the background (dashed lines). A similar test in the (2.0-10.0) keV band showed no variability signs for any sources. This could be caused by the low number of counts for most sources above 2 keV.

In order to search for long term variability we compared our data with the ROSAT HRI observations of the Antennae (Fabbiano et al. 1997). Prior to Chandra these are the highest angular resolution observations in existence. In order to directly compare our data with the HRI measurements, we concentrated on isolated sources (X-3, X-4, X-8, X-11, X-12, X-13 from Fabbiano et al. 1997). We also used exactly the same regions used to measure the HRI luminosities; we measured background from regions outside the galaxy. In order to match the HRI band we used the (0.3-2.0) keV Chandra images to measure the count rates of the sources. Then we converted both the HRI and the Chandra count rates to luminosities in the (0.1-2.5) keV band assuming a 5 keV thermal bremsstrahlung model with Galactic line-of-sight absorption ($`3.4\times 10^{20}`$ cm-2). These results are given in Table 4, where Column (1) gives the HRI source, Column (2) the Chandra source(s) encompassed by the extraction region, Column (3) the HRI luminosity, Column (4) the Chandra luminosity, and Column (5) the significance (in $`\sigma `$) of the difference between the ROSAT and Chandra measurements. These results suggest that sources 16, 42 and the combination of sources 44 and 46 show variability or a factor of $`2`$ in timescales of a few years.

### 2.4 Spectral Analysis

As mentioned in the introduction, there are two major components in the galaxian source populations: X-ray binaries and SNRs. X-ray spectra are a very powerfull tool for distinguishing between these two types of sources. Moreover, from the parameters of detailed spectral fitting it is possible to distinguish between different types of X-ray binaries (e.g. binaries in different states, pulsar from black-hole binaries, see van Paradijs 1999).

We extracted spectra for all the detected sources, using extraction regions defined to include as many of the source photons as possible, but at the same time minimizing contamination from nearby sources and background. The background region was usually a source-free circular or elliptical annulus surrounding each source, in order to take into account the spatial variations of the diffuse emission and to minimize effects related to the spatial variations of the CCD response. Any sources encompassed by the background region were excluded. The regions used to extract the source and background spectra are presented in Fig. 4a and 4b. In order to take into account spatial variations of the detector gain, the spectra were extracted in PI (Pulse Invariant; gain-corrected space.

For each spectrum we created response matrices and ancillary response matrices using the mkrmf and mkarf tools in the CIAO v.2.0 package with the N0001 (December 2000) release of the response calibration files (in the form of FITS Embedded Functions; FEFs). The variations of the charge transfer inefficiency (CTI) in adjacent FEF regions of the ACIS-S3 chip, within which the response is modeled and assumed to remain stable, are very small. For the case of ACIS-S3 chip which was at the focus of the telescope, these regions are 32$`\times `$32 pixels large. For point sources which encompass more than one FEF regions, we extracted response matrices from the FEF which includes most of the source photons. In order to assess the uncertainty of this approximation, we also created an RMF for the source which encompasses the largest number of FEF regions, using the calcrmf tool 555http://asc.harvard.edu . This tool calculates the response matrix for a source by combining the appropriate FEFs corresponding to the detector pixels on which source photons are detected after taking into account the aspect solution. The FEFs are then weighted by the number of counts in each region. We find that the difference in the RMFs calculated by the two different methods are below the 10% level, giving us confidence that the use of a single FEF in the calculation of the RMF does not affect our results. In the case of a few sources which fall at the boundary between two ACIS nodes we used the response matrix for the node onto which most of the events are detected. Since these sources have typically less that 200 counts, this will not affect our spectral fitting results, which are dominated by source statistics.

The spectra were fitted using the XSPEC v10.0 package (Arnaud et al. 1996). In order to use $`\chi ^2`$ statistics, we grouped the data to include at least 15 counts per spectral bin before background subtraction. In the spectral fitting, we excluded any events with energies above 10.0 keV or below 0.3 keV. Although, the calibration at energies below 1.0 keV is uncertain, we retained data in the this range because the statistical error of the data points at these energies is still significantly larger than the expected calibration uncertainty. Moreover, these data are important for constraining the absorbing column density. It is known that the ”energy scale” has a shift of $`20`$ eV 666http://asc.harvard.edu/cal/Links/Acis/acis/Cal\_projects/Energyscale\_120.html. Although this feature is not taken into account in the release of the calibration data we used, it does not affect our analysis as each spectral bin below 1.0 keV covers an energy range much larger than the energy resolution of ACIS-S3.

#### 2.4.1 Single component models

We fitted the spectra of all the sources with more than $`50`$ net counts (5 spectral bins or more) in the 0.3-10.0keV band, with simple one component models: absorbed power-law (PO) and absorbed Raymond-Smith (RS; Raymond & Smith 1977). In both cases, the $`\mathrm{N}_\mathrm{H}`$ was initialy fixed to the Galactic ($`3.4\times 10^{20}`$ cm-2) value. The RS model was suggested by the residuals around 1 keV in the PO fits of some sources, as a thermal plasma component produces the Fe-K$`\alpha `$ blend in this energy region. We then left the $`\mathrm{N}_\mathrm{H}`$ free to vary and found an improved fit in all cases at the 99% confidence level, based on an F-test for one additional parameter. The best fit $`\mathrm{N}_\mathrm{H}`$ is typically higher than the Galactic value. In all the fits with the RS model the abundance was fixed to the solar value. The results of the fits are summarized in Table 5: Column (1) gives the source identifier, Columns (2) and (3) the photon index and the $`\chi ^2`$ (with the number of degrees of freedom; dof) for the PO model with absorption fixed to the galactic. Columns (4)-(6) give the results for the fits with the PO model and free $`\mathrm{N}_\mathrm{H}`$: Col. (4) gives the photon index $`\mathrm{\Gamma }`$, Col. (5) gives the best fit absorbing column density in units of $`10^{22}`$ cm-2, and Col (6) gives the $`\chi ^2`$ with the dof. Columns (7) and (11) give the analogous results for the RS models. All errors are at the 90% confidence level for one parameter of interest. The spectra, together with the best fit PO models and the resulting residuals are presented in Figure 5. In the same figure we present $`\mathrm{\Gamma }\mathrm{N}_\mathrm{H}`$ confidence contours for these fits. The contours are at the 1,2 and 3$`\sigma `$ levels for two parameters of interest.

In most cases the single component PO and RS models give satisfactory fits. Only for 10 sources it is possible to favor one of the two models, based on an F-test. The $`\chi ^2`$ of these sources is shown in bold type in Table 5. We note that for a few sources the PO fits give relatively large best fit absorbing column densities. For some of them this is clearly indicated by the low energy cutoff of their spectra (e.g. sources 12, 24,25, 33, 35; Fig. 5). However, for other sources the high inferred absorption may be due to inadequate fitting of the continuum by the PO model as for example in the 3 sources (src 6, 10, 18) which are better fitted with an RS model, where the derived column density is significantly lower than that obtained by the PO fits. In the next section we will discuss more complex models which give in most cases lower overall absorption but usually require a heavily absorbed hard component.

In Table 6 we give the luminosities of the individual sources based on the best fit PO model with free $`\mathrm{N}_\mathrm{H}`$ except for the sources which can be modeled better or equally well with a complex model (§2.4.2) where we used the latter to determine the luminosities. In this Table Column (1) gives the source number, Column (2) gives the background subtracted count rate in the (0.3-10.0) keV band, Columns (3) and (4) give the observed flux in the soft (0.1-2.5) keV and hard (2.5-10.0) keV bands whereas Columns (5) and (6) give the flux corrected for the line of sight Galactic absorption (but keeping the remaining intrinsic absorption) and Columns (7), (8) give the flux corrected for the total absorption in the same bands. The fluxes are given in units of $`10^{14}`$ erg cm-2 s-1. Columns (9)-(14) give the logarithm of the observed and absorption corrected luminosities in the soft and hard bands (in units of erg s-1). Although these are the most accurate estimates of the source luminosities, they are not available for all the sources; for this reason throughout the paper we will use luminosities from Table 1 unless otherwise stated.

#### 2.4.2 Complex models

From the data residuals relative to the best fit power-law models (Fig. 5), one can see that there is a large number of sources with a soft excess below 2.0 keV. For this reason, the spectra of the sources with adequate numbers of counts were fitted with a double component model consisting of a PO+RS spectrum. The PO is seen through an absorber free to vary, whereas the total spectrum is observed through the Galactic line of sight $`\mathrm{N}_\mathrm{H}`$ (fixed). The results of these fits are presented in Table 7. In this Table Column (1) give the source number, column (2) and (3) give the best fit temperature and photon index for the RS and the PO models respectively and column (4) gives the column density in units of $`10^{22}\mathrm{cm}^2`$. Column (5) gives the $`\chi ^2`$ and the number of dof. In three cases only, is the improvement in the fit over the single component models statistically significant above the 99% confidence level, based on an F-test for two additional parameters (the abundance of the RS was fixed to the solar value as was for the single component models). The $`\chi ^2`$ for these sources are shown in bold. For these sources (src 5,6,29), the resulting photon index is much steeper than with the PO models alone. This is expected since now the soft excess if modeled by the additional thermal component. The temperatures of this thermal component are generally between 0.3 and 1.0 keV, with a few sources having temperatures up to 3.0 keV. Although with the current data we cannot discriminate between one and two component models for most sources, the PO+RS fits indicate that in 5 cases there may be a highly obscured ( $`\mathrm{N}_\mathrm{H}`$$`>0.3`$ cm-2, 10 times higher than the Galactic column) hard component.

In order to constrain any excess absorption for the hard component and assess its significance we fitted the spectra using the following method. First we fitted the spectra above 2.0 keV or 1.5 keV (depending on the number of bins in the spectrum) with a single unabsorbed power-law. Then we fitted the full band spectrum with a PO and RS model with the photon index fixed to the previously obtained value and the other parameters free to vary. The results from this fit for the sources for which it was possible (i.e. those with more than 3 bins above 1.5keV) are presented in Table 7: in Column (6) we give the photon index (determined from the fit above 2.0 keV), in Columns (7) and (8) we give the temperature of the RS component and the overall absorption (in units of $`10^{22}\mathrm{cm}^2`$) and in Column (9) we give the $`\chi ^2`$ and the number of dof for this fit. Then for the sources which gave in the previous fit a reduced $`\chi ^2`$ higher than 1.0, we included a second absorber applying it only to the PO component and we fitted the spectrum again. For this fit the photon index of the PO component was also fixed to the previous values. These results are also presented in Table 7. In Column (10) we give the the best fit temperature for the thermal component and in Columns (11) and (12) the intrinsic absorption of the hard component and the overal absorbing column density. Finally in Column (13) we give the $`\chi ^2`$ and the number of d.o.f for the fits after adding the second absorber. Although generally the best fit absorption of the PO component is higher than the overall absorption, for none of the sources the improvement of the fit is statistically significant. However, the PO model for some sources predicts unrealistically high luminosities ($`\mathrm{L}_\mathrm{X}>10^{42}`$ erg s-1), while the double component models predict more realistic luminosities. This suggests that the latter may be a more realistic model, although keeping in mind that given the limitations of the present data we have no proof of it.

Fits with more complex models such as power-law plus a single or multi-temperature black-body component, appropriate to model the emission of X-ray binaries (e.g. Nagase 1989), did not yield a statistically significant improvement in the fits (above the 99% confidence level based on an F-test).

#### 2.4.3 Line features

The residuals after the PO fits show evidence for emission lines in the spectra of sources 11, 16, 37, 42 and 45 (Fig. 5). We therefore modeled the spectrum of these sources with an absorbed power-law model to which we added narrow gaussians to account for the emission lines. The addition of these Gaussian components is significant at the 99% confidence level in three cases (sources 11, 16 and 42). The parameters of the lines for these three sources are presented in Table 8. In this Table Column (1) gives the source ID, Column (2) gives the energy of the line, Column (3) gives the identification of the line, Column (4) gives the normalization of the Gaussian and Column (5) gives the $`\chi ^2`$ and the number of degrees of freedom. Regarding the identification of the lines, we note that these strongly depend on the exact energy of the line. Therefore, although the emitting element is well identified, there is some uncertainty in the identification of the exact species which is reflected in the ranges given in column (3) of Table 8. Surprisingly, in every spectrum we find emission lines from different species suggesting that the gas parameters and the environment of the sources are very diverse. Inspection of high S/N ACIS-S3 background spectra (Markevitch et al. 2001)777http://hea-www.harvard.edu/ maxim/axaf/bg/index.html show that there are two very strong emission lines at $`1.7`$ keV and $`2.2`$ keV. Only the line we detect in source 16 at 2.17 keV could be due to residual background. All other lines are likely to originate from the source. In Fig. 6 we show the spectra of these three sources together with a PO model of the continuum radiation and the Gaussians used to model the emission lines. The bottom panel of the spectra shows the ratio of the data and the model.

#### 2.4.4 Hardness ratios

There are 19 sources with less than 50 net counts, for which it is not possible to perform any spectral fitting, even with the simplest models. For this reason we tried to determine their spectral parameters by using hardness ratios (HR). We calculated hardness ratios using the (0.3-1.0) keV (soft; S), (1.0-2.5) keV (medium; M) and (2.5-7.0) keV (hard; H) bands. We selected these bands after detailed modeling which showed that these are optimal to distinguish between different spectral models, when combined in HR-HR diagrams. We define three hardness ratios as: $`\mathrm{HR1}=(\mathrm{S}\mathrm{H})/(\mathrm{S}+\mathrm{H})`$, $`\mathrm{HR2}=(\mathrm{S}\mathrm{M})/(\mathrm{S}+\mathrm{M})`$ and $`\mathrm{HR3}=(\mathrm{M}\mathrm{H})/(\mathrm{M}+\mathrm{H})`$. The HRs were calculated by extracting the unbinned spectra and binning them in the S, M, and H energy bands. To subtract the background we followed the same procedure but for annular regions around the sources. Then the background was rescaled to the area of the source extraction cell and sutbracted from the source counts in each band.

To calibrate the HR-HR diagrams, we calculated hardness ratios from simulated absorbed power-law and thermal plasma spectra, for a wide range of parameters. We used the response matrix for an on-axis point source. Based on these model spectra we created grids on plots involving two different hardness ratios. We checked the accuracy of these grids by comparing them with the hardness ratios of sources with good spectral fit results (see Fig. 7). In this figure the red grids correspond to the PO models and the green grids correspond to the RS models. In the sides of the grids we give the values of the model parameters. The column density in both cases ranges between 0.01$`\times 10^{22}`$ cm-2 and 0.5$`\times 10^{22}`$ cm-2 ((0.01, 0.025, 0.05, 0.075, 0.10, 0.12, 0.25, 0.5)$`\times 10^{22}`$ cm-2). The temperature ranges between 2.0 and 10.0 keV (for lower temperatures these grids become degenerate and is not possible to distinguish between different temperature), and the photon index ranges between 0.5 and 2.0. In general the agreement between the spectral parameters derived from spectral fits and those estimated from the HR-HR diagrams is very good. Only weak sources with a soft excess appear to have much lower $`\mathrm{N}_\mathrm{H}`$ than derived by spectral fitting. This is because the determination of $`\mathrm{N}_\mathrm{H}`$ from hardness ratios depends on the S-M ratio. If a source has even a small soft excess, this results in a larger S-M ratio and therefore a smaller estimated $`\mathrm{N}_\mathrm{H}`$. The estimated spectral parameters are consistent between the different diagrams. The HR1-HR3 diagram is best suited for determining the spectral parameters, because of the larger separation of the hardness ratios for different spectral models. This is very important in the case of very faint sources where the error bars in the HRs are very large. Fig. 8 shows the hardness ratio diagrams for the faint sources. We were able to apply these diagrams to 16 sources which were detected in at least two of the three bands. This figure shows that the majority of the sources have soft spectra. Also some of them may have an additional soft thermal component.

### 2.5 Obscured sources

Spectral fits with a single power-law suggested that 12 sources are observed through a column density at least 10 times higher than that of the Galactic line of sight $`\mathrm{N}_\mathrm{H}`$ (Table 5). However, for two of them (src 10, 18) a thermal model gives a better fit with a lower column density (constistent with the Galactic $`\mathrm{N}_\mathrm{H}`$). For three other sources (src 5, 6, 29) a composite thermal/power-law model also gives a better fit with an absorption much lower than that required by the single power-law model. For one additional source (src 40) the composite model also gives a significantly lower $`\mathrm{N}_\mathrm{H}`$ although the improvement in the fit was not statistically significant. We consider this source as non obscured since its spectrum (Fig 5) does not show any significant cutoff below 1.0 keV. For the other sources for which it was possible to perform more complicated fits, the composite models indicate absorption for the hard component slightly lower than the overal absorption measured with the single PO model (Table 7), although the improvent in the fit over the latter is not statistically significant. Therefore, we conclude that we have 6 highly obscured sources (i.e. $`\mathrm{N}_\mathrm{H}`$$`>0.34\times 10^{22}`$ cm-2): sources 12, 24, 25, 34, 35 and 36.

### 2.6 Spectral trends

The spectra of X-ray sources are useful tools for understanding their nature, since different types of objects have different spectral signatures. Pulsar XRBs as well as BH binaries in low state tend to have hard spectra ($`\mathrm{\Gamma }<2.0`$) (e.g. Nagase 1989, Tanaka & Lewin 1995). The spectra of BH binaries in high state and unmagnetised neutron stars are dominated by a multi-temperature disk-BB spectrum (e.g. van Paradijs 1999, Tanaka & Lewin 1995). SNRs typically exhibit soft thermal emission (kT $`0.55.0`$ keV) (Schlegel 1995), although young supernovae detonating in dense environments (cSNRs) can have much harder spectra (Plewa 1995, Terlevich 1994).

From both the spectral fit results and the HR diagrams, it is clear that the discrete sources in the Antennae cover a very wide range of spectral parameters. Of course this is not surprising for such a large number of sources spanning 3 orders of magnitude in luminosity, and for an actively star-forming galaxy, which may contain several generations of stellar objects. We searched for possible trends in the spectral shape, in order to understand which types of objects dominate in each luminosity range. First we plotted the source luminosity (corrected for galactic absorption) against the best fit photon index from the single component model (Fig. 9). We plot the point like sources in red, and the extended sources in blue. The typical uncertainties in the luminosity are $`5\%`$ for the most luminous sources and $`20\%`$ for the faintest. We see a trend for more luminous sources to have harder spectra. When we split the sample in three subsamples containing $`10`$ sources each, we find that the mean photon index for the first bin ($`\mathrm{L}_\mathrm{X}<5\times 10^{38}`$ erg s-1) is $`6.89\pm 2.15`$, for the second ($`5\times 10^{38}<\mathrm{L}_\mathrm{X}<10^{39}`$ erg s-1) $`2.19\pm 0.62`$, and for the third ($`\mathrm{L}_\mathrm{X}>10^{39}`$ erg s-1) $`1.56\pm 0.27`$. These are represented by the black points with the error bars in the bottom of Fig. 9. These results suggest that there may be a marginal statistically significant trend for more luminous sources to have harder spectra ($`2\sigma `$ difference between the faintest and brightest sources).

A plot of the X-ray luminosity against $`\mathrm{N}_\mathrm{H}`$ (also in Fig. 9) shows that there is no correlation between the intrinsic X-ray luminosity and the amount of absorption.

## 3 Summary and Conclusions

In this paper we have described the detailed analysis of the spatial, temporal and spectral properties of the discrete X-ray sources detected with a 72 ks Chandra ACIS-S observation of the Antennae galaxies. Our results are summarized below:

1. We used two different algorithms for source detections: the sliding box celldetect and the wavelet transform wavdetect. We find that both algorithms lead to the detection of 49 X-ray sources. Three additional sources are found with celldetect, but visual inspections of the image shows that they are peaks of the diffuse emission of the Antennae. The measured net source counts from the two detection methods typically agree within 30%, the biggest discrepancies occurring in confuse fields, where image inspection shows that wavdetect performs better than celldetect. For this reason we use the wavdetect results to derive the source fluxes.
2. We detect a total of 49 sources within the optical area of the Antennae. The limiting luminosity is $`10^{38}`$ erg s-1 ($`\mathrm{H}_0=50`$ km s-1 Mpc-1; $`5\times 10^{37}`$ erg s-1, for $`\mathrm{H}_0=75`$ km s-1 Mpc-1). 31 sources (38 for $`\mathrm{H}_0=75`$ km s-1 Mpc-1) have luminosities below $`10^{39}`$ erg s-1 while 18 (11) have luminosities well in the ULX range, i.e. well in excess of the Eddington limit for a neutron star accretor. Based on the expected number of serendipitous sources in the field, only 2-3 sources may not be associated with the Antennae, and the probability of chance association is $`<1`$ for the ULXs (see Paper I).
3. We analyzed the spatial extent of the 20 sources detected with 100 counts or more, by comparing the spread of counts with the Chandra$`+`$ACIS-S PSF. All but 6 sources are point-like. Six sources have extended components, with extent ranging between 3” and 15” ($`4002100`$pc). These extended sources include the two nuclear regions.
4. Two luminous sources are found to vary within the observation with the KS test. By comparing the Chandra data with the ROSAT HRI observations of the Antennae (5” angular resolution) of Fabbiano et al (1997), we find evidence of long-term variability (year timescales) in 3 more sources, one of which includes a Chandra variable source.
5. We extracted ACIS-S spectra for all the 49 detected sources, and derived X-ray colors for all of them. Model fitting was performed on all the 31 sources with more than $`50`$ net counts and for which the data could be grouped in at least 5 spectral bins.
6. We fitted the data with both an absorbed power-law model and a thermal Raymond model with solar abundances and absorption. In both cases absorption columns in excess of the Galactic line of sight $`N_H`$ are favored. Two-component models (power-law plus Raymond) were fitted to the spectra of 10 sources for which a large enough number of counts was detected. With the exception of one source, no evidence was found in these composite models for an intrinsically absorbed power-law. Using more complex models, such as power-law plus multi-temperature black-body, which fit well X-ray binary spectra, did not produce any detectable improvement in the fit. Three sources show significant fit residuals suggesting the presence of emission lines. Six sources appear to have absorbed spectra, with $`\mathrm{N}_\mathrm{H}`$ higher than 10 times the Galactic along the line of sight.
7. We used Hardness Ratio (HR) diagrams, to derive spectral information from all the detected sources. We experimented with different selections of energy boundaries, and chose 0.3-1 keV (Soft), 1-2.5 keV (Medium), and 2.5-7.0 keV (Hard), because these offer the best differentiation of models in an HR-HR diagram. We calibrated these diagrams by simulating both power-law and thermal spectra for a variety of parameters, and compared these results with the results of the spectral fits, with good agreement. By using HR diagrams, we find that the majority of low-luminosity sources have soft spectra.
8. By plotting the power-law photon indices from the single power-law fits against source luminosity, we find an indication of an overall spectral trend, with the most luminous sources exhibiting very hard spectra ($`\mathrm{\Gamma }1.2`$), while softer emission is prevalent at the low luminosities (with values of $`\mathrm{\Gamma }>3`$ in some cases). No luminosity – $`N_H`$ trend is observed.

Our results suggest that a good fraction of the most luminous sources (ULXs) in the Antennae may be compact accretion binaries. This conclusion is supported by the presence of variability in few sources and by the luminosity – photon index trend suggesting that hard sources (consistent with X-ray binary emission) dominate the source population at high luminosities. At lower luminosities, instead, the softer X-ray spectra may suggest the presence of supernova remnant emission as well. We must however remember that these fainter spectra may include a sizeable amount of residual emission from the soft hot ISM in the spectral extraction region (see Paper I).

A detailed discussion of these results, augmented by the X-ray analysis of the average spectral properties in different source luminosity ranges, and supported by multi-wavelength comparison to constrain the nature of the emitting sources, is presented in Paper-III (Zezas et al. 2002).

We thank the CXC DS and SDS teams for their efforts in reducing the data and developing the software used for the reduction (SDP) and analysis (CIAO). We thank Martin Ward, Jeff McClintock, Andrea Prestwich and Phil Kaaret for useful discussions on these results. This work was supported by NASA contract NAS 8–39073 (CXC) and NAS8-38248 (HRC).


text:
# 1 Introduction

## 1 Introduction

Various massive fields like the gravitino, the Polonyi, the moduli, and the dilaton figure in supersymmetric and string theory models. The corresponding particles are long-lived. These fields pose quite serious cosmological problems. If they decay during baryogenesis or nucleosynthesis(BBN), the baryon to photon ratio may be greatly diluted in the first case and nuclear abundances distorted in the second. If they are stable, there is a problem of over-abundance. People have studied the problem in detail and proposed various ranges of masses and decay rates to minimise damage to standard cosmology.

An interesting proposal is that of thermal inflation. Quite apart from primordial inflation, a scalar field, called the flaton, is used to generate inflation at late times so that the potentially dangerous fields mentioned above may be diluted away. Typically, the inflation starts as the thermal energy density falls below the flaton potential at the origin at a temperature of about $`10^7`$ Gev, and stops when the flaton vev reaches $`10^{12}10^{14}`$ Gev at a temperature $``$ the flaton mass, which may be taken to be $`10^3`$ Gev. Such a flaton field will go on decaying into the Mev-scale temperature era, and the parameters in thermal inflation must be such that the entropy generated due to flaton decay may not affect nucleosynthesis.

Recently, the effect of the decay of a massive particle, like the flaton, on nucleosynthesis has been subjected to detailed investigation. As the parameters and decay schemes of such particles are yet to be generally agreed upon, it is useful to address the general problem of the effect of late-time entropy production, during reheating, in the context of BBN. In the present work, neutrino decoupling is studied in the presence of such late-time entropy production, neutrino decoupling temperature being one of the key parameters determining BB nucleosynthesis. Bernstein’s method of ”pseudo chemical potentials” is used to introduce the entropy production rate directly into the Boltzmann equation for the neutrino number density. The decoupling temperature is estimated in the presence of this entropy generation.

The paper is organised as follows. Section 1 is this Introduction. In Section 2, the late-time development of a universe, once dominated by a scalar field, is set out with some modification of the usual formalism. In Section 3, the method of ”pseudo chemical potentials” is applied to modify the Boltzmann equation for the neutrino number density and incorporate into it the late-time entropy production rate directly. In Section 4, the effect of late-time entropy production on neutrino decoupling is studied, using the results of the previous sections. In Section 5, conclusions are stated.

## 2 Reheating at Late Times

In this section, a two-component universe is considered, made up of the scalar field energy density $`\rho _\varphi `$ and a radiation energy density $`\rho _R`$. It is assumed , further, that either the scalar field parameters have not supported resonance, or, in any case, the times considered are late enough to be beyond a possible preheating  interlude.

Taking the scalar decay rate constant to be $`\mathrm{\Gamma }`$ , the equation for the evolution of the scalar energy density is

$$\frac{}{t}\rho _\varphi +3H\rho _\varphi =\mathrm{\Gamma }\rho _\varphi .$$
(1)
Defining $`\mathrm{\Phi }=a^3\rho _\varphi `$ and $`R=a^4\rho _\varphi `$, where $`a`$ is the scale factor, (1) becomes

$$\dot{\mathrm{\Phi }}=\mathrm{\Gamma }\mathrm{\Phi }.$$
(2)
From

$$\frac{}{t}[a^3(\rho _\varphi +\rho _R)]+p_R\frac{}{t}a^3=0,$$

one obtains

$$\dot{R}=a\mathrm{\Gamma }\mathrm{\Phi }.$$
(3)
The Friedmann equation is

$$H^2=\frac{K^2}{a^4}(\mathrm{\Phi }a+R),\text{ with }K^2=\frac{8\pi }{3M_{Pl}^2}.$$
(4)
### 2.1 Incomplete $`\mathrm{\Phi }`$-domination

In this subsection, it is assumed that $`R\mathrm{\Phi }a`$, so that the square and higher powers of $`R/(\mathrm{\Phi }a)`$ can still be neglected at these times. (4) gives

$$\dot{H}=\frac{3}{2}H^2(1+\frac{R}{3(\mathrm{\Phi }a+R)}).$$

With the $`R\mathrm{\Phi }a`$ asssumption, one can write

$$\dot{H}\frac{3}{2}H^2(1+\frac{R}{3\mathrm{\Phi }a}).$$
(5)
For full matter domination, at yet earlier times, $`R`$ is neglected compared to $`\mathrm{\Phi }a`$, and

$$\dot{H}\frac{3}{2}H^2,$$

leading to $`t=\frac{2}{3H}`$. In the era of incomplete $`\mathrm{\Phi }`$-domination, the correction term on the RHS of (5) is not neglected, and may be evaluated to the approximation

$$RR_I=\frac{dR}{da}(aa_I),$$
(6)
where $`t_I`$ refers to some initial epoch such that $`aa_I`$. Also, if it is supposed that the scalar decay produces sufficiently copious radiation, $`R_IR`$. As a correction term is being dealt with, these approximations should not cause much deviation from the actual evolution. Then, for $`t`$ sufficiently later than $`t_I`$, but within the regime under consideration, one may write, in the correction term on the RHS of (5),

$$R\frac{\mathrm{\Gamma }}{H}\mathrm{\Phi }a,$$
(7)
using (3) and (6). (5) becomes

$$\dot{H}=\frac{3}{2}H^2(1+\frac{\mathrm{\Gamma }}{3H}).$$
(8)
Now, evolution will be described by the new variable

$$x=\frac{\mathrm{\Gamma }}{H},$$
(9)
instead of the time $`t`$. (4) can be written as

$`H^2`$ $`=`$ $`{\displaystyle \frac{K^2\mathrm{\Phi }}{a^3}}(1+{\displaystyle \frac{R}{\mathrm{\Phi }a}})`$ (10)
$`=`$ $`{\displaystyle \frac{K^2\mathrm{\Phi }}{a^3}}(1+x),`$ (11)
using (7) in the correction term on the RHS of (10).

Then, (2) leads to the evolution equation for $`\mathrm{\Phi }`$

$$\frac{d\mathrm{\Phi }}{dx}=\frac{2}{3+x}\mathrm{\Phi },$$
(12)
with the solution

$$\mathrm{\Phi }=\mathrm{\Phi }_I\frac{(1+\frac{1}{3}x_I)^2}{(1+\frac{1}{3}x)^2},$$
(13)
in the approximation $`x1`$. To obtain $`R`$, (13) is substituted in (3), $`a`$ being calculated from (11), using (9). The result is

$$\frac{dR}{dx}=\frac{2}{3}(\frac{K}{\mathrm{\Gamma }})^{\frac{2}{3}}x^{\frac{2}{3}}\mathrm{\Phi }_I^{\frac{4}{3}}\frac{(1+\frac{1}{3}x_I)^{\frac{8}{3}}(1+x)^{\frac{1}{3}}}{(1+\frac{1}{3}x)^{\frac{11}{3}}}.$$
(14)
Integration will give $`R`$.

To compare with earlier work, at early times, one goes to the approximation $`x,x_I1`$. (14) becomes

$$\frac{dR}{dx}=\frac{2}{3}(\frac{K}{\mathrm{\Gamma }})^{\frac{2}{3}}\mathrm{\Phi }_I^{\frac{4}{3}}x^{\frac{2}{3}}.$$

On integration, taking $`R_IR`$,

$`R`$ $`=`$ $`{\displaystyle \frac{2}{5}}({\displaystyle \frac{K}{\mathrm{\Gamma }}})^{\frac{2}{3}}\mathrm{\Phi }_I^{\frac{4}{3}}(x^{\frac{5}{3}}x_I^{\frac{5}{3}})`$
$`=`$ $`{\displaystyle \frac{2}{5}}{\displaystyle \frac{\mathrm{\Gamma }}{K}}\mathrm{\Phi }_I^{\frac{1}{2}}(a^{\frac{5}{2}}a_I^{\frac{5}{2}}),`$
which agrees with .

### 2.2 Incomplete Radiation domination

Now, one may consider even later times when $`\mathrm{\Phi }aR`$, such that $`\mathrm{\Phi }a/R`$ cannot be neglected, but its higher powers can. The Friedmann equation reads

$$H^2=\frac{K^2}{a^4}R(1+\frac{\mathrm{\Phi }a}{R}),$$
(15)
which leads to

$$\dot{H}=2H^2[1\frac{\mathrm{\Phi }a}{4(\mathrm{\Phi }a+R)}].$$
(16)
If, well into this epoch, the correction term on the RHS of (16) is neglected, the full radiation domination relations are found:

$`H`$ $`=`$ $`{\displaystyle \frac{1}{2t}},\text{ and}`$
$`a`$ $`=`$ $`At^{\frac{1}{2}},`$ (17)
$`A`$ being a constant.

(2) has, as solution, a falling exponential in $`t`$, viz. $`\mathrm{\Phi }e^{\mathrm{\Gamma }t}`$. The new evolution variable $`x=\mathrm{\Gamma }/H`$ is now sought to be introduced in place of $`t`$. To do this, instead of taking the falling exponential in $`t`$ directly, suitable approximations to the correction terms on the RHS of (15) and (16), are first worked out. Let $`t_0`$ be a sufficiently late epoch, when $`\mathrm{\Phi }=\mathrm{\Phi }_00`$. Then, for use only in the correction terms, one takes

$`\mathrm{\Phi }\mathrm{\Phi }_0=\stackrel{~}{\mathrm{\Phi }}({\displaystyle \frac{1}{t}})\stackrel{~}{\mathrm{\Phi }}({\displaystyle \frac{1}{t_0}})={\displaystyle \frac{d\stackrel{~}{\mathrm{\Phi }}}{d\frac{1}{t}}}|_{t_0}({\displaystyle \frac{1}{t}}{\displaystyle \frac{1}{t_0}}).`$
Neglecting $`\mathrm{\Phi }_0,1/t_0`$ compared to $`\mathrm{\Phi },1/t`$, respectively, an approximation

$$\mathrm{\Phi }\frac{B}{t},$$
(18)
will be used only in the correction terms, i.e. in the correction terms, the falling exponential will be approximated by a rectangular hyperbola. B is a constant. Using this approximation in the correction terms, the final form of $`\mathrm{\Phi }`$ will be found below to be a falling exponential in $`(1/2)x`$, with a pre-exponential correction.

Next, a similar approximation is considered for $`R`$. It ought to be mentioned that $`R`$ refers to the total radiation present, and not only to that produced by decay. However, the change in $`R`$ is due to $`\varphi `$ decay and consequent entropy production . In the absence of this decay, $`\dot{R}=0`$.

Using (17) and (18) in (3), and, integrating, one obtains, for use only in the correction terms,

$$RR_E2AB\mathrm{\Gamma }(t^{\frac{1}{2}}t_E^{\frac{1}{2}}).$$

If $`t_E`$ is sufficiently early compared to $`t`$, though within the regime under consideration, and there is sufficiently copious radiation production since $`t_E`$, it is sufficient to take

$$R2AB\mathrm{\Gamma }t^{\frac{1}{2}}$$

in the correction terms. (17) and (18) are now used to give, in the correction terms,

$`{\displaystyle \frac{R}{\mathrm{\Phi }a}}`$ $``$ $`{\displaystyle \frac{\mathrm{\Gamma }}{H}},\text{ once again, as in (}\text{7}\text{),}`$ (19)
$`=`$ $`x.`$
The approximation is now $`x1.`$ (15) and (16) become

$`H^2`$ $`=`$ $`{\displaystyle \frac{K^2}{a^4}}R(1+{\displaystyle \frac{1}{x}}),`$ (20)
$`\text{and, }\dot{H}`$ $`=`$ $`2H^2{\displaystyle \frac{x+\frac{3}{4}}{x+1}}.`$
This last equation and (2) give

$$\mathrm{\Phi }=\mathrm{\Phi }_E(\frac{4x_E+3}{4x+3})^{\frac{1}{8}}e^{\frac{1}{2}(xx_E)}$$
(21)
One can see that this solution is dominated by the falling exponential. The pre-exponential factor arises because radiation domination is not yet complete, and $`(1/2)x`$ is not quite $`\mathrm{\Gamma }t.`$ For sufficiently late times $`t,t_E`$, one has $`x,x_E1`$, (20) gives $`H=1/(2t)`$, in the usual way, and (21) reduces to

$$\mathrm{\Phi }=\mathrm{\Phi }_Ee^{\mathrm{\Gamma }(tt_E)},$$

which agrees with .

From (20),

$$a(\frac{K}{H})^{\frac{1}{2}}R^{\frac{1}{4}}(1+\frac{1}{x})^{\frac{1}{4}}.$$
(22)
Now, on integrating (3), using (21) and (22), one has

$`R^{\frac{3}{4}}`$ $`=`$ $`R_{\mathrm{}}^{\frac{3}{4}}`$ (23)
$`{\displaystyle \frac{3}{8}}({\displaystyle \frac{K}{\mathrm{\Gamma }}})^{\frac{1}{2}}(x_E+{\displaystyle \frac{3}{4}})^{\frac{1}{8}}\mathrm{\Phi }_Ee^{\frac{1}{2}x_E}{\displaystyle _x^{\mathrm{}}}x^{\frac{3}{8}}e^{\frac{1}{2}x}𝑑x,`$
for $`x1`$, where $`R_{\mathrm{}}`$ is the value of $`R`$ after decay is finished. The integral in (23) is an incomplete Gamma function, and can be numerically evaluated for any value of $`x`$.

The temperature is defined from

$$R=g^{}\frac{\pi ^2}{30}a^4T^4,$$
(24)
where $`g^{}`$ is the effective number of relativistic degrees of freedom.

The entropy production originates from

$$dS_\varphi =\frac{1}{T}d\mathrm{\Phi }.$$
(25)
## 3 Entropy Production and the Boltzmann Equation for the Number Density with Scalar Decay

The effect of entropy production on decoupling is studied here with reference to the particular problem of the electron neutrino, $`\nu ,`$ decoupling from a thermal bath of electrons, positrons, and photons, in the presence of a decaying scalar field. It is assumed that the scalar field has a very small branching ratio into neutrinos. So, the main process which changes the number density $`n`$ of the $`\nu `$ neutrinos is

$$\nu (k)+\overline{\nu }(\overline{k})(p)+\overline{}(\overline{p}),$$
(26)
where $`,\overline{}`$ are fermions, and, $`k`$ is the energy-momentum 4-vector $`(E_k,\stackrel{}{k})`$. Apart from $`e^{},e^+`$, the relevant fermions, at the epoch of $`\nu `$ decoupling, are the neutrinos of the other families.

If there is no other process like scalar decay, it is usual to assume  that, in the absence of Fermi degeneracy, the distribution functions may be written as

$`f_\nu (k)=f^{}(k)`$ $`=`$ $`e^{\alpha ^{}(t)\beta E_k}`$
$`f_e(p)`$ $`=`$ $`e^{\beta E_p}`$
$`f_{\overline{\nu }}(\overline{k})=\overline{f^{}}(\overline{k})`$ $`=`$ $`e^{\alpha ^{}(t)\beta E_{\overline{k}}}`$
$`f_{e^+}(\overline{p})`$ $`=`$ $`e^{\beta E_{\overline{p}}},`$ (27)
where $`T=\frac{1}{\beta }`$ is the temperature of the thermal bath, and $`\alpha ^{}`$ is a time-dependent ”pseudo chemical potential”, introduced to take into account the departure of the decoupling neutrinos from equilibrium. (In (27), it has been assumed that $``$ is the electron.)

Decoupling is, in this case, governed by the integrated Boltzmann equation

$$\dot{n}+3Hn=<\sigma |v|>(n^2n_{EQ}^2),$$
(28)
where $`n_{EQ}`$ is the equilibrium number density,

$$n_{EQ}=\frac{gd^3k}{(2\pi )^3}e^{\beta E_k}=T^3/\pi ^2,$$
(29)
$$nn_{EQ}e^\alpha ^{}\text{[15]},$$
(30)
and

$$<\sigma |v|>=\frac{1}{n_{EQ}^2}𝑑K𝑑\overline{K}e^{\beta E_k}e^{\beta E_{\overline{k}}}I,$$
(31)
$`I`$ being the invariant integral

$$I=𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})||^2.$$
(32)
$`g`$ is the number of spin degrees of freedom, $`dK=gd^3k/[(2\pi )^32E_k]`$ etc., and $`||^2`$ the spin-averaged square of the modulus of the relevant matrix element.

Considering the process (26), and, assuming CP-invariance and the absence of Fermi degeneracy, the Boltzmann equation for the neutrino number density $`n`$ may be written

$$\frac{1}{a^3}\frac{}{t}(a^3n)=\underset{}{}𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})|_{}|^2(f^{}\overline{f}^{}f\overline{f}).$$
(33)
$`f^{}`$ corresponds to the $`\nu `$ neutrino and $`f`$ to the fermion $``$.
Then,

$`{\displaystyle \frac{1}{a^3}}{\displaystyle \frac{}{t}}(a^3n)`$ $`=`$ $`{\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})||^2(f^{}\overline{f}^{}f\overline{f})}`$ (34)
$``$ $`{\displaystyle \underset{i}{}}{\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})|_i|^2(f^{}\overline{f}^{}f_i\overline{f}_i)},`$
where, on the RHS, the first term refers to the process

$$\nu +\overline{\nu }e^{}+e^+,$$
(35)
and the second to the processes

$$\nu +\overline{\nu }\nu _i+\overline{\nu }_i,i=\mu ,\tau ,$$
(36)
$`f_i`$ being the $`\nu _i`$ distribution function.

The injection of entropy will cause changes in the distribution functions. One way to tackle the situation is to use the Boltzmann equation for the $`\nu `$ distribution function, instead of the integrated Bolzmann equation for the number density, include $`\nu +e^{}(e^+)\nu +e^{}(e^+)`$ processes, and proceed numerically . However, the integrated Boltzmann equation for the number density will be used here to study decoupling at late times, when the entropy generation rate is small, in the following way.

It has been assumed, in accordance with current models, that the $`\varphi `$ does not decay into neutrinos. But, ref. finds that, for reheat temperatures $`<7Mev`$, the neutrino distribution function is materially affected by the decay. How can entropy injected into the $`e,e^+,\gamma `$ sector affect neutrino distribution function and, hence, neutrino density? The transfer of entropy must occur via process (35), and this must lead to extra terms on the RHS of (28) if $`n`$ is to be affected. From (33), this can only occur if $`f`$, the electron distribution function, changes due to entropy injection.

The electromagnetic interactions, faster than the expansion, will tend quickly to thermalise the $`e,e^+,\gamma `$ (i.e. establish kinetic and chemical equilibrium in this sector). But recent work has shown that this thermalisation of $`\varphi `$ decay products is not nearly instantaneous. A massive $`\varphi `$, as it goes on decaying, continually injects highly energetic light particles, with energies much above the thermal mean. There is no possibility of inverse decay and so the term from scalar decay in the Botzmann equation for, say, the electron distribution function opposes the equilibriating effect of the electromagnetic interaction. As long as the decay does not become negligible, complete thermalisation of the distribution functions of the decay products is not assured. How quickly complete thermalisation will occur, i.e. the epoch of thermalisation, depends, not only on the relative rapidity of the electromagnetic interactions and the expansion, but also on the parameters of $`\varphi `$. In fact, the requirement of thermalisation before a specific epoch (e.g. BBN) has been used to derive bounds on these parameters.

The thermalisation of rapidly interacting, light decay products, of a heavy boson was considered in ref.. There, the progress, in time, of the light particle distribution function $`g`$ was directly considered. A solution $`g=e^{\alpha \beta E}`$ of the Boltzmann equation was mooted, where $`\alpha `$ was a time-dependent ”pseudo chemical potential”, although, the actual calculation was done with a $`\delta `$ function distribution, because the problem was simplified by assuming that all the heavy bosons decayed at one instant. Following this lead, it will be supposed, in this paper, that the injection of entropy into the bath, due to scalar decay, can be taken care of by introducing a small time-dependent ”pseudo-chemical potential” $`\alpha (t)`$ into the electron and positron distribution functions, in addition to the potential $`\alpha ^{}(t)`$ already introduced in the $`\nu ,\overline{\nu }`$ distribution functions.

Of course, the electromagnetic interactions, being so fast, will cause $`\alpha `$ to be very small.

In effect, then, the scheme in (27) will be changed, in the absence of Fermi degeneracy, to

$`f_\nu (k)=f^{}(k)`$ $`=`$ $`e^{\alpha ^{}\beta E_k}`$
$`f_e(p)=f(p)`$ $`=`$ $`e^{\alpha \beta E_p}`$
$`\overline{f^{}}`$ $`=`$ $`e^{\alpha ^{}\beta E_{\overline{k}}}`$
$`\overline{f}`$ $`=`$ $`e^{\alpha \beta E_{\overline{p}}}.`$ (37)
The entropy density $`s_j`$, contributed by a particle $`j`$ with a distribution function $`h_j(q)`$, may be defined, in the absence of degeneracy, as

$$s_j=g_j\frac{d^3q}{(2\pi )^3}(h_jlnh_jh_j).$$

Following , the covariant divergence of the entropy density current for a component $`j`$, in a process like that described by (35), can be shown to be

$$\frac{1}{a^3}\frac{}{t}(s_ja^3)g_j\frac{d^3q}{(2\pi )^3}\frac{C_j}{E_j}lnh_j,$$

where $`C_j`$ is the relevant collision integral as defined $`^{b,}`$ 333This differs from the definition in  by a factor of $`E_j`$. in . For example, for the electron,

$`{\displaystyle \frac{1}{a^3}}{\displaystyle \frac{}{t}}(s_ea^3)`$ $`=`$ $`{\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})||^2}`$
$`\times `$ $`lnf(f^{}\overline{f^{}}f\overline{f}).`$
It can immediately be seen that the covariant divergence of the total entropy density current for the process (35) is, assuming energy conservation,

$`{\displaystyle \frac{1}{a^3}}{\displaystyle \frac{}{t}}(s_{tot}a^3)`$ $`=`$ $`{\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})||^2}`$ (38)
$`\times 2(\alpha ^{}\alpha )(e^{2\alpha }e^{2\alpha ^{}})e^{\beta (E_k+E_{\overline{k}})}`$
$``$ $`{\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})||^2}`$
$`\times 2\alpha ^{}(1e^{2\alpha ^{}})e^{\beta (E_k+E_{\overline{k}})}`$
$`+{\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})||^2}`$
$`\times (2\alpha )(1e^{2\alpha ^{}}+2\alpha ^{})e^{\beta (E_k+E_{\overline{k}})}.`$
Here, $`\alpha `$, being the effect of entropy injection, is assumed to be small at late times. Squares and higher powers of $`\alpha `$ have been neglected.

The first term on the RHS of (38) is the covariant divergence of the entropy density current for $`\alpha =0`$, and the second term is proportional to $`\alpha `$. So, it can be concluded that the second term measures, approximately, to first order in $`\alpha ,`$ the contribution of $`\varphi `$decay to the covariant divergence of the entropy density current, corresponding to process (35), viz.

$`ϵ{\displaystyle \frac{1}{a^3}}{\displaystyle \frac{S_\varphi }{t}}`$ $`=`$ $`2\alpha (1e^{2\alpha ^{}}+2\alpha ^{}){\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})||^2}`$ (39)
$`\times e^{\beta (E_k+E_{\overline{k}})},`$
where $`S_\varphi `$ is defined in (25) and $`ϵ`$ is the fraction of the entropy current from $`\varphi `$-decay which contributes to process(35).

Again, in the presence of $`\varphi `$decay,

$`{\displaystyle \frac{1}{a^3}}{\displaystyle \frac{}{t}}(a^3n)`$ $`=`$ $`{\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})||^2}`$
$`\times (e^{2\alpha ^{}}e^{2\alpha })e^{\beta (E_k+E_{\overline{k}})}`$
$``$ $`{\displaystyle \underset{i}{}}{\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})|_i|^2(f^{}\overline{f}^{}f_i\overline{f}_i)}`$
$`=`$ $`(n^2n_{EQ}^2)<\sigma |v|>`$
$``$ $`2\alpha {\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})||^2e^{\beta (E_k+E_{\overline{k}})}}`$
$``$ $`{\displaystyle \underset{i}{}}{\displaystyle 𝑑K𝑑\overline{K}𝑑P𝑑\overline{P}(2\pi )^4\delta ^4(p+\overline{p}k\overline{k})|_i|^2(f^{}\overline{f}^{}f_i\overline{f}_i)}`$ (41)
to first order in $`\alpha `$, using (30), (31), (37) and (34). $`<\sigma |v|>`$ refers to (35).

Now, the finer point that, in this epoch, the electron neutrino has both charged and neutral current interactions, while the other neutrinos have only neutral current interactions (muons and taons having already decoupled), will be neglected, and it will be assumed that the charged and neutral current interactions fall out of equilibrium together. Also, neutrinos of all types will be assumed to be massless. Then, apart from their contributions to $`\rho _R`$ as separate species, the neutrinos of the different types may be treated as identical, and it may be supposed that they decouple together. In this approximation, the difference between $`f^{}`$ and $`f_i`$ may be neglected, and the terms within the summation sign dropped, in comparison to the other terms on the RHS of (41).

So, to first order in $`\alpha ,`$ the integrated Boltzmann equation can be written, in the presence of $`\varphi `$decay, as

$`{\displaystyle \frac{1}{a^3}}{\displaystyle \frac{}{t}}(a^3n)`$ $`=`$ $`\dot{n}+3Hn`$ (42)
$`=`$ $`(n^2n_{EQ}^2)<\sigma |v|>+{\displaystyle \frac{ϵ}{(1e^{2\alpha ^{}}+2\alpha ^{})}}{\displaystyle \frac{1}{a^3}}{\displaystyle \frac{S_\varphi }{t}},`$
from (39) and (41). (42) will replace (28) in the presence of scalar decay.

## 4 Effect of Scalar Decay on Neutrino Decoupling

Decoupling will be supposed to set in when the RHS of (42) becomes less than $`3Hn`$. To analyse the two terms on the RHS, radiation domination is first assumed. A change of variables is introduced:

$$x=\frac{\mathrm{\Gamma }}{H},Y=nx^{\frac{3}{2}},Y_{EQ}=n_{EQ}x^{\frac{3}{2}}.$$
(43)
Also, (20) gives, in the era of neutrino decoupling,

$$H=\frac{4.461\times 10^{19}}{Gev}T^2(1+\frac{1}{x})^{\frac{1}{2}},$$
(44)
taking $`g^{}=10.75`$. The LHS of (42), then, becomes $`(2\mathrm{\Gamma }/x^{\frac{3}{2}})dY/dx`$. Near the start of the decoupling regime, one can put

$$Y=Y_{EQ}+\mathrm{\Delta },$$

where $`\mathrm{\Delta }`$ is very small, and it is possible  to set

$`{\displaystyle \frac{d\mathrm{\Delta }}{dx}}`$ $`=`$ $`0,\text{ and,}`$
$`{\displaystyle \frac{dY}{dx}}`$ $`=`$ $`{\displaystyle \frac{dY_{EQ}}{dx}}.`$ (45)
(43), (44) and (29) are used to find that $`dY_{EQ}/dx>0,`$ i.e. the LHS of (42) is positive. Use of (30) shows that each term on the RHS of (42) is of the same sign. So, for the equation to hold,

$$\alpha ^{}>0,$$
(46)
and each of the terms on the RHS of (42) is positive. This means, that when decoupling sets in, each term must be separately less than $`3Hn`$.

The value of $`\alpha ^{}`$ chosen will determine the temperature of decoupling. This choice is fixed as follows. In the absence of $`\varphi `$decay, there is only the first term on the RHS of (42). By setting this term equal to $`3Hn`$, a set of values for $`\alpha ^{}`$ is obtained for neutrino decoupling temperatures $`T_D`$ between $`0.75`$ and $`3`$ Mev, a range which more or less safeguards BBN. It may be mentioned that, at present, observed $`{}_{}{}^{4}He`$ abundance values, $`Y_p,`$ fall between $`0.234`$  and $`0.244`$ . A rough estimate, as indicated in , gives corresponding decoupling temperatures $`1`$ Mev.

Next, for these values of $`\alpha ^{}`$, the full RHS of (42) is put equal to $`3Hn`$, and the decoupling temperatures, corresponding to different values of the scalar decay constant $`\mathrm{\Gamma }`$, are worked out.

The first step uses

$$\frac{(n^2n_{EQ}^2)<\sigma |v|>}{3Hn}1.$$
(47)
As already discussed, in the evaluation of $`<\sigma |v|>`$, it is enough to consider the process (35), in the s-channel with $`Z`$ exchange and the t-channel with $`W`$ exchange. Evaluating (32), and, then, (31), using (29),

$`<\sigma |v|>`$ $`=`$ $`{\displaystyle \frac{8}{\pi }}G_F^2[(C_{Ve}+1)^2+(C_{Ae}+1)^2]T^2`$ (48)
$`=`$ $`{\displaystyle \frac{4.112\times 10^{10}}{Gev^4}}T^2.`$
The mass of the electron has been neglected. $`G_F`$ is the Fermi constant.

Taking $`H=\frac{4.461\times 10^{19}}{Gev}T^2,`$ in the absence of $`\varphi `$-decay, and using (29), (30), and (48), one obtains the values of $`\alpha ^{}`$ shown in Table I.

Now, the second term on the RHS of (42), i.e. the contribution of scalar decay, is considered. Using (21), this term becomes

$`{\displaystyle \frac{ϵ}{(1e^{2\alpha ^{}}+2\alpha ^{})}}{\displaystyle \frac{1}{a^3}}{\displaystyle \frac{S_\varphi }{t}}`$ $`=`$ $`{\displaystyle \frac{ϵ}{(1e^{2\alpha ^{}}+2\alpha ^{})}}{\displaystyle \frac{\dot{\mathrm{\Phi }}}{Ta^3}}`$ (49)
$`=`$ $`{\displaystyle \frac{ϵ}{(1e^{2\alpha ^{}}+2\alpha ^{})}}{\displaystyle \frac{\mathrm{\Gamma }\mathrm{\Phi }_E(\frac{4x_E+3}{4x+3})^{\frac{1}{8}}e^{\frac{1}{2}(xx_E)}}{Ta^3}}.`$
This term, for $`x>>1`$, will be dominated by the exponential. So, rather drastic assumptions will be made to estimate the pre-exponential, in the absence of phenomenological information about $`x_E\text{ and }\mathrm{\Phi }_E`$. $`x_E`$ refers to any epoch sufficiently early in the regime of incomplete radiation domination, and is set equal to 1, so that $`\mathrm{\Phi }_E=R_E/a_E`$. This is quite an approximation, because (21) is faithful to $`\mathrm{\Phi }e^{\mathrm{\Gamma }t}`$ provided $`x_E,x1`$. With this approximation, $`T_E`$ is, then, estimated by setting $`x=1`$ in (44). Next, one puts $`R_E=(\pi ^2/30)g_E^{}a_E^4T_E^4`$ and $`g_{SE}^{}a_E^3T_E^3=g_S^{}a^3T^3`$. This last is the weakest assumption because it is only true if there is conservation of entropy, with $`aT^1`$ behaviour. But, entropy is being produced, and, in fact, in the regime of $`\mathrm{\Phi }`$ domination, $`Ta^{\frac{3}{8}}`$ . In the early part of the regime of incomplete radiation domination, too, one expects a deviation from $`aT^1`$. However, the expectation is that the exponential dominates the entropy current and its effect on $`\nu `$ decoupling, and one can take the view that the pre-exponential is being estimated only to order of magnitude. Neglect of the deviation from $`aT^1`$ behaviour in the pre-exponential, and the other approximations, should not make too much difference to an order of magnitude calculation of lower bounds on the scalar decay constant.

With all these assumptions, the condition of neutrino decoupling in the presence of scalar decay becomes, using (44),

$`{\displaystyle \frac{(n^2n_{EQ}^2)<\sigma |v|>}{3Hn}}+{\displaystyle \frac{ϵ}{(1e^{2\alpha ^{}}+2\alpha ^{})}}{\displaystyle \frac{1}{a^3}}{\displaystyle \frac{\dot{S_\varphi }}{3Hn}}`$ $`=`$
$`{\displaystyle \frac{208.945sinh\alpha ^{}\mathrm{\Gamma }_0^{\frac{3}{2}}}{x^{\frac{1}{4}}(x+1)^{\frac{5}{4}}}}+{\displaystyle \frac{ϵ10.4927x^{\frac{5}{4}}(1+x)^{\frac{1}{4}}e^{\frac{1}{2}(x1)}}{(x+\frac{3}{4})^{\frac{1}{8}}(1e^{2\alpha ^{}}+2\alpha ^{})e^\alpha ^{}}}`$ $``$ $`1,`$ (50)
where

$$\mathrm{\Gamma }_0=\frac{\mathrm{\Gamma }}{10^{22}Gev}.$$

It must be remembered that $`x>1`$ in (50), corresponding to radiation domination.

### 4.1 Numerical Results

Three values of $`\alpha ^{}`$ are used from Table I, namely, those corresponding to decoupling in the absence of scalar decay at $`T_D=1,2,3Mev`$. In each case, the RHS of (50) is put equal to 1, and the resulting equation solved for $`x=x_D`$, taking different values of $`\mathrm{\Gamma }_0`$. The corresponding decoupling temperatures $`T_D^{}`$ are found from (44), using $`x=\mathrm{\Gamma }/H`$ to rewrite it as

$$\frac{T}{Mev}=\frac{14.9718\sqrt{\mathrm{\Gamma }_0}}{(x^2+x)^{\frac{1}{4}}}.$$
(51)
$`ϵ`$ would be $`1`$ if the entire entropy current of decay contributed to process (35). There are two points here. The decay modes of $`\varphi `$ are unknown, and it is assumed that decay products, more massive than the electron, rapidly decouple, and almost all the entropy current of decay ends up first in the $`e,e^+,\gamma `$ sector and, then, goes on to contribute to the entropy current of process (35) via a non-zero $`\alpha `$. Also, $`\alpha 0`$ means that there will be a small contribution of the entropy current from $`\varphi `$decay to the $`e,e^+,\gamma `$ processes also. But, because $`\alpha `$ is so small, $`ϵ`$ has been put equal to $`1`$.

The results are shown in Table II.

The results indicate that, for values of the scalar decay constant $`\mathrm{\Gamma }>10^{22}Gev`$, the neutrino decoupling temperature is not appreciably affected by scalar decay.

It is necessary to check what happens if there is matter domination. Taking $`x_I<<1`$ in (13), extrapolating this equation to $`\mathrm{\Phi }=\mathrm{\Phi }_E`$ at $`x=x_E=1`$, and evaluating $`\mathrm{\Phi }_E`$ with the same set of assumptions as in the case of radiation domination, the criterion of decoupling will become, in place of (50),

$`{\displaystyle \frac{(n^2n_{EQ}^2)<\sigma |v|>}{3Hn}}+{\displaystyle \frac{ϵ}{(1e^{2\alpha ^{}}+2\alpha ^{})}}{\displaystyle \frac{1}{a^3}}{\displaystyle \frac{\dot{S_\varphi }}{3Hn}}`$ $`=`$
$`{\displaystyle \frac{208.945sinh\alpha ^{}\mathrm{\Gamma }_0^{\frac{3}{2}}}{x^{\frac{1}{4}}(x+1)^{\frac{5}{4}}}}+{\displaystyle \frac{17.3934x^{5/4}(1+x)^{1/4}}{(1e^{2\alpha ^{}}+2\alpha ^{})e^\alpha ^{}(1+\frac{1}{3}x)^2}}1.`$ (52)
Matter domination implies that $`x<1`$ in (52). Writing the two terms on the RHS of (52) as A and B, $`A+B1`$. An examination of B will show that as the temporal variable x increases, B increases. So, if B predominates in an epoch, there can be no decoupling in that epoch. A little numerical work shows that for those values of $`\mathrm{\Gamma }`$ (greater than a critical value which depends on $`\alpha ^{}`$) for which A predominates over B, A+B does not fall below 1, for values of $`x<1`$ (matter domination ). So, no decoupling is possible in the epoch of matter domination.

The numerical work shows, therefore, that for neutrino decoupling to proceed without significant change, in the presence of scalar decay, there must be radiation domination, in the sense that $`x=\rho _R/\rho _\varphi >1`$, and the scalar decay constant $`\mathrm{\Gamma }\text{ must be }>10^{22}Gev`$.

This corresponds to reheating temperatures $`>8.5Mev`$, taking for the reheating temperature, the definition of ref. :

$$\frac{T_R}{Gev}=0.554(\frac{\mathrm{\Gamma }}{Gev})^{\frac{1}{2}}\sqrt{\frac{2.4\times 10^{18}}{Gev}}.$$

In , it was found that the neutrino distribution function was distorted, and the effective number of neutrino types, $`N_{eff},`$ started to decrease below three, as the reheating temperature fell below $`7`$ Mev, and standard BBN was endangered. So, the present results show broad agreement with those of  on the question of the minimum reheating temperature which safeguards standard BBN.

## 5 Conclusions

The effect of entropy injection due to scalar decay on neutrino decoupling has been studied here by introducing the entropy production rate term directly into the Boltzmann equation for the neutrino number density. The method adopted was to introduce a small ”pseudo chemical potential” $`\alpha `$ into the electron distribution function, in addition to the standard introduction of such a potential in the neutrino distribution function. As the electromagnetic interactions tend to thermalise the electron distribution function, this $`\alpha `$ is bound to be very small. Its function is to transmit the entropy current, arising from scalar decay, from the $`e,e^+,\gamma `$ sector, to the neutrino sector through processes like $`e^{}+e^+\nu +\overline{\nu }`$.

The conclusion drawn regarding the condition of validity of standard BBN from this study of neutrino decoupling, viz. reheating temperature$`>8.5`$ Mev, agrees broadly with the conclusions drawn in the literature from calculation of the form of the neutrino distribution function and the effective number of neutrino types, by numerical integration of the Boltzmann equation for the neutrino distribution function.

It is found that, in the presence of entropy injection due to scalar decay, the universe must not be matter dominated($`\rho _\varphi >\rho _R`$) in the neutrino decoupling epoch, as in that case there is no decoupling; rather, the universe must have attained radiation domination in the sense that $`\rho _\varphi <\rho _R`$. Consideration of the effect of scalar decay on decoupling temperatures in the regime of radiation domination leads to a lower bound of about $`10^{22}`$ Gev on the scalar decay constant $`\mathrm{\Gamma }.`$

Acknowledgements:
PA wishes to thank Dr. S. Datta for the use of facilities at Presidency College, Calcutta. DRC wishes to thank Rahul Biswas for help in accessing literature.