format: text

text:
# A Chandra X-ray Study of NGC 1068 — I. Observations of Extended Emission

## 1 Introduction

Extended X-ray emission in Seyfert galaxies represents a potentially powerful probe of these active galactic nuclei. Such X-rays could originate from either a hot, collisionally-ionized or a much cooler photoionized gas. Both may plausibly be expected to be present in these objects. The narrow line regions of Seyferts are known to be the sites of mass motions of several hundreds to $`1000\mathrm{km}\mathrm{s}^1`$. Shocks associated with these motions will generate gas with temperatures $`10^6`$$`10^7`$ K. Hot gas may also be present in the form of outflowing winds driven by radiation or radio jets from the nucleus. The compact, hard, UV – soft X-ray nuclear continuum source is believed to photoionize the narrow line region. Lines from highly ionized species, such as He ii, \[Ne v\], \[Fe vii\], \[Fe x\], \[Fe xi\] and \[Fe xiv\], are found in the optical and infrared spectra of these objects. X-ray lines from highly ionized species may, therefore, also be expected. High velocity shocks can also be powerful sources of ionizing radiation and, if present, should provide both collisionally- and photo-ionized gas. By searching for both of these components, X-ray observations can probe whether shocks are significant sources of ionizing radiation for the narrow line region. Extended X-ray emission can also arise through electron scattering and fluorescence of the nuclear radiation in extended gas. Morphological correspondences between the X-ray emitting gas and the optical line and radio continuum structures may provide clues to the nature of the X-ray emission.

For these reasons, we have begun a program of observing Seyfert galaxies with Chandra. Previous X-ray observatories have lacked the high spatial resolution ($`<1\mathrm{}`$), wide energy coverage (0.1 – 10 keV) and good spectral resolution ($`130\mathrm{eV}`$ from the CCD detectors) of Chandra. These capabilities are ideal for investigation of extended gas in Seyferts, and in this paper we present the results of observations of the nearby Seyfert 2 galaxy NGC 1068.

NGC 1068 was observed on three occasions with the Einstein observatory IPC by Monier & Halpern (1987). They found that the 0.1 – 4.5 keV spectrum could be described by a power law with energy index $`\alpha =2.0\pm 0.3`$ and absorbing column density consistent with the Galactic value. Combination of EXOSAT and Einstein IPC observations (Elvis & Lawrence, 1988) showed that the X-ray spectrum can be decomposed into two components — a steep low energy ($`<2\mathrm{keV}`$) part with $`\alpha =3.5`$ and a flat high energy (2 – 10 keV) part with $`\alpha 0.7`$. This was the first detection in a Seyfert 2 of the hard power-law X-ray source known to be present in Seyfert 1s. Any intrinsic absorption was shown to be small ($`N_H<3\times 10^{20}\mathrm{atom}\mathrm{cm}^2`$). Ginga observations (Koyama et al., 1989) confirmed the presence of the hard source above 2 keV, and also discovered an intense iron line with equivalent width $`1.3\pm 0.1\mathrm{keV}`$, in accord with the predictions of Krolik & Kallman (1987). Marshall et al. (1993) obtained a higher spectral resolution (FWHM $`100\mathrm{eV}`$) observation with BBXRT. They found a total equivalent width for the Fe K line of 2.8 keV, and modeled the line profile in terms of three components corresponding to fluorescence of neutral (and lowly ionized) iron and recombination into both He-like and H-like iron. Fe L-shell emission was also detected. In a study with ASCA, Ueno et al. (1994) found emission from the H-like and / or He-like ions of Ne, Mg, Si and S and confirmed the Fe L and Fe K features found by Marshall et al. (1993). The best fit power-law spectrum above 3 keV has energy index $`\alpha 0.3\pm 0.3`$. Ueno et al. (1994) interpreted the steep low energy ($`<3\mathrm{keV}`$) spectrum in terms of a two temperature (0.59 and 0.14 keV), optically-thin coronal plasma model. Iwasawa, Fabian & Matt (1997) reanalyzed the ASCA data, confirming the 3 components to the Fe K line complex which they interpreted with a model of cold and warm reflection. Netzer & Turner (1997) developed both photoionization and hot plasma models, including a hard reflected continuum. Their photoionized gas model is in good agreement with the emission lines observed by ASCA assuming solar metallicity for all elements except iron, which is more than twice solar, and oxygen, which is less than 0.25 solar. NGC 1068 has recently been studied with BeppoSAX. Matt et al. (1997) detected the galaxy in the 20 – 100 keV band, supporting models envisaging a mixture of both neutral and ionized reflections of an otherwise invisible nuclear source. The nucleus is completely obscured at all energies, implying it is “Compton-thick” — i.e. the column density of absorbing matter exceeds $`N_H10^{25}\mathrm{atoms}\mathrm{cm}^2`$ (see also Matt et al. 2000). Soft X-ray spectroscopy with BeppoSAX (Guainazzi et al., 1999) reveals a spectrum rich in K$`\alpha `$ transitions of various elements below 3 keV.

The only X-ray image of NGC 1068 was obtained with the ROSAT HRI by Wilson et al. (1992) with $`4\mathrm{}`$$`5\mathrm{}`$ resolution. The observed X-ray emission could be divided into three components: i) an unresolved (radius $`3\mathrm{}`$) source associated with the Seyfert nucleus, ii) resolved circumnuclear (radius $`<15\mathrm{}`$) emission extending preferentially towards the NE and iii) large-scale ($`15\mathrm{}<\mathrm{radius}<60\mathrm{}`$) emission aligned NE-SW. Wilson et al. (1992) interpreted the circumnuclear emission as thermal emission from a nuclear wind and the large-scale emission as either associated with the disk starburst or an extension of the nuclear wind to larger radii. Further ROSAT HRI observations were reported by Wilson & Elvis (1997).

Very recently, Paerels et al. (2000) have reported a high spectral resolution observation of NGC 1068 with the XMM-Newton Reflection Grating Spectrometer. Their spectrum reveals a large number of emission lines below 2 keV. Based on the narrowness of the recombination radiation continua, the fact that the forbidden lines are stronger than the resonance lines in helium-like triplets, and the weakness of Fe L emission compared to Fe K emission, Paerels et al. (2000) conclude that the emission is dominated by recombination in cool X-ray photoionized gas. The N and O helium-like triplets may contain a $`<30\%`$ component in collisional equilibrium.

The present paper is organized as follows. In section 2 we describe the observations and their reductions. Section 3 provides a brief discussion of the X-ray morphology, while section 4 presents the X-ray spectra obtained for the nucleus and five spatially extended regions. A comparison of the X-ray images with those in other wavebands is presented in section 5, and conclusions are given in section 6. Future papers will address the compact sources in the field and present a more detailed analysis of the extended emission.

NGC 1068 is at a heliocentric redshift of $`z=0.00379`$ (de Vaucouleurs et al., 1991) which, after correction to the frame of the microwave background and assuming $`H_0=50\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and $`q_0=0`$, corresponds to 22.8 Mpc, so $`1\mathrm{}`$ is equivalent to 110 pc. The Galactic column density towards NGC 1068 is $`N_H(\mathrm{Gal})=2.99\times 10^{20}\mathrm{cm}^2`$ (Murphy et al., 1996)

## 2 Observations and Reduction

Since the nucleus of NGC 1068 is known to be a strong X-ray source, we were concerned that “pile-up” could affect the Chandra observations. In order to measure the Chandra count rate from the nucleus and circumnuclear regions, we first obtained a short observation (obsid 343) of NGC 1068 on 1999 December 9 using chip S3 (backside illuminated) of the Advanced CCD Imaging Spectrometer (ACIS; Garmire et al. 2000). Single exposures with a 0.1 s frame-time were alternated with two exposures with a 0.4 s frame-time, the total integration time being $`<1`$ ks. This observation showed that the nucleus is significantly piled up in the 0.4 s frame-time, but not in the 0.1 s frame-time. It also demonstrated that the count rate in the region extending $`6\mathrm{}`$ NE from the nucleus is so high that data obtained from this region with the default frame-time of 3.2 s would be piled up, but data obtained with a 0.4 s frame-time would not be piled up. On this basis we decided to observe NGC 1068 for separate, longer exposures in all of 0.1 s, 0.4 s, and 3.2 s frame-times. The longer 0.1 s, 0.4 s and 3.2 s frame-time observations were thus designed to study primarily the X-ray emission of the nucleus, the extended region to the NE of the nucleus and the larger scale emission, respectively.

These longer integrations were taken with the ACIS instrument in two separate observations. On 2000 February 21 the standard 3.2 s frame-time was used. CCDs I2, I3, S1, S2, S3 and S4 were read out, though all of the detected X-ray emission from NGC 1068 is on S3. On 2000 February 22 the readouts were alternated between one 0.1 s and two 0.4 s frame-times in order to mitigate the effects of pile-up in the inner regions, as discussed above. The data were screened for times of high background count rates and aspect errors in the usual way. The 3.2 s, 0.4 s and 0.1 s frame-time data have total good exposure times of 47016.2 s, 11465.4 s and 1433.2 s, respectively. The “livetimes” for the alternating mode exposures were given incorrectly in the header and had to be computed separately. A corrected formula for the exposure time (G. Allen, private communication) was used, and the resulting livetimes are equal to the frame-time multiplied by the number of exposures taken with that frame-time.

Spectra and instrument responses were generated using ciao v1.1.5 and analyzed with xspec v11.0.1. Spectra were grouped to have at least 25 counts per energy bin to ensure that the $`\chi ^2`$ fitting of the data was valid. When extracting spectra from complex regions the backscal keyword (equal to the fractional area of the chip occupied by the extraction region) was computed by hand, as the values generated by ciao are often incorrect.

In the absence of the gratings, X-ray spectra from CCD detectors such as ACIS are not invertible, i.e. one cannot uniquely determine the X-ray spectrum from the observed counts. To model such data the following technique is used. A parameterized model of the source spectrum is computed and folded through the instrument response. The folded model and data are then compared using, e.g., the $`\chi ^2`$ test. Parameters of the model spectrum are adjusted to minimize the $`\chi ^2`$ difference between the folded model and the data. The error bars for a single parameter are given by the range over which that parameter may be varied with $`\mathrm{\Delta }\chi ^2=2.706`$, keeping the other model parameters fixed. This is equivalent to the 90 per cent confidence interval for a single interesting parameter.

The backside illuminated S3 chip has a high sensitivity to very soft X-rays and significant counts at $`0.1\mathrm{keV}`$ were seen from all regions. However the instrument response is uncertain below 0.50 keV, and this uncertainty increases towards lower energies. We thus restricted our modeling to photon energies above 0.25 keV.

### 2.1 Nucleus

Counts were extracted from a 3$`\stackrel{}{\mathrm{.}}`$75 diameter circle centered on the nucleus. The 3.2 s frame-time data are extremely piled up and cannot be used for spectral fitting. In the 0.4 s frame-time data there are 25251 counts corresponding to a count rate of $`2.20\pm 0.014`$ ct$`\mathrm{s}^1`$. In the 0.1 s frame-time data there are 4760 counts corresponding to a count rate of $`3.32\pm 0.048`$ ct$`\mathrm{s}^1`$. The difference in count rate is due, in part, to the presence of pile-up in the 0.4 s frame-time data. The background is negligible in both this and the NE region.

### 2.2 NE Region

Counts were extracted from a 4$`\stackrel{}{\mathrm{.}}`$1 $`\times `$ 4$`\stackrel{}{\mathrm{.}}`$6 rectangle with longer length in position angle (PA) = $`45^{}`$ centered on the extended emission 4$`\stackrel{}{\mathrm{.}}`$25 to the NE of the nucleus. This region does not overlap with the region used to extract the spectrum of the nucleus. In the 3.2 s frame-time data there are 26871 counts corresponding to a count rate of $`0.57\pm 0.003`$ ct$`\mathrm{s}^1`$. In the 0.4 s frame-time data there are 10069 counts corresponding to a count rate of $`0.88\pm 0.01`$ ct$`\mathrm{s}^1`$. In the 0.1 s data there are 1349 counts corresponding to a count rate of $`0.94\pm 0.03`$ ct$`\mathrm{s}^1`$. Thus, there is significant pile-up in the 3.2 s frame-time data in this region.

### 2.3 Larger Scale Emission

The distribution and spectrum of larger scale emission was investigated with the 3.2 s frame-time data. The large scale region was divided into four quadrants with position angles $`250^{}340^{}`$ (quadrant A), $`340^{}70^{}`$ (quadrant B), $`70^{}160^{}`$ (quadrant C) and $`160^{}250^{}`$ (quadrant D). Position angle $`25^{}`$, the bisector of quadrant B, corresponds to the approximate position angles of the radio emission and narrow line regions to the NE of the nucleus. Extraction regions between nuclear distances of $`5\mathrm{}`$ and $`60\mathrm{}`$ within quadrants A, B, C and D are termed the “West”, “North”, “East” and “South” regions respectively, and are shown in Figure 1. Obvious point sources determined by eye were excluded. Pile-up is not significant in any of these regions.

It is difficult to determine the spectrum of the background for the larger scale emission as there are no source-free regions nearby. To overcome this problem, we used a compilation of observations of relatively blank fields, from which discrete sources have been excised. For a particular source extraction region, a corresponding background region is taken from this compilation. The two regions have the same physical location on the S3 chip. The background images and software111Available at http://hea-www.harvard.edu/$``$maxim/axaf/acisbg/ of Maxim Markevitch were used for this procedure.

## 3 X-ray Morphology

Grey scale images of the 0.1 s and 0.4 s frame-time Chandra X-ray observations of NGC 1068 are shown in Figures 2 and 3 respectively. The brightest region, which we refer to as the nucleus, extends $`1\stackrel{}{\mathrm{.}}5`$ to the NE. Somewhat fainter X-ray emission extends $`6\stackrel{}{\mathrm{.}}5`$ NE of the nucleus in PA $`35^{}`$$`40^{}`$, and perpendicular to this direction from 2$`\stackrel{}{\mathrm{.}}`$5 SE to 2$`\stackrel{}{\mathrm{.}}`$5 NW of the nucleus. There is a second peak of emission 3$`\stackrel{}{\mathrm{.}}`$5 to the NE of the nucleus, with a brightness equal to approximately 15 per cent of the nucleus; this peak lies within the “NE region” (see section 2.2).

A grey scale image of the 3.2 s frame-time observation is shown in Figure 4. This image agrees well with the ROSAT HRI image (Wilson et al., 1992) when allowance is made for the lower resolution ($`4\mathrm{}`$$`5\mathrm{}`$) of the latter. Two bright “prongs” of X-ray emission extend from the bright region around the nucleus. One extends from the northern tip along PA $`330^{}`$ for 7$`\stackrel{}{\mathrm{.}}`$5, and the second extends from the southern tip along PA $`90^{}`$ for 6$`\stackrel{}{\mathrm{.}}`$5 before turning to follow PA $`5^{}`$ for 4″. Significant extended X-ray emission is seen out to at least $`60\mathrm{}`$ to the NE, $`50\mathrm{}`$ to the SW, $`20\mathrm{}`$ to the NW and $`30\mathrm{}`$ to the SE of the nucleus. The most obvious large scale structures appear to be “spiral arms” that curve to lower PA with increasing radius from the nucleus, i.e. they are “trailing” arms. They may also be seen in an image of the 3.2 s frame-time data that has been smoothed by a Gaussian of $`\sigma =0\stackrel{}{\mathrm{.}}5`$ (Figure 5). The images also show many point sources of X-ray emission associated with NGC 1068.

The X-ray morphologies in the 0.25 – 0.80 keV, 0.80 – 2.00 keV and 2.00 – 7.50 keV bands are shown in Figure 6, for the 0.1 s, 0.4 s and 3.2 s frame-time observations. The 0.25 – 0.80 keV X-ray emission has the greatest spatial extent, covering the range described above. The nucleus, NE region, the two “prongs” and the trailing arms are clearly seen. Only a few point sources, however, appear in these soft X-ray images. In the 0.80 – 2.00 keV energy band, there is still significant emission out to tens of arc-seconds, and many more point sources are visible. In the hardest energy band, 2.00 – 7.50 keV, emission extends at least $`13\mathrm{}`$ to the NE and SW of the nucleus and at least $`4\mathrm{}`$ to the NW and SE. The two “prongs” (to the N and E) do not stand out in hard X-rays, perhaps because of the lower sensitivity in this band. As the point spread function (PSF) at the location of the nucleus is circularly symmetric the observed, elongated extended hard X-rays are real. In addition, many hard point sources are clearly seen $`20\mathrm{}`$ from the nucleus.

## 4 X-ray Spectra

### 4.1 Nucleus

#### 4.1.1 Thermal Plasmas

Initially the X-ray spectrum of the nucleus was modeled as an absorbed hot plasma with solar metal abundances. We used the “mekal” model (Mewe et al., 1985, 1986; Kaastra, 1992; Liedahl et al., 1995; Arnaud & Rothenflug, 1985; Arnaud & Raymond, 1992), which is included in the X-ray spectral fitting software xspec. The 0.1 s frame-time data were used as these do not suffer from pile-up. The model does not fit the data (see table 1), so a single hot plasma may be ruled out. Such plasmas overproduce iron L emission lines (i.e. transitions to the $`n=2`$ level) between 0.7 – 1.2 keV, and underproduce the high energy continuum above a few keV. A better fit is obtained if the metal abundances are allowed to vary, but the best fit value of zero is implausible for the nucleus of a spiral galaxy.

A model consisting of two hot plasmas, with temperatures of 0.17 keV and 0.73 keV and metal abundances of 3.0 per cent of solar and 2.7 per cent of solar respectively, provides an adequate description of the spectrum below 2 keV ($`\chi ^2=139`$ for 74 degrees of freedom (d.o.f.)) Again, such low metal abundances are unlikely to be found in the inner regions of a spiral galaxy, so multiple hot plasma models are implausible for the nucleus of NGC 1068.

#### 4.1.2 Bremsstrahlung Plus Power Law Plus Emission Lines

Many lines are seen below 1 keV in the XMM-Newton RGS spectrum (Paerels et al 2000), but these are confused by the limited spectral resolution of the Chandra ACIS spectrum reported here. Individual lines are sufficiently narrow to be unresolved by ACIS. We therefore decided to take a more phenomenological approach, modeling the data by a continuum plus individual emission lines. Firstly the continuum was described by a smooth model consisting of a bremsstrahlung plus a power law, both absorbed by the Galactic column. The residuals between the data and this smooth model were inspected by eye. Emission lines were then added to the model at energies expected to reduce the value of $`\chi ^2`$. The energy and flux of these lines were free parameters. This process was repeated, adding lines until either an acceptable fit was achieved, or no further improvement was possible. Finally, our observed line energies were compared with those of the XMM-Newton RGS spectrum (a hard copy of which was kindly provided by F. Paerels) and, if agreement found, identified with lines therein.

The 0.1 s frame-time data are well described by a 0.45 keV bremsstrahlung plasma plus a hard power law of photon index $`\mathrm{\Gamma }=1.01_{1.15}^{+0.86}`$ plus a number of narrow emission lines, all absorbed by the Galactic column (last model in table 1). The spectrum and best fitting model are shown in Figure 7 and the emission lines are listed in table 2. Radiative recombination continua (RRC) are due to free-bound transitions from energy levels just above the ionization threshold; hence these lines start at the ionization threshold energy, and extend to higher energies (e.g. Paerels et al. 2000).

Our spectrum below 2 keV is generally consistent with lines seen in the XMM-Newton RGS spectrum (Paerels et al., 2000). The lines near and above 2 keV cannot be uniquely identified but are well described by neutral, hydrogen-like or helium-like lines of Si, S or Ar, with the exception of lines between 6.40 and 6.97 keV attributed to Fe. The Fe line complex is well described by a neutral Fe line with an equivalent width of 2.24 keV, although it should be noted that the true level of the continuum at the iron line energy is difficult to determine. The iron line is better resolved in the 0.4 s frame time data (see below). In view of the low signal-to-noise ratio, the line flux densities are quite uncertain.

The best fitting model has an unabsorbed 0.5 – 2.0 keV flux of $`5.8\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ corresponding to an unabsorbed rest frame 0.5 – 2.0 keV luminosity of $`3.6\times 10^{41}\mathrm{erg}\mathrm{s}^1`$ and an unabsorbed 2.0 – 10.0 keV flux of $`3.8\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ corresponding to an unabsorbed rest frame 2.0 – 10.0 keV luminosity of $`2.3\times 10^{41}\mathrm{erg}\mathrm{s}^1`$. The 2.0 – 10.0 keV flux and luminosity will have fractional errors comparable to the power law component of the model, listed in table 1.

The 0.4 s frame-time observation of the nucleus suffer from pile-up and it is difficult to determine the true shape of the continuum and the individual line fluxes. A recent model by Davis (2001), which is included as part of the Interactive Spectral Interpretation System (ISIS)222http://space.mit.edu/ASC/ISIS, includes effects of pile up in modeling spectra. We fitted the same spectral model to the 0.4 s frame-time data as was used to describe the 0.1 s frame-time data, with the continuum parameters allowed to vary and the emission line parameters fixed. The residuals between the data and this model were examined by eye and any clear discrepancies noted. The fluxes of model lines at energies near such discrepancies were then allowed to vary in an attempt to improve the fit. If this did not provide a satisfactory improvement, the line energy was also allowed to vary. If no model line existed near an obvious feature a new line was added. The reduced $`\chi ^2`$ values for these overall fits were not as good as those for the 0.1 s frame-time data because we kept most of the line parameters fixed, even though very small changes may have improved the agreement between data and model. The pile-up fraction was estimated to be 22 per cent, and the continuum parameters ware close to those found for the 0.1 s frame-time data; the bremsstrahlung component was 0.1 keV warmer, and the power law component was slightly harder (but within the 0.1 s frame-time error bars) than the parameters derived from the 0.1 s frame-time observations. These small discrepancies are probably a result of imperfect modeling of the pile-up.

Emission-line fluxes from the 0.4 s frame-time observation are listed in the lower portion of Table 2. Note that the line fluxes tabulated for the 0.4 s frame time are the new *total* fluxes of those lines that have changed or been added, and are *not* the difference in line flux between the 0.1 s and 0.4 s frame-time observations. The 0.4 s differ from the 0.1 s results in that: i) the O viii RRC line is stronger (significantly above the 0.1 s frame-time error bars), ii) an Al line at 1.72 keV is probably present, iii) an “absorption” feature is seen at 1.54 keV that coincides with a sudden change in the effective area of the telescope, and is therefore probably caused by pile-up, iv) the Ar line better describes the data if it is shifted to 2.97 keV, and v) the Fe line is resolved into a neutral (6.4 keV) and ionized component ($`6.80`$ keV), the normalization of their sum being consistent with the single Fe line seen in the 0.1 s frame-time data. The fluxes of the neutral and ionized components are comparable.

### 4.2 NE Region

The 0.4 s frame-time data do not suffer from pile-up in the NE region, and have been used to determine both the shape of the continuum and the emission lines present. The 3.2 s data suffer from pile-up but have a higher signal-to-noise ratio and have been used where possible to investigate the X-ray emission lines, especially at higher energies.

#### 4.2.1 Thermal Plasmas

The 0.4 s frame-time observations of the NE region were initially modeled by an absorbed hot thermal plasma with solar metal abundances, using the mekal model in xspec. The mekal model does not provide a good description of the data (see table 3) and may be ruled out. The model overproduces emission lines between 0.7 and 1.2 keV due to transitions to the L-shell of iron ions, and underproduces the high energy continuum above a few keV. The addition of a second hot thermal plasma with solar metal abundance improves the quality of the fit but again does not provide a good description of the spectrum (table 3).

#### 4.2.2 Bremsstrahlung Plus Emission Lines

To further investigate the emission lines in the spectrum, the continuum was modeled by a smooth function consisting of two thermal bremsstrahlung components absorbed by the Galactic column. No physical meaning should be attributed to such a model, as an absorbed bremsstrahlung plus a power law model could have been used instead. Repeating the technique used to fit the spectrum of the nucleus, narrow emission lines were then added to the model. The XMM-Newton RGS spectrum (Paerels et al. 2000) was again used as a guide for the identification of the soft X-ray lines. As noted above, the spectral resolution of the ACIS is too low to provide unique identifications of the soft X-ray emission lines, many of which are blended together in our data below 1 keV. Two bremsstrahlung components with temperatures of 0.39 keV and 2.84 keV plus the stronger emission lines seen by Paerels et al. (2000) and a line at 2.38 keV, all absorbed by the Galactic column provides a good description of the spectrum with a $`\chi ^2`$ of 78 for 62 d.o.f. (bottom model, table 3). The emission lines are listed in table 4, and the spectrum and model are shown in Figure 8.

A single data bin at 0.5 keV (see Figure 8) contributes significantly to the $`\chi ^2`$ of the fit and, if removed, the quality of fit is improved to $`\chi ^2=65`$ for 61 d.o.f. This data point is $`9\sigma `$ higher than its two neighbors over an energy range of 15 eV (Figure 8). This is much narrower than the instrumental response and the reason for this feature is unknown.

The line observed in the ACIS spectrum at 0.68 keV falls in the gap between strong lines of O viii Ly $`\alpha `$, Fe xvii and the radiative recombination continua of N vii and O vii seen in the RGS spectrum, and has been attributed to a blurred combination of these. The line at an observed energy of 0.76 keV is close to O viii Ly $`\beta `$ seen in the RGS spectrum, and has been attributed to this transition. The small discrepancy in line energy may result from blending with the RRC of O vii at 0.74 keV. The line at 1.00 keV falls between strong lines of Fe xx and Ne x seen in the RGS spectrum and may represent a combination of these. The weak line at 2.38 keV corresponds to the K$`\alpha `$ transition of sulphur in the ionization range S x – S xiv.

The best fit model to the 0.4 s frame-time data has an unabsorbed 0.5 – 2.0 keV flux of $`1.6\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, corresponding to an unabsorbed rest frame 0.5 – 2.0 keV luminosity of $`1.0\times 10^{41}\mathrm{erg}\mathrm{s}^1`$, and an unabsorbed 2.0 – 10.0 keV flux of $`2.0\times 10^{13}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, corresponding to an unabsorbed rest frame 2.0 – 10.0 keV luminosity of $`1.2\times 10^{40}\mathrm{erg}\mathrm{s}^1`$.

The 3.2 s frame-time data suffer from pile-up so we followed the same procedure as was used for the 0.4 s frame-time observation of the nucleus. The pile-up fraction was estimated to be 46 per cent, but this is uncertain since the extraction region is not a point-source, as was assumed in modeling the pile-up. The continuum parameters are different to those found for the 0.4 s frame data; the cooler bremsstrahlung component is 0.1 keV warmer, and the hotter bremsstrahlung is 6.5 keV warmer. The significantly higher temperature of the latter component may be a result of the larger energy range of the 3.2 s frame-time data — i.e. the continuum model is required to produce significant X-ray flux at higher energies. The emission line fluxes from the 3.2 s frame-time data are listed at the bottom of Table 4. Again, note that the tabulated line fluxes are the *total* line fluxes and not the differences in flux between the 3.2 s and 0.4 s frame-time data. The changes are as follows: i) the N vii line has a preferred energy of 0.54 keV which corresponds to a sudden change in the effective area of the telescope, and is therefore probably caused by pile-up, ii) the Ne ix triplet is stronger, iii) the Ne ix He $`\beta `$ line is weaker, and iv) an Fe line is tentatively detected, but its energy is poorly constrained,

### 4.3 Larger Scale Emission

#### 4.3.1 Thermal Plasmas

We followed a similar procedure to that used for the nucleus and NE region, applying it to the West, North, East and South regions in turn. The 3.2 s frame-time observations were used as they have the longest exposure time and do not suffer from pile-up. A hot plasma absorbed by the Galactic column does not provide a good description of the data from any of these regions, even if the metal abundances are allowed to vary (see table 5). The preferred values of the metal abundances are very low (less than 12 per cent solar), and the thermal model does not fit the significant hard X-ray emission seen from each region above 2 keV. We conclude that a single temperature thermal plasma model, absorbed by the Galactic column, cannot describe the spectrum in any of the larger scale extended regions.

#### 4.3.2 Bremsstrahlung Plus Power Law Plus Emission Lines

The phenomenological approach described in section 4.1.2 was again used to study the emission lines from each of these regions. The continuum was described by a smooth model consisting of a bremsstrahlung plus power law, both absorbed by the Galactic column, to which individual emission lines were then added. The parameters of the continuum fits are listed in table 5, and the emission lines listed in tables 6 (West), 7 (North), 8 (East) and 9 (South). The spectrum and models are shown in Figures 9 (West), 10 (North), 11 (East) and 12 (South). Good fits were obtained with bremsstrahlung components with temperatures in the range 0.36 – 1.06 keV and power laws with photon indices in the range $`\mathrm{\Gamma }=0.21`$ to $`0.92`$, with the harder (lower) photon indices in the South and East regions. There is significant hard X-ray emission above 2 keV from each region, and an iron line is seen from the West, North and South regions. The unabsorbed broad-band fluxes and luminosities are listed in Table 10.

The energy dependent PSF of the Chandra mirrors will scatter a fraction of the flux from the nucleus and NE region over tens of arc-seconds and this can significantly contaminate the spectra of the larger scale extended emission regions. In particular, the wings of the PSF are stronger at higher energies and it is important to assess the reliability of the hard X-ray emission seen on the larger scales. Firstly, the spectra show hard ($`>3\mathrm{keV}`$) X-ray fluxes differing by a factor of $`34`$ between the East (the weakest) and South (the strongest) regions (see Figures 11 and 12 respectively). In the regions between $`23\stackrel{}{\mathrm{.}}5`$$`60\mathrm{}`$ a similar difference in hard X-ray flux is also seen (Figure 13), with stronger emission from the West and South. As the PSF is essentially circularly symmetric these must be due to differences in the intrinsic spectra of these regions. Secondly, images of the extended hard X-ray emission (e.g. Figure 6), are clearly not circularly symmetric with the strongest emission in the large-scale regions being seen in the South and West sectors.. This strongly suggests that scattering of nuclear emission by the circularly symmetric PSF cannot account for all of this emission. An “iron line” image was constructed by taking the 6.3 – 6.5 keV iron line band and subtracting normalized 5.8 – 6.3 and 6.5 – 7.0 keV continuum band images. The resulting distribution of discrete counts from the iron line extends $`20\mathrm{}`$ to the NE and SW, and at lower surface brightness $`50\mathrm{}`$ to the West and South. This agrees well with the hard X-ray flux and iron lines seen in the spectra of the extended and far-extended regions (Figures 912 and Figure 13).

To quantify the effects of scattering by the PSF the following approach was used. Four images were constructed in which the central circular region of radius $`5\mathrm{}`$ contained the observed distribution of emission in the bands 0.25 – 0.40, 1.10 – 1.90, 3.56 – 5.46 and 5.46 – 7.34 keV. Outside this radius (which corresponds to the inner border of the large-scale regions), the emission was set to zero. These images were then convolved with PSFs computed at 0.23, 1.50, 3.51 and 6.40 keV respectively, using the ciao tool mkarf. From the convolved images, the predicted count rate per keV due to telescope scattering from the inner bright region may be obtained for each of the larger scale extended regions. These are the four horizontal bars plotted in Figures 9 to 13. There is clear evidence of an excess of hard X-ray emission from the North and South regions, and tentative evidence from the West, whereas the hard X-ray flux from the East region may have been scattered from the nucleus by the telescope mirrors. Restricting our attention to the more distant regions between $`23\stackrel{}{\mathrm{.}}5`$$`60\mathrm{}`$ (Figure 13) we find there to be evidence of excess hard X-ray emission from the West and South regions, and possibly the North.

### 4.4 The Spatial Distribution of the Absorbing Column

Spectra containing at least 800 counts were extracted from square regions of the 3.2 s frame-time Chandra data, allowing the absorbing column density to each region to be determined. Each spectrum was modeled by an absorbed mekal plasma with variable metal abundance over the energy range 0.25 – 3.00 keV. The resulting map of absorbing column density over the face of the galaxy is shown in Figure 14. The SW side of the galaxy shows absorption in excess of the Galactic column ($`N_H(\mathrm{Gal})=2.99\times 10^{20}\mathrm{cm}^2`$), whereas the NE side shows no significant excess absorption. This result shows that the extended X-ray emission to the NE of the nucleus is on the near side of the galactic disk of NGC 1068, while that to the SW is either in the galactic disk or on its far side. This geometrical arrangement agrees with indications from both the high excitation ionized gas seen at optical wavelengths and radio continuum observations. The ionized gas to the NE is much stronger than that to the SW, a difference usually attributed to excess obscuration on the SW side (e.g. Baldwin, Wilson & Whittle 1987; Evans et al. 1991). Further, the NE radio lobe is strongly polarized while that to the SW is unpolarized (Wilson & Ulvestad 1983, hereafter WU), consistent with the SW lobe being depolarized by gas in the galaxy disk. Lastly, H i 21 cm absorption observations indicate unambiguously that the NE radio lobe is located in front of the galaxy disk and the SW radio lobe is behind it (Gallimore et al., 1994).

## 5 Comparison with Observations in Other Wavebands

### 5.1 Registration of Images

NGC 1068 contains considerable structure within the central arc second. There are several radio components arranged in a “bent linear” structure (e.g. WU; Ulvestad, Neff & Wilson 1987; Muxlow et al. 1996; Gallimore et al. 1996, Gallimore, Baum & O’Dea 1997; Roy et al. 1998). There is also structure in optical emission line and continuum images obtained with HST (e.g. Evans et al. 1991; Macchetto et al. 1994), and in both near- and mid-infrared band images obtained at the diffraction limit of large ground-based telescopes (e.g. Weinberger, Neugebauer & Matthews 1999; Bock et al. 2000).

The location of the true nucleus (meaning the putative black hole and its accretion disk) has been much discussed. Radio source S1 has a flat or inverted spectrum (Gallimore et al. 1996), which may originate through thermal emission from the inner edge of the obscuring torus (Gallimore, Baum & O’Dea 1996; Roy et al. 1998). VLBA imaging (Gallimore, Baum & O’Dea 1997) shows that this source is a parsec-sized disk of ionized gas. Source S1 is also the location of water maser emission, indicative of dense molecular gas (Claussen & Lo 1986; Greenhill et al. 1996). To within measurement errors, radio source S1 is coincident with the peaks of mid-infrared (Braatz et al. 1993) and near-infrared (Thatte et al. 1997) light, and probably with the center of symmetry of the polarization pattern of scattered UV light (Capetti et al. 1995; Kishimoto 1999). There is thus a general consensus that these features mark the location of the true nucleus.

There is, however, no particular reason why the peak of X-ray emission should coincide with this location. The column density, NH, to the X-ray nucleus of NGC 1068 exceeds $``$ 1025 cm-2 (e.g. Matt et al. 2000), so all of the X-ray emission we observe must be either intrinsically extended or scattered nuclear light. The Chandra astrometry is not currently accurate enough for a precise location of the nucleus with respect to the sub arc-second scale features described in the previous paragraph. The position of the peak of the broad-band X-ray emission, as measured by the 0.4s frame-time exposure, is: $`\alpha _\mathrm{x}`$(J2000) = 02h 42m 40$`\stackrel{\mathrm{s}}{\mathrm{.}}`$71, $`\delta _\mathrm{x}`$(J2000) = –00 00 47$`\stackrel{}{\mathrm{.}}`$7. This position is remarkably close to that of the brightest radio source in the nucleus at 5 GHz (Muxlow et al. 1996): $`\alpha _\mathrm{r}`$(J2000) = 02h 42m 40$`\stackrel{\mathrm{s}}{\mathrm{.}}`$715s, $`\delta _\mathrm{r}`$(J2000) = –00 00 47$`\stackrel{}{\mathrm{.}}`$64. The difference in each coordinate is $`<`$ 0$`\stackrel{}{\mathrm{.}}`$1. This agreement is considerably better than expected given Chandra’s current absolute astrometric errors and may thus be coincidental. Nevertheless, we decided to align the peak in the Chandra image with the peak in a 5 GHz radio map with 0$`\stackrel{}{\mathrm{.}}`$4 resolution (WU), which is similar to, though somewhat better than, the resolution of the Chandra image. Simulations showed that smoothing the 60 mas resolution 5 GHz map of Muxlow et al. (1996) to the 0$`\stackrel{}{\mathrm{.}}`$4 resolution of the WU map moves the radio peak slightly to the NE (by $``$ 0$`\stackrel{}{\mathrm{.}}`$1 in each coordinate) to: $`\alpha _{\mathrm{rs}}`$(J2000) = 02h 42m 40$`\stackrel{\mathrm{s}}{\mathrm{.}}`$721, $`\delta _{\mathrm{rs}}`$(J2000) = –00 00 47$`\stackrel{}{\mathrm{.}}`$55. This position agrees to $``$ 0$`\stackrel{}{\mathrm{.}}`$1 with the actual position of the peak in the WU map after precessing the original B1950.0 coordinates of that map to J2000.0 using the facility in the ds9 display program.

Thus, for purposes of comparison with observations in other wavebands, the peak in the broad-band Chandra image was assigned the coordinates $`\alpha _{\mathrm{rs}}`$, $`\delta _{\mathrm{rs}}`$ given above. In order to extend the astrometric comparison to optical images, we used the alignment between the optical and radio emissions obtained by Capetti, Macchetto & Lattanzi (1997). These authors obtained absolute astrometry for the HST images by taking three images: a photographic refractor plate, a ground-based CCD image and an HST image. The plate covered astrometric standard stars from which astrometric positions were obtained for many fainter, secondary astrometric reference stars some of which are on the CCD image. Finally, the astrometric scale of the HST image was registered to the ground-based CCD image by a cross-correlation technique. From this procedure, the location of the peak of emission in the HST optical continuum (filter F547M) image is found to be: $`\alpha _\mathrm{o}`$(J2000) = 02h 42m 40$`\stackrel{\mathrm{s}}{\mathrm{.}}`$711, $`\delta _\mathrm{o}`$(J2000) = –00 00 47$`\stackrel{}{\mathrm{.}}`$81 (Capetti, Macchetto & Lattanzi 1997). This position is also very close to the Chandra-measured location of the X-ray peak: $`\alpha _\mathrm{x}`$$`\alpha _\mathrm{o}`$ = 0$`\stackrel{}{\mathrm{.}}`$015, $`\delta _\mathrm{x}`$$`\delta _\mathrm{o}`$ = 0$`\stackrel{}{\mathrm{.}}`$15, which represents agreement given the Chandra errors.

We emphasize that this astrometric registration between the X-ray image on the one hand and the radio and optical images on the other is not precise and may change as the Chandra astrometry improves. However, this uncertainty only affects the alignment on sub arc-second scales, such as is needed to compare the Chandra image with HST and radio interferometric data. Spatial comparisons on larger scales (i.e. $`>`$ 1′′) are robust. In the following, we compare the X-ray morphology with that in the optical line, optical continuum and radio continuum.

### 5.2 Comparison with Optical Emission-Line Images

Figure 15 shows the X-ray emission (contours) superposed on an HST image taken through filter F502N; this image is dominated by \[O iii\] $`\lambda `$5007. There is remarkable agreement between the structures in the two wavebands. To the S of the nucleus (assumed to be radio source S1, marked by the cross) and within 2′′ of it (Figure 15, right), the structure in both wavebands is aligned approximately N-S (PA $``$ 13) and there are detailed correspondences between the two images. To the N of the nucleus, both structures bend towards higher PA, reaching PA $``$ 35 some 4′′ from the nucleus. It is notable that the bright extension of the X-rays corresponds closely to the bright emission-line clouds which extend from the nucleus to $``$ 2′′ NE of it (Figure 15, right). There is also a close correspondence between the brightnesses slightly further (6′′) from the nucleus, where the X-ray contours take a V-shaped form, with apex towards the NE (see also Figures 17 and 18).

Correspondences are also seen on the larger scale shown in Figure 15, left. Some 7′′ NE of the nucleus, there is a prominent narrow feature in the HST image that first extends N and then curves towards the NW, ending up extending in PA $``$ 287 when 10′′ from the nucleus. The X-ray ridge, although observed with lower resolution, follows this feature closely. An arm-like feature in the HST image begins some 4′′ SE of the nucleus and curves towards lower PA, extending towards the NE and ending up some 6 – 7 ′′ E of the nucleus; the X-ray contours follow it closely. Various other close correspondences between the faint \[O iii\] and X-ray emission are also visible in Figure 15, left.

Similar correspondences are seen between the X-ray emission and the HST image through filter F658N (Figure 16), which is dominated by H$`\alpha `$ \+ \[N ii\] $`\lambda \lambda `$6548, 6583. The relative intensities of the various line features are different in this band (e.g. the narrow feature beginning 7′′ NE is much weaker here than in \[O iii\]) but there is still a close association between the structures in the two wavebands.

A comparison between \[O iii\] $`\lambda `$5007 and X-ray emission on a larger scale is shown in Figure 17. The lower panel shows again the extension from the nucleus to the NE which is common to the two images. The overall shape of the outer X-ray contours is similar to the \[O iii\] (top panel). To the SW, there are two spiral arms that curve outwards in a clockwise sense (i.e. they trail the galactic rotation). The inner, brighter one is some 18′′ and the outer, fainter one some 28′′ from the nucleus. In both cases, there is a remarkably precise alignment with the X-ray contours (top panel). The stronger, X-ray point sources are not seen in the \[O iii\] image; this is not surprising if they are X-ray binaries in NGC 1068.

### 5.3 Comparison with Optical Continuum Images

Figure 18 shows the Chandra X-ray contours superposed on a red continuum image (from Pogge & De Robertis 1993). Features common to both images include a protrusion 7 – 8′′ SW of the nucleus and, as described in the previous subsection, the two spiral arms SW of the nucleus and the arm-like feature SE and E of the nucleus. However, on the largest scales on which X-ray emission is seen, the X-ray image is much more elongated than the optical. There is little or no X-ray emission seen more than 30′′ from the nucleus in the SE or NW directions, while the X-rays extend much further to NE and SW. The impression gained from these comparisons with optical data is that the X-rays correlate more strongly with the optical line than the optical continuum emission.

This absence of a close correlation between the X-ray and optical continuum images strongly suggests that the starburst disk is not the dominant source of the extended X-ray emission.

### 5.4 Comparison with Radio Continuum Images

As noted in Section 5.1, radio continuum images with sub arc-second resolution reveal a “bent linear” structure extending over 1$`\stackrel{}{\mathrm{.}}`$5. Although the Chandra image of the central regions (e.g. Figure 15) does not have good enough resolution for a meaningful comparison, the extension of the X-ray emission immediately to the S of the nucleus is more or less N-S, like the radio emission. North of the nucleus, the X-ray ridge line “bends” towards the NE, as does the radio emission. Capetti, Macchetto & Lattanzi (1997) show that this arc second-scale radio structure lies in a region of relatively low optical line emission and is surrounded by emission line clouds. As argued by many authors (e.g. Wilson & Willis 1980; Wilson, Ward & Haniff 1988; Evans et al. 1991; Macchetto et al. 1994; Gallimore, Baum & O’Dea 1996), the structure of the narrow line region implies compression of interstellar gas by the radio ejecta. Thus, given the close association between the X-rays and the narrow line region (Section 5.2), the alignment between the X-rays and the arc second-scale radio structure is expected.

A comparison between the X-ray (grey scale) and radio (contours) images on a larger scale is shown in Figure 19. The V-shaped, edge-brightened radio lobe to the NE, which may represent a blast wave driven into the interstellar gas by the radio ejecta (Wilson & Ulvestad 1987), is seen to be associated with bright X-ray emission. The X-ray image also exhibits a V-shaped structure with apex some 6′′ from the nucleus (see also Figures 15, 17). This structure is obviously connected with the radio “V”, but has a wider opening angle. This last difference may result, in part, from the lower resolution of the X-ray image. To the SW of the nucleus, the X-rays are much weaker, consistent with the higher absorbing column to the X-ray emitting gas found here (Section 4.4). Nevertheless, the brighter regions of radio emission in the SW lobe do appear to be associated with brighter X-ray emission. There is little or no enhancement of the X-ray emission at the radio “hot spot” 4$`\stackrel{}{\mathrm{.}}`$0 SW of the nucleus.

## 6 Conclusions

We have obtained the highest resolution ($`<1\mathrm{}`$) image to date of the X-ray emission of NGC 1068. This image shows i) a bright nucleus close to or coincident with the brightest radio and optical emission, ii) bright emission extending $`5\mathrm{}`$ (550 pc) to the NE, iii) large-scale structure reaching at least $`1\mathrm{}`$ (6.6 kpc) to the NE and SW, including X-ray emission from spiral arms, and iv) numerous point sources, most of which are likely to be associated with NGC 1068.

X-ray spectra have been obtained for ten regions – the nucleus, the bright region several arc-seconds to the NE, four $`90^{}`$ sectors between $`5\mathrm{}`$ and $`60\mathrm{}`$ from the nucleus and four similar sectors between $`23\stackrel{}{\mathrm{.}}5`$ and $`60\mathrm{}`$ from the nucleus. Hot plasma models are poor descriptions of the spectra and the best approximations require implausibly low ($`<0.1Z_{}`$) abundances. We have constructed models of the spectra comprising two smooth continua (a bremsstrahlung plus a power-law or two bremsstrahlung spectra) plus emission lines. The lines cannot be uniquely identified with the present spectral resolution but are generally consistent with the stronger lines seen in the XMM-Newton RGS spectrum below 2 keV (Paerels et al., 2000). Above 2 keV we observe K$`\alpha `$ transitions of sulfur, argon and iron. Hard X-ray emission, including an iron line, is seen extending $`20\mathrm{}`$ (2.2 kpc) NE and SW of the nucleus. Lower surface brightness hard X-ray emission, with a tentatively detected iron line, extends $`50\mathrm{}`$ (5.5 kpc) to the west and south. Taken together with the XMM-Newton RGS spectrum, our results suggest photoionization and fluorescence of gas by radiation from the Seyfert nucleus to several kilo-parsecs from it.

An investigation of the distribution of absorbing column shows that the emission to the NE of the nucleus is absorbed by only the Galactic column and is thus on the near side of the disk of NGC 1068. The X-ray emission to the SW of the nucleus suffers greater absorption and must originate in or behind the disk of NGC 1068. This geometry agrees with that inferred for both the narrow line region observed at optical wavelengths and the radio ejecta.

We have compared the X-ray image with optical emission-line, optical continuum and radio continuum images. There is a very strong correlation between the X-ray emission and high excitation optical line emission (as traced by \[O iii\] $`\lambda 5007`$), consistent with photoionization in both wavebands. To the SW, the X-ray spiral arms correlate closely with those seen in \[O iii\] $`\lambda 5007`$. The NE radio lobe some $`6\mathrm{}`$ (660 pc) from the nucleus shows a strong morphological correlation with the X-ray emission, suggestive of interstellar gas being compressed by the radio ejecta. The alignment of the large-scale (arc min) and small scale (few arc sec) X-ray emissions and the lack of a close correlation between the larger scale X-ray emission and the optical continuum strongly suggests that the starburst is not the dominant source of these X-rays.

We plan to present a discussion of the compact X-ray sources seen in NGC 1068 and a more detailed analysis of the extended emission in future papers.

We thank the Chandra Science Center, especially Dan Harris and Shanil Virani, for assistance with the observations. Glenn Allen gave valuable advice on the alternating exposure mode. We are grateful to Gerald Cecil and Frits Paerels for providing data in advance of publication. We thank John Houck for assistance with ISIS. This research was supported by NASA grant NAG 81027.


text:
# Quantum Fractal Fluctuations

## Abstract

We numerically analyse quantum survival probability fluctuations in an open, classically chaotic system. In a quasi-classical regime, and in the presence of classical mixed phase space, such fluctuations are believed to exhibit a fractal pattern, on the grounds of semiclassical arguments. In contrast, we work in a classical regime of complete chaoticity, and in a deep quantum regime of strong localization. We provide evidence that fluctuations are still fractal, due to the slow, purely quantum algebraic decay in time produced by dynamical localization. Such findings considerably enlarge the scope of the existing theory.

The study of mesoscopic conductance fluctuations has gained important insights from the analysis of the underlying classical dynamics . In particular, it has been found that some statistical properties of quantum fluctuations are semiclassically determined by the classical law which rules the decay in time of the survival probability inside the open chaotic system. Building on this general relation, it was eventually predicted that, if the classical decay law is algebraic, due to the presence of residual stable islands in the classical phase space, then the quantum fluctuation pattern should display a fractal structure with a fractal dimension directly related to the classical decay exponent . Such fractal structure of quantum fluctuations has been interpreted as a quantal manifestation of classical fractality associated with hierarchical structures of surviving stable islands. This theoretical prediction has received experimental support, and has been numerically confirmed on a dissipative model system, where coupling to continuum is simulated by complete absorption of that part of the wave packets which propagates outside the interaction region.

In such a theoretical frame, quasi-classicality of the quantum dynamics appears essential in order that slow classical decay be quantally reproduced over sufficiently long time scales. In contrast, in a deep quantum regime, localization effects would become dominant, rapidly effacing any memory of the underlying classical dynamics.

In this Letter we provide evidence that fractality of quantum fluctuations nevertheless survives even in such strongly non-classical situations. In the model systems presently to be described, the classical phase space has no significant stable islands, the classical motion is diffusive, and the classical decay is exponential; therefore, in the quantum weakly localized regime, one observes non-fractal, Lorentz-correlated conductance fluctuations. In contrast, we work in a strongly localized regime, and observe fractal fluctuations, apparently unrelated to any classical phase-space structure. Decisive in this respect is the fact that the quantum decay, albeit highly non-classical, is algebraic due to purely quantum localization effects, as recently predicted. As a matter of fact, we find that the fractal dimension of quantum fluctuations is connected to the algebraic decay exponent in exactly the same way predicted by the semiclassical theory . It therefore appears that the quantum fluctuation pattern bears the same relation to the quantum probability decay, which was originally predicted by semiclassical arguments, even though such decay is not by any means quasi-classical.

We consider the kicked rotator model  with absorbing boundary conditions. Classical dynamics is described by the standard map:

$$\overline{I}=I+k\mathrm{sin}\theta ,\overline{\theta }=\theta +T\overline{I},$$
(1)
with absorption for $`\overline{I}<0`$ and for $`\overline{I}>N`$. Classical dynamics only depends on the scaling parameter $`K=kT`$. The corresponding quantum map is obtained by the substitution $`I\mathrm{}\widehat{n}=i\mathrm{}/\theta `$. We denote $`\psi _t`$ the state vector immediately before the $`t`$th kick, $`T`$ the kicking period, and we set $`\mathrm{}=1`$. Throughout this Letter, we will treat the kick counter $`t`$ as a discrete time variable in units of $`T`$. The discrete-time evolution from time $`t`$ to time $`t+1`$ is described in the classical case by (1) and in the quantum case by :

$$\psi _{t+1}=\widehat{P}\stackrel{~}{\psi }_{t+1}=\widehat{P}e^{iT(\widehat{n}+\varphi )^2/2}e^{ik\mathrm{cos}\widehat{\theta }}\psi _t.$$
(2)
The state $`\stackrel{~}{\psi }_{t+1}`$ is obtained from $`\psi _t`$ by means of the unitary one-period propagator of the kicked rotor. The operator $`\widehat{P}`$ then projects $`\stackrel{~}{\psi }_{t+1}`$ over states in a fixed interval of $`n`$ values ($`0nN`$), so it describes complete deletion (absorption) of the part of the wave packet which propagates outside the given interval. The parameter $`\varphi `$ can be interpreted as an Aharonov-Bohm flux through the ring parametrized by the coordinate $`\theta `$.

We consider the case $`K=kT=7`$, $`k=5`$, with absorption for $`n<0`$ and $`n>N=250`$. With such parameter values, the classical phase space has no significant island of stability . We consider a statistical ensemble of orbits starting at $`t=0`$ with $`I=0`$ and randomly distributed phases $`\theta `$. The survival probability $`P(t)`$ is the fraction of the ensemble which hasn’t been absorbed at time $`t`$; equivalently, $`P(t)`$ is the integrated distribution of exit times of orbits in the ensemble.

The classical decay of $`P(t)`$ is well described by a diffusion process, practically insensitive to the chosen value of $`\varphi `$. As shown in Fig.1, after a diffusion time $`t_DN^2/DN^2/k^2`$ ($`D`$ diffusion rate), the survival probability exponentially decays with time: $`P(t)\mathrm{exp}(\gamma t)`$, with $`\gamma 1/t_D2.2\times 10^3`$. In contrast, quantum dynamics is strongly localized, as the localization length $`\mathrm{}k^2`$ is much less than the “sample size” $`N`$. Such strong localization causes a slow decay of the quantum survival probability $`P_\varphi (t)`$ The latter is defined as the total probability on states $`0<n<250`$ at time $`t`$, given that at $`t=0`$ the state was $`n_0=0`$. After the “Heisenberg time” $`t_H`$, it decays proportional to $`1/t`$ up to a time $`t_{max}\mathrm{exp}(2N/\mathrm{})`$ , as shown in Fig.1 (in this case, $`t_{max}5\times 10^8`$).

Though the exponent $`1`$ of algebraic quantum decay is independent of the value of $`\varphi `$, the quantum survival probability $`P_\varphi (t)`$ at a fixed time $`t`$ is sensitively dependent on $`\varphi `$. The fluctuation pattern exhibited by the graph of $`P_\varphi (t)`$ versus $`\varphi `$ at fixed $`t`$ (two examples are given in Fig.2) is precisely the object of our analysis, at values of $`t`$ up to $`10^4`$, and for $`10^4`$ values of $`\varphi `$ in the interval $`[0,0.1]`$. First we have computed autocorrelation functions:

$$C(\delta \varphi )=<P_\varphi (t)P_{\varphi +\delta \varphi }(t)>_\varphi .$$
(3)
Such correlations play a key role in the existing semiclassical theory of fractal fluctuations, which consists of three steps:

(i) correlation functions are related via Fourier transform to the classical decay of the canonical variable conjugate to $`\varphi `$, notably, to the function $`P_{cl}(\mathrm{\Theta })`$ yielding the probability distribution of the angle $`\mathrm{\Theta }`$ accumulated by the rotor at the time of exit from the finite sample (“accumulated” means multiples of $`2\pi `$ included). In particular, if $`P_{cl}(\mathrm{\Theta })`$ decays as $`\mathrm{\Theta }^\alpha `$ at large $`\mathrm{\Theta }`$, then

$$C(\delta \varphi )C(0)\text{const}|\delta \varphi |^\alpha $$
(4)
at small $`\delta \varphi `$.

(ii) It is physically intuitive, and numerically confirmed , that the decay exponent $`\alpha `$ coincides with the decay exponent of the classical survival probability.

(iii) a signal-theoretic argument  allows to conclude from (i) and (ii) that the fluctuation graph has a fractal dimension given by

$$D=2\frac{\alpha }{2}(\text{for }\alpha 2),$$
(5)
where $`\alpha `$ is the time decay exponent of the classical survival probability. Step (iii) actually invokes certain statistical properties of fluctuations, which were assumed to be ensured by the chaotic nature of classical dynamics.

In our case, semiclassical arguments are invalid, because of strong quantum localization. In addition, there is no classical decay exponent $`\alpha `$, because classical decay is exponential. So the above theory doesn’t apply; still, numerical results to be described below demonstrate that (i), (ii) and (iii) are still legitimate in our highly non semiclassical case, too, provided one replaces the classical distribution functions $`P(t)`$,$`P_{cl}(\mathrm{\Theta })`$ by the corresponding quantum distributions $`P_\varphi (t)`$,$`P_q(\mathrm{\Theta })`$, and $`\alpha `$ by the exponent of the quantum algebraic decay. In order to compute $`P_q(\mathrm{\Theta })`$ we perform a lifting of our model (2), that is, we translate the model (2) of a kicked particle on a circle into the model of a kicked particle on the line, described by the coordinate $`\mathrm{\Theta }`$. To simulate the corresponding dynamics we restrict $`\mathrm{\Theta }`$ inside a large, yet finite box, with periodic boundary conditions. At every evolution step $`t`$ we compute the function $`|\psi _t(\mathrm{\Theta })\stackrel{~}{\psi }_t(\mathrm{\Theta })|^2`$, which yields the probability loss occurring at time $`t`$ and position $`\mathrm{\Theta }`$ due to absorption. Summing this quantity over $`t`$ gives the total probability lost at position $`\mathrm{\Theta }`$, that is, the probability distribution of the angle accumulated at the exit time. The corresponding integrated distribution is shown in Fig.3 and displays a decay $`\mathrm{\Theta }^1`$. The faster decay in the rightmost part of Fig.3 is due to the finite total time and to the finite $`\mathrm{\Theta }`$box used in the integration. Thus, at any given time $`t`$, the distribution of accumulated exit angles up to that time can be assumed to decay proportional to $`\mathrm{\Theta }^1`$ in a range $`\mathrm{\Theta }_{min}<\mathrm{\Theta }<\mathrm{\Theta }_{max,t}`$, with $`\mathrm{\Theta }_{max,t}`$ increasing with $`t`$ as long as $`t<t_{max}`$. One can estimate $`\mathrm{\Theta }_{min}`$ as the angle accumulated until the time $`t_H\mathrm{}`$. This gives $`\mathrm{\Theta }_{min}_0^{t_HT}\sqrt{(D/T)\tau }𝑑\tau (2/3)t_H^2T`$. Assuming a maximal momentum $`\mathrm{}`$ for the rotor started at $`n_0=0`$, a rough estimate for $`\mathrm{\Theta }_{max,t}`$ is provided by $`\mathrm{}tT`$ ($`\mathrm{\Theta }_{min}2\times 10^2`$ and $`\mathrm{\Theta }_{max,t}2.1\times 10^4`$ for the case of Fig.3).

Numerically computed correlations (3) shown in the inset of Fig.4, at two different times, indicate that relation (4) is still valid in our strongly localized case. We have computed the fractal dimensions of the corresponding graphs of $`P_\varphi (t)`$ versus $`\varphi `$ shown in Fig.2 using the box-counting algorithm described in . The results are presented in Fig.4 and demonstrate that the graph has a fractal dimension $`1.5`$ over a significant range of $`\delta \varphi `$ scales, as predicted by equation (5). This range roughly lies in between scales $`\delta \varphi _{min,t}`$, $`\delta \varphi _{max}`$. The latter scale is inverse to $`\mathrm{\Theta }_{min}`$. The former scale can be estimated by $`\mathrm{\Theta }_{max,t}^1`$ (hence it decreases as $`t`$ increases), as long as it remains above a minimal $`\delta \varphi _02\pi /(N^2T)`$ scale. This scale is imposed by the finiteness of the sample, and can be explained as follows. The evolution operator (2) has complex eigenvalues inside the unit circle. Varying $`\varphi `$ causes these eigenvalues to move, and their maximal shift under a change $`\delta \varphi `$ is estimated by $`TN\delta \varphi `$. Any time one moving eigenvalue comes close to the unit circle, a local maximum in the graph of $`P_\varphi (t)`$ vs $`\varphi `$ appears. The number of peaks produced by a single moving eigenvalue as $`\varphi `$ changes from $`0`$ to $`1`$ is at most $`TN/(2\pi )`$. As there are $`N`$ eigenvalues, the total number of peaks is $`TN^2/(2\pi )`$, with an average $`\varphi `$ spacing $`\delta \varphi _0=(2\pi )/(N^2T)`$. As the fractal structure of the graph is due to the superposition of many tiny peaks originated in this way, no fractality can be expected below the $`\delta \varphi _0`$ scale.

Further confirmation for the validity of the above illustrated scenario was obtained by analyzing the following variant of the basic model (2):

$$\psi _{t+1}=\widehat{P}e^{iT(\widehat{n}+\varphi )^2/2}e^{ik|\mathrm{cos}\widehat{\theta }|}\psi _t.$$
(6)
which differs from (2) because the kicking potential has now a discontinuous derivative. Localization is in this case algebraic, as eigenfunctions decay like $`1/n^2`$. The argument used to predict the $`t^1`$ decay of the survival probability in the case of exponential localization  yields an asymptotic decay $`t^{3/4}`$. The numerical analysis illustrated above for the case of exponential localization was replicated for this model, too, yielding the results shown in Fig.5. Fractal analysis of the fluctuation graph yields good agreement over a broad range with the theoretical $`D`$ value of $`13/8`$ predicted using equation (5) and the decay exponent $`\alpha =3/4`$.

In conclusion, we have provided numerical evidence for fractal fluctuations of the quantum survival probability in a classically chaotic system in the regime of strong localization, and in the absence of significant classical critical structures. Whereas the original theory  predicted fractal conductance fluctuations, our present results are about fluctuations of a different quantity, namely survival probability. As both quantities ultimately reflect the fluctuations of the scattering matrix, their fluctuations should be similar in nature. As a matter of fact, the very same semiclassical arguments used to predict conductance fluctuations could be adapted to the survival probability, too . The main difference is that a scan of $`P_\varphi (t)`$ at finite however large $`t`$ cannot achieve the $`\varphi `$resolution exhibited by conductance at fixed (quasi-)energy. On account of the effect of increasing $`t`$ reported earlier in this Letter, such coarsening is unlikely to “fractalize” otherwise non-fractal fluctuation patterns. In any case, survival probability is a meaningful quantity in a broader class of problems than electronic transport in semiconductor structures. The scenario we have analyzed in this Letter can also be exported to other realistic problems, where ’conductance’ is instead a problematic concept. As an example we quote microwave ionization of Rydberg atoms . This may open a way for experimental observation of localization-induced fractal fluctuations in atomic physics.

Our results show that quantum fluctuations may be fractal even in situations where the existing semiclassical theory does not apply. This does not command a reinterpretation of existing experimental data, which were obtained in situations far from the strongly localized regime considered in this Letter. Nevertheless it signals that the current understanding of the deep quantum mechanisms responsible for fractal fluctuations is still far from complete.

Support from the Progetto Avanzato INFM “Quantum transport and classical chaos” is gratefully acknowledged.


text:
# Production and collective behavior of strange particles in Au + Au collisions at 2–8 AGeV

## Abstract

The E895 experiment at the AGS measured strange particle production and collective behavior in Au+Au collisions between 2–8 AGeV. The production of $`\mathrm{\Lambda }`$ Baryons and K0 Mesons as function of energy rises smoothly and exhibits a nonlinear impact parameter dependence. Neutral and positively charged Kaons exhibit a strong anti-flow behavior. $`\mathrm{\Lambda }`$ Baryons show a smaller flow signal than protons.

The production of strange particles in relativistic heavy ion collisions is an important probe for high density nuclear matter. Of particular interest are suggestions that strange particle yields can be used to investigate the nuclear equation of state  or possible signatures of a change of the in-medium mass . The flow effects of strange particles offer a good probe for their in-medium potentials .

Motivated by these suggestions the E895 collaboration has performed an extensive set of measurements at the Alternating Gradient Synchrotron (AGS) at the Brookhaven National Laboratory. These measurements were performed for Au+Au collisions in a beam energy range of 2–8 AGeV. Charged reaction products were detected in the EOS Time Projection Chamber (TPC)  situated in a uniform magnetic field. The TPC provides continuous 3D tracking and particle identification for charged particles (-1 $``$ Z $``$ 6) with full azimuthal coverage. The method of identification via the rigidity of the track and the energy loss in the TPC gas leads to ambiguities for certain rigidities, thus restricting the identification of charged kaons to low transverse momenta at backward rapidities. This restriction does not affect our neutral strange particles measurements ($`\mathrm{\Lambda }`$ Baryons and K$`{}_{s}{}^{}{}_{}{}^{0}`$ Mesons) which were reconstructed from their charged decay products by means of a neural network .

Fig. 1 shows the mean multiplicities of $`\mathrm{\Lambda }`$ baryons and K0 Mesons in central events corrected for acceptance and reconstruction efficiency. The yields show a monotonic increase with beam energy and are roughly reproduced by the RQMD  Model. They exhibit no obvious features which would support a threshold phenomenon. RQMD fails to reproduce the impact parameter dependence at 6 AGeV depicted in Fig. 2. The data exhibit a strong nonlinear dependence which was previously observed for charged Kaons for Au + Au collisions in the same energy range.

The sideward flow analysis used the procedure outlined in  to reconstruct the reaction plane for each event with the necessary corrections for the reaction plane dispersion. Subsequently the projection of the momentum of each strange particle on its associated reaction plane vector was obtained. The average momentum transfer ($``$px$``$) into the reaction plane for protons and $`\mathrm{\Lambda }`$ Baryons is shown in Fig. 3. The $`\mathrm{\Lambda }`$ Baryons exhibit a smaller sideward flow than protons which where sampled from the same events in which a $`\mathrm{\Lambda }`$ Baryon was reconstructed. This condition results in a bias of the proton sideward flow towards more central events than reported in . A comparison to a previous measurement in the Ni+Cu System at 2.0 AGeV  shows an increase of the flow with system size similar to protons . If flow effects of $`\mathrm{\Lambda }`$ Baryons are only due to interactions of its non-strange quark constituents, one would expect a 2/3 ratio of the strength of the $`\mathrm{\Lambda }`$/proton flow  shown in Fig. 3. This is indeed the case for the 2 AGeV data but it deviates significantly at higher energies, indicating that there may be additional contributions by the strange quark.

The K$`{}_{s}{}^{}{}_{}{}^{0}`$ Mesons depicted in Fig. 4 exhibit a strong anti-flow signal which increases strongly with beam energy. The RQMD model, which does not contain a kaon potential, can serve as a baseline for the effects of scattering. Its failure to reproduce even the sign of the flow is a strong hint for the existance of a repulsive Kaon potential. Invoking a strong density dependence of the Kaon potential, a calculation with ART is able to reproduce the sign and the strength of the K$`{}_{s}{}^{}{}_{}{}^{0}`$ Meson flow . Since the K$`{}_{s}{}^{}{}_{}{}^{0}`$ Mesons at these energies originate mainly from K0 Mesons – according to RQMD the K̄0/K0 ratio is smaller than 0.1 – one would expect a similar flow pattern for the K+ Mesons. Indeed, even though the acceptance for charged Kaons in E895 is very limited (\[yLab/yc.m.-1\] $`<`$-0.6) with a rapidity dependent low and high pt cut, the K+ Meson shows a clear anti-flow signal displayed in Fig. 4 which also increases with beam energy. The flow signal of the K- Mesons is compatible with no flow.

The E895 collaboration has measured neutral strange particle production and flow at AGS energies. The results indicate that production yields are consistent with existing systematics and show no apparent dramatic increase with beam energy. The transverse momentum analysis shows that $`\mathrm{\Lambda }`$ Baryons flow with protons but with smaller magnitude. K0 and K+ Mesons follow an anti-flow pattern indicating a repulsive Kaon potential.


text:
# On the Creation of Universe out of Nothing

## I  <br>Introduction

The science of cosmology has progressed up to a point where it is possible to make a valid picture of the first moment of the Universe. It is widely accepted that the creation of the Universe should be attributed to Quantum fluctuations of the vacuum. As a consistent Quantum gravity theory is not yet available, cosmologists employ General Relativity, or other classical theories, in order to study the evolution of the Universe after Planck’s time ( $`t10^{43}s`$ $`)`$ $`.`$

The present authors will try to offer now a Classical picture covering the creation of the Universe. The only thing that we need to admit, is that Einstein´s field equations yield the average values for the quantities that, in the Quantum Universe, when $`t<10^{43}s`$ , fluctuate quantum-mechanically around those average values, somehow like the Path Integral theory of FeynmanFey admits paths that fluctuate around the average trajectories given by the Classical theory. Throughout this paper, when we refer to the Quantum Universe, all quantities should be understood as given by the most probable values of those quantities, even if this is not explicitly stated.

Along a rather similar situation, Kiefer, Polarski and Starobinsky Staro  have shown that although according to the inflationary scenario for the very early Universe, all inhomogeneities in the Universe are of genuine quantum origin, clearly no specific quantum mechanics properties are observed when looking at this inhomogeneities and measuring them. It looks like we can get along without a Quantum gravity theory, and still ”explain” the origin of the Universe.

Berman Berman2009 , shows how to calculate by means of pseudo-tensors, the zero-energy of the Universe. His calculations date the year 1981 Berman1981 . References on zero energy calculations for the Universe, are given by Berman2009 . It was also shown by Cooperstock and Israelit, later, CooIs  that the total energy of the Robertson-Walker’s Universe is zero, without employing pseudo-tensors. The total energy is the sum of the positive matter energy plus the negative energy of the gravitational field, and it is zero. Thus, in the very creation moment there was no supply of energy to the Universe, because it was not needed (see Feynman10 )

## II The EoS for the very early Universe

We suppose that the Universe obeyed Einstein’s field equations, and Robertson-Walker’s metric. As the field equations yield an expansion, the most probable value of $`R`$ (the scale factor) increases, beginning with zero value. We may suppose that the Universe obeyed the relation

$$GM=\gamma R$$
(1)
where $`\gamma `$ is a constant and $`M`$ is the mass of a large ”sphere” of radius $`R`$. When $`t=0`$ , $`R=0`$ and, thus, $`M=0`$ . As stated in next section, when relation (1) is plugged with Planck ’s length $`R_{pl}`$ and mass $`M_{pl}`$ , we find $`\gamma =1`$ $`.`$ Note that, when $`\gamma 1`$ , we have a Machian UniverseWein . So, our Universe needs no initial mass, at least if it is Machian.

Let us consider the balance of energy equation Adler :

$$p=\frac{dM}{dV}$$
(2)
Using (1) and (2) Berman and Marinho Jr. Berman  have shown how the very early evolution of such a Universe obeyed the equation of state

$$p=\frac{1}{3}\rho $$
(3)
where $`V=\alpha R^3`$ , ($`\alpha `$=const) , $`p`$ and $`\rho M/V`$ stand for cosmic pressure and energy density.

We obtain relation (3) above, by finding, from (1), and (2) successively:

$$p=\frac{dM}{dV}=\frac{dM}{dR}\frac{dR}{dV}=\frac{1}{3\alpha R^2}\frac{dM}{dR}=\frac{\gamma }{3\alpha G}R^2$$

while

$$\rho =\frac{M}{V}=\frac{\gamma R}{G\alpha R^3}=\frac{\gamma }{\alpha G}R^2$$

and then, by comparing the expressions for $`p`$ and $`\rho ,`$ relation (3) is obtained. Berman Berman2006 Berman2006a , have tied the zero total energy of the Universe to the Machian calculation of energy density and cosmic pressure done above.

We shall now show that the solution of the Einstein’s field equations, in the $`\mathrm{\Lambda }0`$ case, yields the exponential law of expansion, which in turn accounts for the “inflationary” phase devised by Guth Guth  and others.

## III  Solution for the Scale-Factor

Berman and Marinho Jr. Berman  admitted that, for the very early Universe, the Machian hypothesis $`\frac{GM}{R}=1`$ , which is exactly valid, when mass $`M`$ and scale factor $`R`$ are plugged with Planck ’s mass and length ASTRO  , could point out to an equation of state for the very early Universe, in the classical domain, and thus, not be just a coincidence , only valid for the Planck ’s Universe. The resulting equation of state, that would be valid in the classical Universe domain, just after the Quantum State phase, was derived as equation (3), along with:

$$\rho =(\alpha GR^2)^1.$$
(4)
where $`\alpha =const`$ , that relates volume to $`R^3`$ , by means of

$$V=\alpha R^3$$
(5)
In (4), $`G`$ stands for Newton ’s gravitational constant, while, in (3), $`p`$ and $`\rho `$ stand respectively for cosmic pressure and energy density.

We now propose that these results be accepted also as yielding the most probable values of the quantities they represent, even in the Quantum domain.

In Berman , Berman and Marinho Jr. considered the $`\mathrm{\Lambda }=0`$ case, obtaining, from Einstein ’s field equations, and Robertson-Walker’s metric, the relation

$$R=\sqrt{\left(\frac{8\pi }{3\alpha }k\right)}.t$$
(6)
where $`k=0,`$ $`\pm 1`$ is the tricurvature.

For the $`\mathrm{\Lambda }0`$ case, Einstein’s equations read, for a perfect fluid:

$$\frac{8\pi G}{3}\rho +\frac{\mathrm{\Lambda }}{3}=+\frac{k}{R^2}+\frac{\left(\dot{R}\right)^2}{R^2}$$
(7)
and also

$$8\pi Gp\mathrm{\Lambda }=\frac{2\ddot{R}}{R}\frac{(\dot{R})^2}{R^2}\frac{k}{R^2}$$
(8)
By assuming equation of state (3), and relation (4), we find the following solution

$$R(t)=e^{\beta t}e^{\beta t}$$
(9)
where $`\beta =const`$ (to be determined) .We can check that $`R(0)=0`$ $`.`$

On plugging back into equations (7) and (8), we find

$$\beta =\frac{1}{2}\sqrt{\left(\frac{8\pi }{3\alpha }k\right)}.$$
(10)
We can also obtain a value for $`\mathrm{\Lambda },`$

$$\mathrm{\Lambda }=\frac{3}{4}\left(\frac{8\pi }{3\alpha }k\right)$$
(11)
When $`\mathrm{\Lambda }`$ goes to zero, (or $`\beta t`$ is small, in (9)) we recover result (6). Though we worked with $`\gamma =1,`$ (see equation 1), similar results are obtained with $`\gamma 1`$ $`,`$ if $`\gamma `$ is still a constant. Notice that we have a graceful entrance into exponential inflationary phase, just as we enter the Classical domain. In fact, the negative exponential part of the scale-factor , can be neglected with the passage of time. For inflationary scenario, we suggest the updated book by Weinberg Wein2008 .

## IV Discussion of Results and Conclusion

Our result for the scale factor, was adopted by Berman Berman2009a Berman2009b  in a different way: here, we produce it as a solution of the field equations with the given equation of state; there, he began to write the scale-factor equation, and , then , he found the most general equation of state that could be found for that scale-factor. It turned that the equation of state (3) is a particular one in the new solution.

As far we know, our result is so good that it reduces to a known result (relation (6)) in the limit $`\mathrm{\Lambda }=0`$ $`,`$ so we can be confident on it. We have also shown that the equation of state (3) is valid not only for $`\gamma =1`$ , as in reference , but for any finite fixed value of $`\gamma `$ $`.`$We remember that $`\mathrm{\Lambda }`$ models were considered earlier by Berman and Som Som , and othersWein2008 .

We explained how, in the origin of the Universe, a ”Classical”, picture, given by the most probable values for the Quantum, picture, coinciding with their classical functions, does not require input of initial energy nor initial mass; just, we need valid the Einstein field equations and R.W’s metric . Thus, we have a consistent picture of the creation of the Universe, out of ”nothing”.

On the other hand we have found an equation of state and a solution for the scale factor in the $`\mathrm{\Lambda }`$ $`0`$ case that implies the inflationary phase, automatically, for the very early Universe.

A first version of this paper was originally published in 2001 (Berman and Trevisan, 2001) BermanTrevisan2001 .

## V Acknowledgments

M.S.B thanks support of Geni, Albert and Paula; and the careful typing by Marcelo Fermann Guimarães.

## VI References


text:
# Classification of the pairing transition in finite fermionic systems

## Abstract

In this work we employ a simple pairing interaction model in order to study and classify an eventual pairing phase transition in finite fermionic systems. We show that systems with as few as $`1016`$ fermions can exhibit clear features reminiscent of a phase transition. To classify the nature of the transition we apply two different numerical methods, one based on standard thermodymanics, and another based on a recently proposed scheme by Borrmann et al. \[Phys. Rev. Lett. 84, 3511 (2000)\]. The transition is shown to be of second order, in agreement with results for infinite fermionic systems.

PACS number(s): 05.20.-y, 05.70.Ce, 05.70.Fh, 21.10.Ma

The standard BCS theory has been widely used to describe systems with pairing correlations and phase transitions to a superconducting phase for large systems, from the solid state to nuclear physics, with neutron stars as perhaps the largest object in the universe exhibiting superfluidity in its interior. An eventual superfluid phase in a neutron star will condition the neutrino emission and thereby the cooling history of such a star, in addition to inducing mechanisms such as sudden spin ups in the rotational period of the star; see, for example, Ref.  for a recent review. For an infinite system, such as a neutron star, the nature of the pairing phase transition is well established as second order.

When a system of correlated fermions such as electrons or nucleons is sufficiently small, the fermionic spectrum becomes discrete. If the spacing approaches the size of the pairing gap, superconductivity is expected to break down ; however, recent experiments on superconducting ultrasmall aluminum grains by Tinkham et al.  revealed the existence of a spectroscopic gap larger than the average electronic level density. This feature was interpreted as a reminiscence of superconductivity and renewed the interest  in studies of what is the lower size limit for superconductivity.

Other finite fermionic systems such as nuclei are expected to exhibit a variety of interesting phase-transition like phenomena, like the disappearence of pairing at a critical temperature $`T_c0.51`$ MeV or the nuclear shape transitions of deformed nuclei associated with the melting of shell effects at $`T_c112`$ MeV. In recent theoretical and experimental studies  of thermodynamical properties of finite nuclei, the heat capacity has been found to exhibit a non-vanishing bump at temperatures proportional to half the pairing gap. These bumps were interpreted as signs of the quenching of pair correlations, representing in turn features of the pairing transition for an infinitely large system.

In the study of phase transitions in e.g., solid state, nuclear and high energy physics, it is important to know whether a given transition really is of first order, discontinuous, or if there is a continuous change in a physical quantity like the mean energy, as in phase transitions of second order. If one works in the canonical or grand canonical ensembles, for finite systems it is rather difficult to decide on the order of the phase transition. This is due to the fact that in ensembles like the canonical, any anomaly is smeared over a temperature range of $`1/N`$, $`N`$ being the number of particles. In the analysis of finite systems, both a $`\delta `$-function peak and a power law singularity sharpen as the number of particles is increased, making it difficult to distinguish between the two cases, see, for example, Ref. .

In addition, first order phase transitions in finite systems have recently been inferred, theoretically and experimentally, from observed negative heat capacities that are associated with anomalous convex intruders in the entropy versus energy curves, resulting in backbendings in the caloric curves; see, for example, Refs. . Negative heat capacities are often claimed to appear only in calculations done in the microcanonical ensemble and are thought to vanish in the canonical or grand-canonical ensembles.

We aim in this Letter to identify the nature of the pairing transition, if any, in a finite fermionic system. Since we are dealing with pairing correlations, our Hamiltonian is

$$H=\underset{i}{}\epsilon _ia_i^{}a_iG\underset{ij}{}a_i^{}a_{\overline{ı}}^{}a_{\overline{ȷ}}a_j,$$
(1)
where $`a^{}`$ and $`a`$ are fermion creation and annihilation operators, respectively. The indices $`i`$ and $`j`$ run over the number of levels $`L`$, and the label $`\overline{ı}`$ stands for a time-reversed state. The parameter $`G`$ is the strength of the pairing force while $`\epsilon _i`$ is the single-particle energy of level $`i`$.

We assume that the single-particle levels are equidistant with a fixed spacing $`d`$. Moreover, in our simple model, the degeneracy of the single-particle levels is set to $`2J+1=2`$, with $`J=1/2`$ being the spin of the particle. Introducing the pair-creation operator $`S_i^+=a_{im}^{}a_{im}^{}`$, one can rewrite the Hamiltonian in Eq. (1) as

$$H=d\underset{i}{}iN_iG\underset{ij}{}S_i^+S_j^{},$$
(2)
where $`N_i=a_i^{}a_i`$ is the number operator. Seniority $`𝒮`$ is a good quantum number and the eigenvalue problem can be block-diagonalized in terms of different seniority values. Loosely speaking, the seniority quantum number $`𝒮`$ is equal to the number of unpaired particles.

For systems with less than $`1618`$ particles, this model can be diagonalized exactly, and we can obtain all eigenstates. In our studies below, we will always consider the case of half-filling, i.e., equally many particles and single-particle levels. This case has the largest dimensionality: for 16 particles in 16 doubly degenerate single-particle shells, we have a total of $`4\times 10^8`$ states. We choose units MeV for the energy and set $`G=0.2`$ MeV in all calculations while we let $`d`$ vary.

Through diagonalization of the above Hamiltonian we can define exactly the density of states $`\mathrm{\Omega }_N(E)`$ for an $`N`$-particle system with excitation energy $`E`$. An alternative to the exact diagonalization, would be to use Richardson’s well-known solution , however, we are interested in all eigenstates, and the amount of numerical labor will most likely be similar. The density of states is an essential ingredient in the evaluation of thermal averages and for the discussion of phase transitions in finite systems. For nuclei, experimental information on the density of states is expected to reveal important information on nuclear shell structure, pair correlations and other correlation phenomena in the nucleonic motion.

The density of states $`\mathrm{\Omega }_N(E)`$ is the statistical weight of the given state with excitation energy $`E`$, and its logarithm

$$S_N(E)=k_B\mathrm{ln}\mathrm{\Omega }_N(E),$$
(3)
is the entropy (we set Boltzmann’s constant $`k_B=1`$) of the $`N`$-particle system. The density of states defines also the partition function in the microcanonical ensemble and can be used to compute the partition function $`Z`$ of the canonical ensemble through

$$Z(\beta )=\underset{E}{}\mathrm{\Omega }_N(E)e^{\beta E},$$
(4)
with $`\beta =1/T`$ the inverse temperature. With $`Z`$ it is straightforward to generate other thermodynamical properties such as the mean energy $`E`$ or the specific heat $`C_V`$.

The density of states can also be used to define the free energy $`F(E)`$ in the microcanonical ensemble at a fixed temperature $`T`$ (actually an expectation value in this ensemble),

$$F(E)=T\mathrm{ln}\left[\mathrm{\Omega }_N(E)e^{\beta E}\right].$$
(5)
Note that here we include only configurations at a particular $`E`$.

The above free energy was used by e.g., Lee and Kosterlitz , based on the histogram approach for studying phase transitions developed by Ferrenberg and Swendsen , in their studies of phase transitions of classical spin systems. If a phase transition is present, a plot of $`F(E)`$ versus $`E`$ will show two local minima which correspond to configurations that are characteristic of the high and low temperature phases. At the transition temperature $`T_C`$ the value of $`F(E)`$ at the two minima equal, while at temperatures below $`T_C`$, the low-energy minimum is the absolute minimum. At temperatures above $`T_C`$, the high-energy minimum is the largest. If there is no phase transition, the system developes only one minimum for all temperatures. Since we are dealing with finite systems, we can study the development of the two minima as function of the dimension of the system and thereby extract information about the nature of the phase transition. If we are dealing with a second order phase transition, the behavior of $`F(E)`$ does not change dramatically as the size of the system increases. However, if the transition is first order, the difference in free energy, i.e., the distance between the maximum and minimum values, will increase with increasing dimension.

To elucidate the nature of the pairing transition we employ two different approaches, which both rely on our ability of computing the exact density of states $`\mathrm{\Omega }(E)`$.

First, we calculate exactly the free energy $`F(E)`$ of Eq. (5) through diagonalization of the pairing Hamiltonian of Eq. (1) for systems with up to 16 particles in $`16`$ doubly degenerate levels. For $`d/G=0.5`$ and 16 single-particle levels, we develop two clear minima for the free energy. This is seen in Fig. 1 where we show the free energy as function of excitation energy using Eq. (5) at temperatures $`T=0.5`$, $`T=0.85`$ and $`T=1.0`$ MeV. The first minimum corresponds to the case where we break one pair. The second and third minima correspond to cases where two and three pairs are broken, respectively. When two pairs are broken, corresponding to seniority $`𝒮=4`$, the free energy minimum is made up of contributions from states with $`𝒮=0,2,4`$. These contributions serve to lower the free energy. Similarly, with three pairs broken we see a new free energy minimum which receives contributions from $`𝒮=0,2,4,6`$. At higher excitation energies, population inversion takes place, and our model is no longer realistic.

We note that for $`T=0.5`$ MeV, the minima at lower excitation energies are favored. At $`T=1.0`$ MeV, the higher energy phase (more broken pairs) is favored. We see also, at $`T=0.85`$ MeV, that the free-energy minima where we break two and three pairs equal. Where two minima coexist, we may have an indication of a phase transition. Note however that this is not a phase transition in the ordinary thermodynamical sense. There is no abrupt transition from a purely paired phase to a nonpaired phase. Instead, our system developes several such intermediate steps where different numbers of broken pairs can coexist. At e.g., $`T=0.95`$ MeV, we find again two equal minima. For this case, seniority $`𝒮=6`$ and $`𝒮=8`$ yield two equal minima. This picture repeats itself for higher seniority and higher temperatures.

If we then focus on the second and third minima, i.e., where we break two and three pairs, respectively, the difference $`\mathrm{\Delta }F`$ between the minimum and the maximum of the free energy, can aid us in distinguishing between a first order and a second order phase transition. If $`\mathrm{\Delta }F/N`$ remains constant as $`N`$ increases, we have a second order transition. An increasing $`\mathrm{\Delta }F/N`$ indicates a first order phase transition. In Table I we display $`\mathrm{\Delta }F/N`$ for $`N=10`$, 12, 14 and 16 at $`T=0.85`$ MeV. It is important to note that the features seen in Fig. 1, apply to the cases with $`N=10`$, 12 and 14 as well, where $`T=0.85`$ MeV is the temperature where the second and third minima equal. This means that the temperature where the transition is meant to take place remains stable as function of number of single-particle levels and particles. This is in agreement with the simulations of Lee and Kosterlitz . We find a similar result for the minima developed at $`T=0.95`$ MeV, where both $`𝒮=6`$ and $`𝒮=8`$ coexist. However, due to population inversion, these minima are only seen clearly for $`N=12`$, $`14`$ and $`16`$ particles.

Table I reveals that $`\mathrm{\Delta }F/N`$ is nearly constant, with $`\mathrm{\Delta }F/N0.5`$ MeV, indicating a transition of second order. This result is in agreement with what is expected for an infinite system. It is also easy to see from Fig. 1, that the entropy in the microcanonical ensemble can be convex for certain excitation energy ranges, resulting in eventual negative heat capacities, as inferred from the authors of Refs. . The analysis above however, does not lend support to interpreting this as a sign of a first order phase transition; see also the recent work of Moretto et al. .

We note the important result that for $`d/G>1.5`$, our free energy, for $`N16`$, developes only one minimum for all temperatures. That is, for larger single-particle spacings, there is no sign of a phase transition. This means that there is a critical relation between $`d`$ and $`G`$ for the appearance of a phase transition-like behavior, being a reminiscence of the thermodynamical limit. This agrees also with e.g., the results for ultrasmall metallic grains .

We next compute the distributions of zeros (DOZ) of the canonical partition function and corresponding poles of the specific heat in the complex temperature plane, following a scheme outlined by Lee and Yang , Grossmann et al.  and Borrmann et al. . Following Ref. , we restrict our discussion to the canonical ensemble and denote complex temperatures by $`=\beta +i\tau `$. The partition function of Eq. (4) is now a function of $``$. The authors of Ref.  showed that the distributions of zeros are able to reveal the thermodynamic secrets of small systems in a distinct manner. Major contributions to the specific heat come from zeros close to the real axis, and a zero approaching the real axis infinitely closely causes a divergence in the specific heat. To characterize the DOZ close to the real axis, one assumes that the zeros lie approximately on a straight line. Defining the three parameters $`\tau `$, $`\gamma `$ and $`\alpha `$, see Ref.  for details, one can define entirely the nature of the phase transition. For a true phase transition in the Ehrenfest sense we have $`\tau 0`$. For this case it has been shown  that a phase transition is completely classified by $`\alpha `$ and $`\gamma `$. In the case $`\alpha =0`$ and $`\gamma =0`$ the specific heat $`C_V(\beta )`$ exhibits a $`\delta `$-peak corresponding to a phase transition of first order. For $`0<\alpha <1`$ and $`\gamma =0`$ (or $`\gamma 0`$) the transition is of second order. A higher order transition occurs for $`1<\alpha `$ and arbitrary $`\gamma `$.

In Fig. 2 we show contour plots of the specific heat $`C_v()`$ in the complex temperature plane for $`N=11`$ (a), $`14`$ (b), and $`16`$ (c) particles at normal pairing $`d/G=0.5`$ and the $`N=14`$ (d) in the weak pairing limit, $`d/G=2`$. The poles are at the center of the dark contour regions. We see evidence of two phases in these systems. The first phase, labeled $`I`$ in Fig. 2, is a mixed seniority phase while the second phase, $`II`$, is a paired phase with zero seniority and exists only in even-$`N`$ systems. No paired phase exists in the $`N=11`$ system and no clear boundaries are evident in the weak pairing case. We find that for (b) and (c) the DOZ are apparently distributed along two lines where the intersection occurs at $`\tau _1`$, which is the pole closest to the real axis. As the pairing branch (for $`\beta >\beta _1`$) only encompasses two points, we are unable to precisely determine $`\alpha `$ along this branch while $`\gamma >0`$. Based on our free energy results discussed above, we believe $`\alpha `$ along this branch will be positive. In the mixed phase branch (for $`\beta <\beta _1`$) we find $`\gamma <0`$, and $`\alpha <0`$ in all normal-pairing cases. The parameter $`\tau _1`$, which is a measure of discreteness shows a $`\tau _1N^{1.4\pm 0.12}`$ dependence.

We note several significant results of this work. From two independent methods we find that the transition from the paired seniority zero ground state to a mixed phase state is second order. The free-energy analysis also demonstrates that each transition in seniority phases in the microcanonical ensemble is of second order. The strength of the pairing in these systems determines the nature of the phase transitions. In particular, for a weakly paired system, we found no evidence for two phases, while normal pairing strengths, such as those found in nuclei, may well exhibit the paired-phase and mixed seniority phases that we demonstrated in this model. We will include more realistic interactions to investigate this point in future work. We also found, using Auxiliary Field Monte Carlo computations for this system  together with the histogram method of Refs. , that the energy fluctuations in the canonical ensemble make it rather difficult to extract useful information on the nature of the phase transitions from these techniques.

Oak Ridge National Laboratory is managed by UT-Battelle, LLC for the US Department of Energy under Contract No. DE-AC05-00OR22725. We are indebted to M. Guttormsen and B. Mottelson for many suggestions and comments.