text:
# X-ray luminous galaxies I. Chandra Observations of IRAS00317-2142

## 1 Introduction

Deep ROSAT (Boyle et al. 1995) and Chandra surveys (Mushotzky et al. 2000, Brandt et al. 2000) have discovered large numbers of apparently “normal” galaxies. The optical spectra of these galaxies do not show evidence for either broad or high excitation lines. However, in many cases the high X-ray luminosities observed ($`L_x>10^{42}\mathrm{erg}\mathrm{s}^1`$) suggest the possible presence of an AGN. The nature of these galaxies remains unknown although they are generally believed to be associated with an obscured active nucleus. Unfortunately, no conclusive answers can be given as these are faint and therefore there is a pressing need to discover their counterparts, if any, in the nearby Universe. It is often claimed (Barger et al. 2001) that these sources may resemble the nearby galaxy NGC6240. This galaxy offers one of the few examples in the nearby Universe where although the optical spectrum (Veilleux et al. 1995) gives an unambiguous LINER classification, the X-ray spectrum reveals the presence of a heavily absorbed AGN (Vignati et al. 1999). Another class of nearby galaxies which has many similarities to those discovered in deep ROSAT and Chandra surveys is the ’composite’ class of objects discovered by Veron et al. (1981) and Moran et al. (1996). Their main characteristic is that their \[OIII\] lines are broader than all other narrow lines forbidden or permitted (Moran et al. 1996). The diagnostic emission line diagrams of Veilleux & Osterbrock (1987) classify these as star-forming galaxies. Yet, some of the composites present weak $`H_\alpha `$ broad wings and their X-ray luminosities usually exceed $`10^{42}`$ $`\mathrm{erg}\mathrm{s}^1`$ (Moran et al. 1996) which is considered as the upper limit for the X-ray emission of bona-fide star-forming galaxies (Halpern et al. 1995).

One of the most luminous “composite” objects in the Moran et al. (1996) sample is IRAS00317-2142 with a luminosity of $`10^{43}`$ $`\mathrm{erg}\mathrm{s}^1`$at a redshift of $`z=0.0268`$. This galaxy belongs to the small group of galaxies HCG4. The optical spectrum of IRAS00317-2142 has been discussed in detail by Coziol et al. (1993) and Moran et al. (1996). The diagnostic emission line ratios ($`H_\alpha /[NII]`$ versus $`H_\beta /[OIII]`$) classify it as a star-forming galaxy. Nevertheless, its high X-ray luminosity as well as a faint $`H\alpha `$ wing would classify this as an AGN. Georgantopoulos (2000) presented ASCA hard X-ray observations of this object. The spectrum could be well fitted by a single power-law with $`\mathrm{\Gamma }1.8`$ while the ratio of the X-ray to the broad $`H_\alpha `$ luminosity is $`\mathrm{L}_{(210\mathrm{k}\mathrm{e}\mathrm{V})}/\mathrm{L}_{\mathrm{H}\alpha }2`$. Both the above are typical of Seyfert nuclei and therefore suggest the presence of an unobscured AGN.

In this paper, we present Chandra ACIS-S observations of IRAS00317-2142. The excellent spatial resolution of Chandra combined with its high effective area and good energy resolution can be used to shed more light on the nature of this enigmatic AGN and elucidate the origin of the disagreement between the X-ray and the optical line ratio classification.

## 2 Observations and Data Reduction

IRAS00317-2142 was observed using the Advanced CCD Imaging Spectrometer, ACIS-S, (Garmire et al. 2001) onboard Chandra (Weisskopf et al. 1996). The observation date was 31-08-2001 (Sequence Number 700382). We use the cleaned event file provided by the standard pipeline processing. The resulting exposure time is 19723 sec. The target coordinates are (J2000) $`\alpha =00h34m13.8s`$, $`\delta =21^{}26m21s`$. The whole galaxy falls on only one chip (S3). Charge Transfer Inefficiency (CTI) problems do not affect our observations as S3 is a back-illuminated chip. Each CCD chip subtends an 8.3 arcmin square on the sky while the pixel size is $`0.5`$ arcsec. The spatial resolution on-axis is 0.5 arcsec FWHM. The ACIS-S spectral resolution is $``$100 eV (FWHM) at 1.5 keV. We observed in $`1/2`$-subarray mode in order to minimize pile-up problems. For our source’s count rate the expected pile-up is only $`4\%`$. We used only Grade 0,2,3,4 and 6 events in the analysis. Images, spectra an lightcurves have been created using the CIAO v2.2 software. The imaging and timing analysis were performed using the SHERPA software. For the spectral fitting we use the XSPEC v11 software package. Throughout this paper, we adopt a Hubble constant of $`H_o=65\mathrm{k}\mathrm{m}\mathrm{s}^1\mathrm{Mpc}^1`$.

## 3 Analysis

### 3.1 Imaging

In Figure 1 we show the X-ray contours in the total band (0.3-8 keV) overlayed on a Digital Sky Survey image. The X-ray contours used are those from the adaptively smoothed image created using the CIAO csmooth routine. Most of the X-ray emission emanates from the source coinciding with the optical nucleus. The coordinates of the X-ray nuclear source are $`\alpha =00h34m13.6s`$, $`\delta =21^{}26m18s`$. Its count rate is 0.054 $`\mathrm{cts}\mathrm{s}^1`$ in the 0.3-8 keV band, corresponding to a flux of $`3.2\times 10^{13}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$in the same band. The nucleus contributes the majority of the galaxy counts in the 0.3-8 keV band within 5 arcsec (about 80 %). The unobscured flux is $`3.9\times 10^{13}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$corresponding to a luminosity of $`\mathrm{L}(0.38\mathrm{k}\mathrm{e}\mathrm{V})7\times 10^{41}`$ $`\mathrm{erg}\mathrm{s}^1`$. For the conversion from count rate to flux we used a a power-law spectrum with $`\mathrm{\Gamma }=1.9`$ and a column density of $`N_H=8\times 10^{20}`$ $`\mathrm{cm}^2`$ (see section 3.3).

An X-ray nebulosity can be clearly seen extending out to at least 5 arcsec corresponding to $`6`$ kpc. The extended emission within a circular region of 5 arcsec radius (excluding the nuclear emission i.e. a circular region of 1.5 arcsec radius) has a count rate of 0.014 $`\mathrm{cts}\mathrm{s}^1`$ in the 0.3-8 keV band. The extended X-ray emission in the hard 2-8 keV band is very small (0.001 $`\mathrm{cts}\mathrm{s}^1`$ in the above region). The luminosity of the extended X-ray emission is $`\mathrm{L}(0.38\mathrm{k}\mathrm{e}\mathrm{V})7\times 10^{40}`$ $`\mathrm{erg}\mathrm{s}^1`$, assuming a Raymond-Smith spectrum with a temperature of kT=0.3 keV. We use the SHERPA software to fit the radial profile of the extended emission. We find that a Gaussian with a FWHM=3.5 arcsec fits well the extended emission profile.

Finally, two point sources eastwards of the nucleus are detected at distances of about 10 and 16 kpc respectively. Their coordinates are $`\alpha =00h34m14.1s`$, $`\delta =21^{}26m24s`$ and $`\alpha =00h34m14.4s`$, $`\delta =21^{}26m26s`$ (J2000). Their count rates are $`6.6\times 10^4`$ and $`1.3\times 10^3`$ $`\mathrm{cts}\mathrm{s}^1`$ (0.3-8 keV) corresponding to luminosities of $`L_x6\times 10^{39}`$ $`\mathrm{erg}\mathrm{s}^1`$and $`L_x1\times 10^{40}`$ $`\mathrm{erg}\mathrm{s}^1`$respectively. These ultra-luminous X-ray sources are most likely associated with X-ray binaries (see Fabbiano, Murray & Zezas 2001 and references therein).

### 3.2 Timing

We first examine whether our source presents any long term variability. IRAS00317-2142 has been first observed in 1990 during the ROSAT all-sky survey (exposure 340 s). Subsequently, it has been observed in a ROSAT PSPC pointed observation, between June 22 and 23 1992, with an exposure time of 9.3 ks. The 1-2 keV unabsorbed fluxes obtained in these two ROSAT observations are $`2.0\pm 0.46\times 10^{12}`$ and $`2.1\pm 0.05\times 10^{12}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$implying that the source exhibited no variability in this two year period. The 1-2 keV flux has decreased by a factor of about 25 within a period of 10 years ie between the ROSAT and the Chandra observation.

ASCA observed IRAS00317-2142 between December 11 and 12 1995. The 2-10 keV X-ray flux was $`8\times 10^{13}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. Hence, the 2-10 keV flux has decreased by a factor of about 5 between the ASCA and the Chandra epoch. The long term light curve of IRAS00317-2142 is shown on Figure 2. We plot the unabsorbed flux in the 1-2 keV energy range as this is the common energy band between ROSAT, ASCA and Chandra. The conversion between count rate to flux for all four points has been performed assuming a single power-law spectrum with a photon index $`\mathrm{\Gamma }=1.9`$ absorbed by a column density of $`N_H=8\times 10^{20}\mathrm{cm}^2`$ (the best fit Chandra spectral model).

Although the temporal sampling is sparse the data are consistent with an AGN component that has almost faded away by the time of the Chandra observation. A similar behaviour has been observed in the case of NGC4051 (Guainazzi et al. 1998) and NGC2992 (Gilli et al. 2000). Note that the extraction radius in the case of the Chandra data is much smaller that those used in the ROSAT and ASCA analysis (1 and 2 arcmin radius respectively). However, this cannot explain the observed decrease in flux between the ASCA and Chandra epoch. Indeed, even when we use a 1 arcmin extraction radius we obtain an unabsorbed flux of only $`1\times 10^{13}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$(1-2 keV).

Next, we check for the presence of short term variability in our Chandra data. We extract the light curve using a 1 arcsec radius region. This encompasses $``$90% of the photons from a point source at an energy of 1.5 keV (on-axis). We group our data into 1 ksec time bins. The light curve is shown in Figure 3. The errors correspond to the 68% confidence level. There is no evidence for variability. Indeed, when we fit a constant to the observed count rate we obtain a $`\chi ^2`$ of 23.1 for 20 degrees of freedom (dof). This implies that the no-variability hypothesis cannot be rejected at a statistically significant level (less than 90%).

### 3.3 Spectral Fitting

We group the data so that there are at least 15 counts per energy bin. The quoted errors to the best fitting spectral parameters correspond to the 90 % confidence level for one interesting parameter. We discard data below 0.3 keV and above 10 keV due to both the low response and uncertainties in the calibration matrix. We present spectral fits for the following regions:

a) the nuclear region. Here, we use a 1 arcsec region around the central nuclear pixel. In order to take account of the diffuse background originating from the galaxy itself, we use background from an annulus with an inner and outer radius of 2 and 4 arcsec respectively. We fit the data using a single power-law model absorbed by two columns: one fixed to the Galactic ($`\mathrm{N}_\mathrm{H}=1.5\times 10^{20}`$ Dickey & Lockman 1990) and one intrinsic which is left free to vary. We find that this model provides an excellent fit ($`\chi ^2=47.3/49`$) to the data, see table 1 where all spectral fits results are listed. The power-law has the “canonical” slope, $`\mathrm{\Gamma }=1.91_{0.15}^{+0.17}`$, while the intrinsic absorbing column is $`8\pm 3\times 10^{20}\mathrm{cm}^2`$. Unfortunately, there is not enough signal above 6 keV to provide strong constraints on the Fe line. Actually, the addition of a line at 6.4 keV (rest-frame) results in no change in the $`\chi ^2`$. The 90 % upper limit on the equivalent width of the Fe line is $``$1.4 keV.

We finally note that as it is possible that the above selected background region includes some residual light from the nuclear source, as a cross-check we have also fitted the spectrum using source-free background regions outside the galaxy but still on the S3 chip. The background subtracted spectra give identical fits as the point source overwhelms both the diffuse emission and any small amounts of scattered light from the wings of the PSF. In particular, the intrinsic column density is $`N_H=8_2^{+3}\times 10^{20}\mathrm{cm}^2`$ while the power-law slope is $`\mathrm{\Gamma }=1.94_{0.16}^{+0.17}`$ in excellent agreement with the fits presented above.

b) The “diffuse” emission. Here we use an annulus with inner and outer radius of 2 and 5 arcsec respectively. Again we use source-free regions for the background subtraction. We fit a Raymond-Smith (RS) model to the data absorbed by two column densities, the Galactic and an intrinsic one which is let free to vary. The fit yields $`\chi ^2=12.2/7`$. Some residuals above 2 keV could be either due to a second harder component from X-ray binaries or from residual light from the nuclear point source. We find $`\mathrm{kT}=0.59_{0.42}^{+0.18}\mathrm{keV}`$ for the temperature of the thermal component and $`N_H=1.5_{0.8}^{+2}\times 10^{21}`$ $`\mathrm{cm}^2`$. The abundance is abnormally low $`0.04_{0.03}^{+0.06}`$ possibly suggesting the presence of a multi-temperature medium (see eg Buote 2000).

c) Combined diffuse and nuclear emission. Although we have already analysed the emission from the nucleus and the diffuse emission separately, for the sake of completeness we treat here the total galaxy emission as well. We use a source region of 5 arcsec radius which encompasses both the nuclear and the extended emission. Note however, that this extraction radius does not include the two point sources (ULXs) eatsward of the nucleus. The background is taken from source-free regions on the S3 chip. We use a model consisting of a power-law and a thermal (RS) component. We assume that both are absorbed by the Galactic column density while the power-law component is obscured by an additional intrinsic absorber. In Figure 4 we plot the spectral data, together with the best-fit model and $`\chi ^2`$ residuals. We find $`\mathrm{kT}=0.30_{0.04}^{+0.04}\mathrm{keV}`$, $`\mathrm{\Gamma }=2.01_{0.10}^{+0.10}`$, $`N_H=0.11_{0.02}^{+0.07}\mathrm{cm}^2`$ and an abundance $`1_{0.7}^{+1.5}`$. The fit is good yielding $`\chi ^2=62.9/62`$. When we assume that both the power-law and the thermal component are obscured by the intrinsic column density we obtain a similar fit with $`\chi ^2=62.7/62`$. Therefore, on the basis of the spectral fitting alone it is difficult to decide on the location of the absorbing screen.

## 4 Discussion

The Chandra observations allowed to us to disentangle the nuclear from the extended star-forming component in this peculiar AGN. It becomes obvious that the nucleus contributes the bulk of the X-ray emission, producing over 80% of the photon counts in the total 0.3-8 keV band. The nuclear spectrum is well fit by a single power-law with $`\mathrm{\Gamma }1.9`$ absorbed by an intrinsic column density of $`N_H8\times 10^{20}`$ $`\mathrm{cm}^2`$. Although one has to take into account the use of different spectral models and the degeneracy between the power-law and the column density, this is somewhat higher than the column observed at the ROSAT 1992 observation, $`N_H2\times 10^{20}`$, (Georgantopoulos 2000). Interestingly, the observed column in the Chandra observation is more consistent with the column derived from the Balmer decrement of the narrow line region, $`N_H10^{21}`$ $`\mathrm{cm}^2`$(Georgantopoulos 2000), suggesting that the absorbing screen is located further away from the nucleus.

It is well known that star-forming regions co-exist with AGN (eg Levenson et al. 2001, Della Ceca et al. 2001, Pappa et al. 2002). IRAS00317-2142 is probably the extreme case of this class where the emission from star-forming processes almost completely masks the active nucleus in the optical. The optical spectrum shows a low degree of ionization. The ratio of the optical emission lines ($`H_\alpha /[NII]`$ against $`H_\beta /[OIII]`$) place IRAS00317-2142 in the star-forming galaxy regime. However, the detection of a weak broad emission component to the $`H\alpha `$ line suggests the presence of an AGN. Note that this weak broad emission would hardly be detected at lower fluxes. Actually 3 out of the 7 composite galaxies in Moran et al. (1996) do not present evidence for a broad line. We discuss below a few models for the origin of the X-ray and optical emission in these objects, starting from the most possible one.

A strong star-forming component Moran et al. (1996) favoured a scenario in which the star-forming component is masking the weak nucleus in the optical while in the X-rays it is the nuclear emission which dominates. The Chandra observations have shown that in the X-ray band this is clearly the case. In the optical a strong star-forming component is observed. The integrated galaxy light has $`M_B21`$ and very blue colours with U-B=-0.6 (Coziol et al. 1993). This component is also detected in X-rays having $`L_x7\times 10^{40}`$ $`\mathrm{erg}\mathrm{s}^1`$. The 2 arcsec slit used by Moran et al. (1996) corresponds to a large area on the galaxy ($`2.5kpc`$) and thus the masking of a ’normal’ AGN by the powerful star-forming component offers the most plausible scenario.

An obscured nucleus. The subtle AGN features in the optical spectrum could be due to a weak UV ionizing radiation field. For example circumnuclear obscuration could cause this effect. However, the observed X-ray column ($`10^{21}\mathrm{cm}^2`$) corresponds to only a couple of magnitudes of absorption in the U band. Therefore one would need a more contrived model where the optical is obscured while the X-rays are not (eg a warm absorber model). According to this the X-ray emission emerges relatively unattenuated through an ionised column density while ’pockets’ of dust in this absorber are responsible for the partial extinction of the optical radiation. Still, the X-ray spectrum reveals no absorption edges due to ionized matter. Moreover, in the standard AGN unification model, the Narrow Line Region, above the torus should be ionised by the nuclear radiation and hence we would expect our galaxy to be classified as a Seyfert in the standard line ratio diagnostic diagrams. Only in a geometry where the nucleus is spherically covered by the obscuring screen, we should expect the suppression of high excitation lines. However, an argument against such a scenario comes from the $`\mathrm{L}_{(210\mathrm{k}\mathrm{e}\mathrm{V})}\mathrm{L}_{\mathrm{H}\alpha }`$ relation. Ward et al. (1988) found that in the case of unobscured AGN, the above two luminosities are tightly correlated. In our case, the 2-10 keV X-ray luminosity observed by ASCA (December 1995) and the optical broad $`H\alpha `$ luminosity (obtained one year earlier) follow the Ward et al. (1988) relation (Georgantopoulos 2000). This implies that the optical radiation is not heavily attenuated, provided of course that IRAS00317-2142 follows the same spectral energy distribution as the AGN in the sample studied by Ward et al. (1988). Again the fact that the optical and the X-ray observations were not simultaneous introduces some uncertainty on the actual value of the $`L_x/L_{H\alpha }`$ ratio.

X-ray outbursts Large amplitude X-ray outbursts have been reported in a few galaxies in the last few years. For example the optically quiescent galaxies IC3599 and NGC5905 have exhibited outbursts of a factor of 100 (see Komossa & Bade 1999 and references therein). The peak luminosities exceeded $`L_X=10^{42}`$ $`\mathrm{erg}\mathrm{s}^1`$while the spectra were very soft ($`\mathrm{\Gamma }4`$). These outbursts were interpreted as tidal disruptions of stars by supermassive black holes. NGC5905 and our object bear some similarities, namely the absence of high excitation lines and the very large variability amplitude. On the other hand the “X-ray outburst galaxies” such as NGC5905 show very rapid decay in their flux (Komossa & Bade 1999) while IRAS00317-2142 has remained constant in flux between the two ROSAT observations. Furthermore, our X-ray spectrum is not as steep as these of Komossa & Bade (1999) but is more like that of a typical AGN spectrum, probably ruling out such a model for the composite objects.

Regardless of the nature of the optical emission in composite objects, it is important to understand why the optical spectrum only subtly indicates the presence of the X-ray components. This is probably not uncommon in the Low Luminosity AGN regime. The “composite” galaxies represent only one example of this. Other even more striking examples are those of the passive galaxies detected in deep Chandra and XMM surveys (Mushotzky et al. 2000, Comastri et al. 2002). Most of these objects present no optical emission lines while their X-ray emission rivals that of AGN. Finally, the most well known case of the optical/X-ray controversy classification is probably the Compton-thick Seyfert-2 galaxy NGC6240 which appears to have a LINER optical spectrum. All the above examples cast severe doubt on the use of optical spectroscopy alone, in obtaining the identifications of the sources in deep surveys.

## 5 Summary

We have presented Chandra observations of the enigmatic galaxy IRAS00317-2142. This object belongs to the ’composite’ class of ROSAT all-sky survey sources studied by Moran et al. (1996). These galaxies have optical spectra dominated by star-forming galaxy features while in contrast their powerful X-ray emission suggests the presence of an AGN. In this respect the composite class of galaxies has many similarities to the numerous NLXGs detected in deep ROSAT and Chandra surveys. Chandra helped to shed more light on the X-ray emission mechanisms in these objects. The Chandra observations spatially resolve the nuclear emission from the surrounding extended component. They show that the X-ray emission of IRAS00317-2142 is dominated by a nuclear point source having a luminosity of $`6\times 10^{41}`$ $`\mathrm{erg}\mathrm{s}^1`$. The luminosity of this object has rapidly declined (by a factor of about 25) in a period of about 10 years since the first ROSAT observation. This definitively argues for the presence of an AGN in this object. The nuclear component is absorbed by a column density of $`8\times 10^{20}`$ $`\mathrm{cm}^2`$a factor of a few higher than the Galactic column. Thus the X-ray spectrum is suggestive of a type-1 AGN. Extended X-ray emission has been observed around the nuclear source up to scales of about 6 kpc contributing less than 20 % of the total counts. The spectrum of this diffuse component is soft (described by a RS model with a temperature of kT$``$0.6 keV) consistent with the diffuse gas temperatures encountered in nearby star-forming galaxies (eg Read, Ponman & Strickland 1997).

It remains then unclear why the optical spectrum only subtly gives away the AGN features while the X-ray emission is clearly dominated by the active nucleus. There are several models which can account for such a behaviour. The most plausible scenario is that the star-forming component is very powerful and masks the AGN component in the optical. Our observations have implications for the identifications of galaxies detected in deep ROSAT and Chandra surveys. It appears that at least some fraction of the galaxies with apparently ’normal’ star-forming galaxy optical spectra can harbour an obscured nucleus.

We are grateful to the referee for many comments and corrections. This work has been supported by the NASA grant G01-2120X.


format: text

text:
# Echoes in classical dynamical systems

## 1 Introduction

Echoes arise when through suitable manipulations in a system the dynamics is reversed and a more or less complete recovery of the initial state is achieved. Acoustical echoes arise from reflections of sound at walls, spin echoes from reversals of magnetic fields (Hahn 1950, Carr and Purcell 1954), current echoes through a sequence of suitable electromagnetic pulses (Niggemeier et al 1993) and Loschmidt echoes from a reversal of momenta in a Hamiltonian system (Loschmidt 1876). That echoes can also appear in many particle systems is at first surprising since it seems to be in conflict with the irreversibility implied by the second law of thermodynamics. Closer inspection shows, however, that the recovery of the initial state is not perfect, and studies of the deviations tell a lot about the mechanisms that break reversibility.

Several aspects of echo phenomena in dynamical systems have recently been studied in connection with Loschmidt echoes in quantum systems (triggered by Pastawski et al 1995 and Levstein et al 1998). An initial state $`|0`$ is propagated forward in time with Hamiltonian $`H`$ and then back in time with a slightly different Hamiltonian $`H^{}`$ (as suggested by Peres 1984). The loss of coherence is measured in terms of the fidelity $`0|e^{iH^{}t/\mathrm{}}e^{iHt/\mathrm{}}|0`$. This is the same as calculating the overlap between the states $`|t=e^{iHt/\mathrm{}}|0`$ and $`|t^{}=e^{iHt^{}/\mathrm{}}|0`$, obtained by propagating the same initial state $`|0`$ for the same time $`t`$ but with two different Hamiltonians. The decay of the overlap as a function of time and difference in Hamiltonian and the various time regimes have been the subject of several recent papers, e.g. (Jalabert and Pastawski 2001, Jacquod, Silvestrov and Beenakker 2001, Tomsovic and Cerruti 2002, Benenti and Casati 2002, Wisniacki and Cohen 2001, Prosen and Znidaric 2002, Prosen and Seligman 2002).

The present paper is devoted to Loschmidt echoes in classical dynamical systems. If there is no difference between forward and backward equations of motion, the initial state is recovered perfectly. But what happens if there are small differences or if the system is exposed to noise? And what are the differences between echo experiments in integrable and chaotic systems? These questions will be addressed for Gaussian densities in linearized flows: they provide a convenient and sufficiently general class of densities in which a large variety of dynamical behaviour can be identified.

Besides the obvious connection to the quantum echo experiments, the calculations are of some relevance for two other, directly classical situations: numerical trajectory reversals and reversibilty in advection.

When the equation of motion are reversed in numerical calculations trajectories will typically not return to their starting point. For chaotic systems, the inherent sensitivity to initial conditions suggest an exponentially large deviation. A lack of growth has been used as an indicator for quantum regularity (Casati et al 1986). Obviously, no such problems should arise for perfect reversals and perfect numerical integrators since the solutions to the equations of motion are uniquely specified by the initial conditions (baring singular points in the differential equations). The fact that trajectories do not return to their starting points thus reflects numerical inaccuracies from finite time steps and limited resolution. Using Gaussians to represent a cloud of initial conditions, noise to reflect truncation errors and differences between forward and backward integration algorithm this numerical reversibility experiment can be connected to the problem considered here.

A popular and impressive demonstration of echoes in classical systems is provided by flow reversals in viscous liquids: a blob of dye can be stretched out until it is barely visible but upon reversal of the flow it reforms almost completely! The multimedia fluid mechanics CD (Homsey et al 2001) contains several demonstrations in laminar flows. The connection to chaos comes through experiments on chaotic advection (Aref 1984, 2002, Ottino 1990). For instance, in their experiments on chaotic advection, Chaiken et al (1986) noted that if the dye passed through a chaotic region the recovery was less perfect, but they did not investigate this in detail. Since many of the typical flows can be realized experimentally, also the dynamics of the Gaussians discussed below should be experimentally accessible.

The types of system considered here are classical. They may be exposed to additive white noise and forward and backward motion may differ. The state of the system is characterized by a smooth density in phase space, and the evolution equation is the Fokker-Planck equation with appropriate drift term. The discussion will be limited to densities that are Gaussian in shape, for which a reasonably complete general discussion is possible. Superpositions of such Gaussians can be used to approximate other densities. General expressions in arbitrary dimensions will be given, but their implications are most evident in 2-d conservative flows: this is also the case that is accessible in hydrodynamical systems.

The outline of the paper is as follows. In the next section the dynamics of the center of mass and the variances for Gaussian densities will be discussed. In section 3 this information is applied to the discussion of echoes in integrable systems (section 3a) and chaotic ones (3b). Section 4 contains a discussion of the results and some remarks on experimental tests in 2-d advection systems.

## 2 Gaussian densities

### 2.1 Outline of echo experiments

The calculation of an echo naturally divides into two steps: the forward evolution up to some time $`T`$ under one set of equations, followed by the backward evolution under a perhaps slightly different set of equations for the same time interval. If a mapping of initial conditions under a class of time evolutions can be found, say $`\rho _f=U_1(T)\rho _i`$ for the forward evolution under flow $`𝐮_1`$, then we can write for the backward evolution with flow $`𝐮_2`$ the formal expression $`\rho _b=U_2(T)\rho _f`$, so that the mapping to the echo state $`\rho _e`$ becomes

$$\rho _e=U_2(T)U_1(T)\rho _i.$$
(1)
Thus, on the technical level it suffices to find $`U`$ for the forward dynamics of a sufficiently large class of densities and systems. Note that in the presence of noise or in a dissipative system there is a difference between (i) comparing the initial state with the state obtained by propagating an initial condition over the complete cycle of forward and backward evolution, and (ii) comparing the states obtained by propagating the same initial condition forward in time with two different flows: dissipation and noise break the reversibility that permitted the change in protocol in the quantum case.

The evolution equation for the densities $`\rho (𝐱,t)`$ is the Fokker-Planck equation,

$$\dot{\rho }=(𝐮\rho )+D\mathrm{\Delta }\rho $$
(2)
where $`D`$ is the molecular diffusion constant. Echoes are most easily identified when initial and final density are sufficiently similar, as is the case for strongly localized objects. More complicated initial densities may be approximated by superpositions of localized ones. The dynamics for localized densities splits, in leading order in moments, into two parts, the motion of the center of mass and the changes in shape and size, as measured by the variances. This expansion may be extended to higher order moments of the density, but the equations become too cumbersome to analyze.

### 2.2 Center of mass motion

The center of mass of a localized density follows a classical trajectory $`𝐱_P(t)`$, where

$$\dot{𝐱}_P(t)=𝐮(𝐱_P(t),t).$$
(3)
The notation here is borrowed from hydrodynamic advection (Aref 1984, 2002), where $`𝐮`$ is a velocity field and $`𝐱_P`$ the trajectory of a particle advected by the fluid. In other situations the velocity field $`𝐮(𝐱,t)`$ has to be replaced by the right hand side $`𝐟`$ of an evolution equation $`\dot{𝐱}=𝐟(𝐱,t)`$.

For the backward propagation, where a small modification of the flow field is permitted, the trajectory may differ from the one during forward evolution and we need to estimate the differences between the two. Let $`𝐱_P(t)`$ be a trajectory in a velocity field $`𝐮`$ and $`𝐱_P(t)+𝐪(t)`$ one in the perturbed velocity field $`𝐮+\delta 𝐮`$. To first order in $`𝐪`$ the equation becomes

$$\dot{𝐪}(t)=A(t)𝐪(t)+\delta 𝐮(𝐱_P(t),t),$$
(4)
where $`A`$ is the linearization of the full velocity field $`𝐮+\delta 𝐮`$ at the trajectory $`𝐱_P`$,

$$A_{ij}=\frac{(u_i+\delta u_i)}{x_j}(𝐱_P(t),t).$$
(5)
With the help of the monodromy matrix $`M(t)`$, the solution to

$$\dot{M}=AM$$
(6)
with initial condition $`M(0)=1`$, a formal solution can be given,

$$𝐪(t)=M(t)\left(q(0)+_0^t𝑑\tau M^1(\tau )\delta 𝐮(𝐱_P(\tau ),\tau )\right).$$
(7)
This is the general solution for the displacement in a perturbed velocity field. A discussion of specific examples will be deferred to section 3 below.

### 2.3 Variances

With $`\stackrel{~}{𝐱}=𝐱𝐱_P(t)`$ the coordinates relative to the center of mass, the density can be written as $`\rho (\stackrel{~}{𝐱},t)=\rho (𝐱𝐱_P(t),t)`$ and the Fokker-Planck equation becomes

$$\dot{\rho }(\stackrel{~}{𝐱},t)=(𝐮(𝐱_P+\stackrel{~}{𝐱},t)\rho )+(𝐮(𝐱_P,t))\rho +D\mathrm{\Delta }\rho ,$$
(8)
where now all spatial derivatives are with respect to the relative coordinate $`\stackrel{~}{𝐱}`$. The localization of the densities allows a linearization of the velocity field near the trajectory, i.e.

$$𝐮(𝐱,t)=𝐮(𝐱_P(t),t)+A(t)\stackrel{~}{𝐱}$$
(9)
with the derivative matrix $`A`$, eq. (5). Then

$$\dot{\rho }=(\mathrm{t}rA)\rho ((\stackrel{~}{𝐱}^TA^T))\rho +D\mathrm{\Delta }\rho .$$
(10)
For conservative Hamiltonian systems and incompressible flows, the trace of $`A`$ vanishes; the discussion will henceforth be limited to that case. The next step in the analysis is to note that the equations are second order with a linear position dependence at most, so that a solution in terms of Gaussian densities is possible (Eckhardt 1990). In an $`n`$-dimensional phase space with variance matrix $`\mathrm{\Gamma }(t)`$ they are given by

$$\rho (\stackrel{~}{𝐱},t)=\pi ^{n/2}\left(det\mathrm{\Gamma }\right)^{1/2}e^{\stackrel{~}{𝐱}^T\mathrm{\Gamma }^1\stackrel{~}{𝐱}}$$
(11)
where $`\mathrm{\Gamma }`$ satisfies

$$\dot{\mathrm{\Gamma }}=2D+A\mathrm{\Gamma }+\mathrm{\Gamma }A^T$$
(12)
(note that $`\mathrm{\Gamma }`$ as a kernel for a quadratic form is symmetric). With the help of the monodromy matrix $`M(t)`$ also a closed expression for the variance matrix can be given,

$$\mathrm{\Gamma }(t)=M(t)\left(\mathrm{\Gamma }_i+2D_0^t\left(M(\tau )^TM(\tau )\right)^1𝑑\tau \right)M(t)^T.$$
(13)
This is the central formula for the dynamics of the variances on which the calculation of the various situations can be based. For the specific cases studied in the next section a direct solution of (12) was found to be simpler.

## 3 Echoes

### 3.1 Classical fidelity

The analysis of echoes now proceeds as follows: We start with an initial density $`\rho _i`$ with variance matrix $`\mathrm{\Gamma }_i`$. This density evolves under the influence of a velocity field $`𝐮`$ and additive white noise for some time $`T`$. At the end of this time interval the position of the density is $`𝐱_P(T)`$ and the variance is $`\mathrm{\Gamma }_f(T)`$, as given by (13). Then the field is reversed. If the reversal is perfect, the new velocity field $`𝐮^{}`$ equals $`𝐮`$ and the center of mass returns exactly. If there is a slight deviation the center of mass will move along a trajectory with a small displacement $`𝐪(t)`$ according to (7). In both cases, however, the variance does not return exactly, unless all noise is suppressed, i.e. $`D=0`$. This holds quite generally, the only approximation being the Gaussian shape of the density and linerization of the velocity field near the trajectory.

A general Gaussian echo may then be displaced with its center of mass by $`𝐪`$ and will have a variance matrix $`\mathrm{\Gamma }_e`$. For a comparison between initial and echo density of states the positions and variances can be used directly. But it is also possible to mimick the quantum fidelity expression and to introduce an overlap between classical phase space densities. The main difference is that quantum wave functions are normalized within the $`L^2`$ norm and classical densities are not. The proper definition of the fidelity as the cosine of the angle between the two densities in Hilbert space then is

$`𝒪`$ $`=`$ $`{\displaystyle \frac{𝑑𝐱\rho _i(𝐱)\rho _e(𝐱)}{\left(𝑑𝐱\rho _i(𝐱)^2\right)^{1/2}\left(𝑑𝐱\rho _e(𝐱)^2\right)^{1/2}}}.`$ (14)
Without the normalization the overlap could change even if both densities evolve in the same way. The definition given by Prosen and Znidaric (2002) thus has to be modified. For the case of two Gaussians, an initial density $`\rho _i`$ centered at zero and with variances $`\mathrm{\Gamma }_i`$, and an echo density

$$\rho _e=\pi ^{n/2}\left(det\mathrm{\Gamma }_e\right)^{1/2}e^{(𝐱𝐪)^T\mathrm{\Gamma }_e^1(𝐱𝐪)}$$
(15)
the overlap becomes

$`𝒪`$ $`=`$ $`2^{n/2}\sqrt{{\displaystyle \frac{\sqrt{det\mathrm{\Gamma }_i}\sqrt{det\mathrm{\Gamma }_e}}{det(\mathrm{\Gamma }_i+\mathrm{\Gamma }_e)}}}e^{𝐪^T\left(\mathrm{\Gamma }_i^1+\mathrm{\Gamma }_e^1\right)𝐪}.`$ (16)
The overlap integral indicates two very different kinds of contributions: The prefactor measures the reduction in overlap by spreading of the density. The exponential factor accounts for the rapid drop off in overlap when the centers of the Gaussians are separated; the Gaussian form is clearly connected to the Gaussian tails in the density, and would be different, e.g., for exponential tails in the density.

### 3.2 Shear flows

In laminar flows neighboring trajectories see slightly different velocities and separate linearly in time. When combined with noise a cubic growth of the variance results. Specifically, consider a 2-d shear flow

$$𝐮=\left(\begin{array}{c}\alpha y\\ 0\end{array}\right)$$
(17)
of shear rate $`\alpha `$. The associated monodromy matrix is

$$M(t)=\left(\begin{array}{cc}1& \alpha t\\ 0& 1\end{array}\right).$$
(18)
During forward evolution the variances become

$`\mathrm{\Gamma }_{11}^{(f)}`$ $`=`$ $`\mathrm{\Gamma }_{11}^{(i)}+2\alpha \mathrm{\Gamma }_{12}^{(i)}T+\alpha ^2\mathrm{\Gamma }_{22}^{(i)}T^2+2DT+{\displaystyle \frac{2}{3}}\alpha ^2DT^3`$ (19)
$`\mathrm{\Gamma }_{12}^{(f)}`$ $`=`$ $`\mathrm{\Gamma }_{12}^{(i)}+\alpha \mathrm{\Gamma }_{22}^{(i)}T+\alpha DT^2`$ (20)
$`\mathrm{\Gamma }_{22}^{(f)}`$ $`=`$ $`\mathrm{\Gamma }_{22}^{(i)}+2DT.`$ (21)
The $`T^3`$ contribution to the variances has also been discussed by Rhines and Young (1983) in the context of fluid mixing.

For the reversal we allow for a different shear rate $`\alpha ^{}`$ and some perturbation in the velocity field. If $`𝐪_f`$ is an initial displacement in the trajectory and $`\delta 𝐮=(u_1,u_2)`$ a constant perturbation to the velocity field, the displacement of the echo will be

$`q_{1,e}`$ $`=`$ $`q_{1,f}(u_1\alpha ^{}q_{2,f})t\alpha ^{}u_2t^2/2`$ (22)
$`q_{2,e}`$ $`=`$ $`q_{2,f}\alpha ^{}u_2t\alpha ^{}u_2t^2/2.`$ (23)
The variances become

$`\mathrm{\Gamma }_{11}^{(e)}`$ $`=`$ $`\mathrm{\Gamma }_{11}^{(i)}+2(\alpha \alpha ^{})\mathrm{\Gamma }_{12}^{(i)}T+(\alpha \alpha ^{})^2\mathrm{\Gamma }_{22}^{(i)}T^2`$ (24)
$`+4DT+{\displaystyle \frac{1}{3}}(2\alpha ^2+8\alpha ^26\alpha \alpha ^{})DT^3`$
$`\mathrm{\Gamma }_{12}^{(e)}`$ $`=`$ $`\mathrm{\Gamma }_{12}^{(i)}+(\alpha \alpha ^{})\mathrm{\Gamma }_{22}^{(i)}T+(\alpha 3\alpha ^{})DT^2`$ (25)
$`\mathrm{\Gamma }_{22}^{(e)}`$ $`=`$ $`\mathrm{\Gamma }_{22}^{(i)}+4DT.`$ (26)
This result for the variances may be verified for a few limiting cases: (i) Without shear $`\alpha =\alpha ^{}=0`$ the diagonal elements increase like $`4DT`$ as for regular diffusion over a time interval $`2T`$. (ii) Without diffusion (D=0) the determinant of the matrix does not change. (iii) Without diffusion (D=0) and equal shear in the forward and backward direction ($`\alpha =\alpha ^{}`$) the reversal is perfect and the initial variances are restored. (iv) The parameters $`\alpha ^{}=\alpha `$ correspond to the situation that the backward integration is just a continuation of the forward integration and the expressions (24)-(26) agree with (19)-(21) for an evolution time of $`2T`$. (iv) For equal shear $`\alpha =\alpha ^{}`$ in forward and backward direction the variances are

$`\mathrm{\Gamma }_{11}^{(e)}`$ $`=`$ $`\mathrm{\Gamma }_{11}^{(i)}+4DT+{\displaystyle \frac{4}{3}}\alpha ^2DT^3`$ (27)
$`\mathrm{\Gamma }_{12}^{(e)}`$ $`=`$ $`\mathrm{\Gamma }_{12}^{(i)}2\alpha DT^2`$ (28)
$`\mathrm{\Gamma }_{22}^{(e)}`$ $`=`$ $`\mathrm{\Gamma }_{22}^{(i)}+4DT;`$ (29)
the prefactor of the cubic term in $`\mathrm{\Gamma }_{11}^{(e)}`$ is smaller than would be obtained from (24) for a time $`2T`$, indicating that the reversal of the shear induced broadening is partial and not complete.

The different time regimes in the variances are easily identified. Consider the terms linear in time first: diffusion will be noticable on a time scale $`T_D1/D`$, the differences in the shear rates on a time scale $`T_\delta 1/|\alpha \alpha ^{}|`$. Thus, a large diffusion can swamp the effects from the difference between the two Hamiltonians. For the nonlinear terms, the one with the difference in shear rates appears around $`T_\delta `$, and the one with the cubic term in diffusion near $`1/\alpha `$. In typical applications the latter should be smaller than $`T_\delta `$.

The classical fidelity contains the determinants of the variances. In the absence of diffusion, $`D=0`$, we have $`det\mathrm{\Gamma }_e=det\mathrm{\Gamma }_i`$ and

$$det(\mathrm{\Gamma }_i+\mathrm{\Gamma }_e)=4det\mathrm{\Gamma }_i+(\alpha \alpha ^{})^2T^2\mathrm{\Gamma }_{22}^{(i)}$$
(30)
so that the fidelity has a $`1\text{const}T^2`$ behaviour for short times and a $`1/(|\alpha \alpha ^{}|T)`$ decay for times longer than $`T_\delta `$. With diffusion and equal forward and backward shear, $`\alpha =\alpha ^{}`$, we have

$`det\mathrm{\Gamma }_e`$ $`=`$ $`det\mathrm{\Gamma }_i+4(\mathrm{\Gamma }_{11}^{(i)}+\mathrm{\Gamma }_{22}^{(i)})DT`$ (32)
$`+16D^2T^2+4\mathrm{\Gamma }_{12}^{(i)}\alpha DT^2+{\displaystyle \frac{4}{3}}\mathrm{\Gamma }_{12}^{(i)}\alpha ^2DT^3+{\displaystyle \frac{4}{3}}\alpha ^2D^2T^4`$
and

$`det\left(\mathrm{\Gamma }_i+\mathrm{\Gamma }_e\right)`$ $`=`$ $`4det\mathrm{\Gamma }_i+8(\mathrm{\Gamma }_{11}^{(i)}+\mathrm{\Gamma }_{22}^{(i)})DT`$ (34)
$`+16D^2T^2+8\mathrm{\Gamma }_{12}^{(i)}\alpha DT^2+{\displaystyle \frac{8}{3}}\mathrm{\Gamma }_{12}^{(i)}\alpha ^2DT^3+{\displaystyle \frac{4}{3}}\alpha ^2D^2T^4`$
so that a short-time behaviour $`1\text{const}t`$ and a long-time behaviour of $`1/(\sqrt{\alpha D}T)`$ follow.

The contributions from displacements enter in the exponent and can in principle introduce rapid decays. Consider, e.g., the case $`\alpha =\alpha ^{}`$ and weak diffusion, weak shear and short times, so that $`DT`$ and $`\alpha T`$ are smaller than the initial variances. Then $`\mathrm{\Gamma }^{(e)}\mathrm{\Gamma }^{(i)}`$ and the only time dependence will come from the perturbations $`𝐪_0`$ and $`\delta u`$: the exponent will contain a polynominal of fourth order in time, and this leads to a rapid and faster than exponential decay. If diffusion is added the increase in variance can compensate part of the displacement growth, but for strong diffusion and large times an exponential decay will remain: the square of $`𝐪^2`$ will increase like $`T^4`$ and the variance increases like $`T^3`$, so that the ratio increases linearly, giving an exponential decay.

Thus, without displacement the overlap integral decays algebraically as contained in the prefactor, but with displacement the decrease can be dramatic when the densities are separated by more than their widths.

### 3.3 Chaotic systems

As a model for a chaotic system take a simple hyperbolic motion, $`𝐮=(\lambda x,\lambda y)`$ and assume an initial density aligned with the unstable ($`x`$-) and stable ($`y`$-) direction, i.e. $`\mathrm{\Gamma }_{12}^{(i)}=0`$. After forward evolution we have

$`\mathrm{\Gamma }_{11}^{(f)}`$ $`=`$ $`\mathrm{\Gamma }_{11}^{(i)}e^{2\lambda T}+{\displaystyle \frac{D}{\lambda }}(e^{2\lambda t}1)`$ (35)
$`\mathrm{\Gamma }_{22}^{(f)}`$ $`=`$ $`\mathrm{\Gamma }_{22}^{(i)}e^{2\lambda T}+{\displaystyle \frac{D}{\lambda }}(1e^{2\lambda t}).`$ (36)
Thus, there is an exponential contraction down to the limit set by diffusion. If the backward integration has a slightly different stretching rate $`\lambda ^{}`$, then for the echo

$`\mathrm{\Gamma }_{11}^{(e)}`$ $`=`$ $`\mathrm{\Gamma }_{11}^{(i)}e^{2(\lambda \lambda ^{})T}+{\displaystyle \frac{D}{\lambda }}(e^{2\lambda T}1)e^{2\lambda ^{}T}+{\displaystyle \frac{D}{\lambda ^{}}}(1e^{2\lambda ^{}T})`$ (37)
$`\mathrm{\Gamma }_{22}^{(e)}`$ $`=`$ $`\mathrm{\Gamma }_{22}^{(i)}e^{2(\lambda \lambda ^{})T}+{\displaystyle \frac{D}{\lambda }}(1e^{2\lambda T})e^{2\lambda ^{}T}+{\displaystyle \frac{D}{\lambda ^{}}}(e^{2\lambda ^{}T}1).`$ (38)
As in the previous case we can study various limiting situations, such as $`\lambda ,\lambda ^{}0`$, where linear diffusion results, or $`\lambda =\lambda ^{}`$, where the expressions (35) and (36) for times up to $`2T`$ are recovered. When forward and backward stretching rates are the same, $`\lambda =\lambda ^{}`$, the variances become

$`\mathrm{\Gamma }_{11}^{(e)}`$ $`=`$ $`\mathrm{\Gamma }_{11}^{(i)}+2{\displaystyle \frac{D}{\lambda }}(1e^{2\lambda T})`$ (39)
$`\mathrm{\Gamma }_{22}^{(e)}`$ $`=`$ $`\mathrm{\Gamma }_{22}^{(i)}+2{\displaystyle \frac{D}{\lambda ^{}}}(e^{2\lambda ^{}T}1).`$ (40)
Note that the variance in $`x`$ has hardly changed whereas the one in $`y`$ grows exponentially. The reason is that the $`x`$-variance grows during the forward integration and collapses then during the backward evolution, down to a limit set by the diffusional broadening. The exponential growth and contraction is thus almost perfectly compensated. For the $`y`$-variance we have first the contraction, down to the limit set by diffusion. The expansion during the backward evolution then starts from this finite amplitude, and not from the exponentially small contraction of the deterministic evolution of the initial variance. As a result, the echo is stretched out along the direction that was the stable one during forward evolution. If forward and backward evolution are interchanged, then so is the orientation of the spreading of the density: it will then point in the $`x`$-direction. For this growth to be noticable the time evolution has to be followed for times longer than about $`(\mathrm{ln}(\lambda \mathrm{\Gamma }^{(i)}/D))/(2\lambda )`$.

The behaviour of the classical fidelity for short times is quadratic or linear, for $`D=0`$ and $`D0`$, respectively. On longer times there is an exponential decay, like $`\mathrm{exp}(|\lambda \lambda ^{}|T)`$ without diffusion and like $`\mathrm{exp}(\lambda ^{}T)`$ with diffusion.

The appearance of differences between Lyapunov exponents reflects a relation between forward and backward flow: the stable and unstable manifolds are aligned. A more general difference between forward and backward flow will break this alignment and introduce exponentials in $`\lambda ^{}`$.

Small perturbations in position and in the velocity fields will grow exponentially. Eq. (7) gives

$`q_{1,e}`$ $`=`$ $`(q_{1,f}u_1/\lambda ^{})e^{\lambda ^{}T}+u_1/\lambda ^{}`$ (41)
$`q_{2,e}`$ $`=`$ $`(q_{2,f}+u_2/\lambda ^{})e^{\lambda ^{}T}u_2/\lambda ^{}.`$ (42)
The contributions from the displacement to the decay are weaker than in the linear shear flow, since the stretching of the variances is in the same directions as the separation of trajectories, so that an exponential increase in $`q_2`$ can be compensated by an exponential increase in variances. However, in the absence of diffusion and for the same Lyapunov exponents in forward and backward direction an exponentially growing displacement can lead to a drastic drop off, like $`\mathrm{exp}(\mathrm{exp}\lambda t)`$, in overlap, simply because the Gaussians are shifted relative to each other.

## 4 Final remarks

Already the simple examples in the previous section show a wide range of dynamical behaviour in classical echo experiments. There are two ingredients: the variation in variances and a displacement between initial and echo density. Without displacement the decay in the integrable system is slower than in the chaotic one, and the same applies when there is a displacement in both systems. However, the overlap can drop off faster in an integrable system with displacement compared to a chaotic one without.

The drop off from the displacement is connected with the Gaussian shape of the densities: if the tails fall off more slowly then also the overlap will decay more slowly as a function of displacement. The shape dependence should be stronger in integrable systems than in chaotic ones, since the stretching of variances and shapes goes in parallel with the exponential growth of displacement.

The calculations are based on Gaussian densities and linearizations of the flow fields near trajectories. This becomes questionable if the densities spread out too far, a problem that occurs more likely and more quickly in chaotic systems than in integrable systems. In integrable systems the dangerous terms come from quadratic dependencies of the winding frequencies on action and from curvatures introduced when mapping the tori onto position space. The density will then coil up in whirls. In a 2-d of freedom system this can be expected to happen linearly in time. In chaotic systems the linear approximation is applicable if the stable and unstable manifolds are close to straight lines within the area covered by the density and if the variations in stretching and contraction rate are small. Typically, this will limit the time interval to a logarithmically short one. For longer times the density develops tendrils that follow the manifolds as they wiggle through phase space (for an early discussion of such effects within the quantum maps, see Berry et al, 1979). In both cases the degree to which the Gaussian approximation breaks down depends on details of the nonlinear contributions and has to be considered for specific models.

The transformation to a comoving frame eliminates the center of mass motion and emphasizes the linearized dynamics near the trajectory. The same discussion thus applies to stationary points. Two examples, shear flow or parabolic type and hyperbolic type, have been discussed here, the third class, elliptic type, has periodically oscillating variances.

Two-dimensional flows provide an ideal testing ground for the results presented here. Localized spots of dye can be prepared as initial densities and their motion can be followed in various 2-d flows (Aref 1984, 2002, Ottino 1990, Homsey 2001, Chaiken 1986). Elliptic and hyperbolic points can be realized most easily in cellular flows (Jütner et al 1997, Williams et al 1997, Rothstein et al 1999) Molecular diffusivities are fixed by the selection of dye and solvent, but variations of shear rate through the amplitude of the velocity fields and initial variances through the size of the spot provide enough degrees of freedom to explore the full range in behaviour. In particular the 2-d Lorentz force driven flows (Jütner et al 1997, Williams et al 1997, Rothstein et al 1999) should have enough flexibility to study echoes in flows and to provide quantitative tests of the various expressions derived here. It should be possible to see the cubic laws in the variances in integrable systems, the chaos assisted spreading in chaotic systems and the dependence on the order of forward/backward propagation.

It would also be of interest to reconsider the experimental protocol of Chaiken et al (1986): they followed a closed line of dye. After reversal the line was displaced left and right of the original trace and also had slightly different widths. This suggests different histories of the center of mass of small line elements and passage through regions with different degrees of chaos. Using the complete tracings of the flow field as in (Voth, Haller and Gollub 2002) it should be possible to pin down the different regions and interactions and to characterize the reversal completely.

## Acknowledgment

I would like to thank Th. Seligman and Th. Gorin for stimulating discussions on the relation between classical and quantal fidelity and their symmetries.

## References


text:
# Is Earthquake Triggering Driven by Small Earthquakes?

Now at \]Institute of Geophysics and Planetary Physics,
University of California Los Angeles

## Abstract

Using a catalog of seismicity for Southern California, we measure how the number of triggered earthquakes increases with the earthquake magnitude. The trade-off between this relation and the distribution of earthquake magnitudes controls the relative role of small compared to large earthquakes. We show that seismicity triggering is driven by the smallest earthquakes, which trigger fewer events than larger earthquakes, but which are much more numerous. We propose that the non-trivial scaling of the number of triggered earthquakes emerges from the fractal spatial distribution of seismicity.

Large shallow earthquakes are always followed by aftershocks, that are due the stress change of the mainshock. The number $`n(M)`$ of aftershocks of a mainshock of magnitude $`M`$ has been proposed to scale with $`M`$ as U69 ; KK87 ; K91 ; R85 ; R99 ; SS88 ; O88 ; RJ89 ; YS90 ; DF91 ; MD92 ; HZK00 ; DL01 ; F02

$$n(M)10^{\alpha M}.$$
(1)
This relation accounts for the fact that large earthquakes trigger many more aftershocks than small earthquakes. A similar relation holds for the distribution of earthquake magnitudes $`P(M)`$ GR49  given by

$$P(M)10^{bM},$$
(2)
with $`b`$ typically close to $`1`$, which implies that small earthquakes are much more frequent than large earthquakes.

Because large earthquakes release more energy and trigger more aftershocks than smaller earthquakes, it is usually accepted that interactions between earthquakes and earthquake triggering are dominated by the largest earthquakes. However, because they are much more frequent that larger earthquakes, small earthquakes are also just as important as large earthquakes in redistributing the tectonic forces if $`b=1`$ H92 . Other quantities, such as the Benioff strain $`ϵ10^{0.75M}`$, are dominated by small earthquakes if $`b>0.75`$.

The $`\alpha `$-exponent is an important parameter of earthquake interaction that is used in many stochastic models of seismicity or prediction algorithms KK87 ; K91 ; R85 ; R99 ; O88 ; RJ89 ; CM01 ; F02 . This parameter controls the relative role of small compared to large earthquakes. While there is a significant amount of literature on the $`b`$-value, very few studies have measured accurately the $`\alpha `$ exponent in real seismicity data. Many studies use $`\alpha =b`$ without justification KK87 ; RJ89 ; DF91 ; CM01 ; F02 . In this case, small earthquakes are just as important as larger ones for the triggering process. Using (1) and (2), the global number $`N(M)`$ of aftershocks triggered by all earthquakes of magnitude $`M`$ scales as

$$N(M)=P(M)n(M)10^{(\alpha b)M},$$
(3)
and is indeed independent of $`M`$ in the case $`\alpha =b`$. In the case $`\alpha <b`$, aftershock triggering is controlled by the smallest earthquakes, while the largest earthquakes dominate if $`\alpha >b`$.

A few studies measured directly $`\alpha `$ from aftershocks sequences, using a fit of the total number of aftershocks as a function of the mainshock magnitude SS88 ; YS90 ; MD92 ; DL01 . These studies yield $`\alpha `$-value close to 1, but the limited range of the mainshock magnitude considered and the large scatter of the number of aftershocks per mainshock do not allow an accurate estimation of $`\alpha `$. The case $`\alpha =b`$ also explains another well documented property of aftershocks, known as Bath’s Law B65 ; DL01 ; F02 , which states that the difference between the mainshock magnitude and its largest aftershock is on average close to 1.2, independently of the mainshock magnitude. Again, the limited range of mainshock magnitudes used in these studies and possible biases of data selection V69  does not allow to test the dependence of the magnitude difference as a function of the mainshock magnitude.

Other studies measured $`\alpha `$ indirectly using a stochastic triggering model called “Epidemic Type Aftershock Sequence” model (ETAS) O88 ; K91 ; GO97  based only on the Gutenberg-Richter and Omori laws KK87 ; O88 . This model assumes that each earthquake above a magnitude threshold $`m_0`$ can trigger direct aftershocks, with a rate that increases as $`10^{\alpha M}`$ with its magnitude, and decays with time according to Omori law O94 . The average total number of aftershocks $`n(M)`$, including the cascades of indirect aftershocks, has the same dependence $`10^{\alpha M}`$ with the mainshock magnitude $`M`$ as the number of direct aftershocks. Using this model, $`\alpha `$ can be measured using a maximum likelihood method O88 ; K91 ; GO97 . For instance, GO97  analyzed 34 aftershock sequences in Japan and measured $`\alpha `$ in the range $`[0.21.9]`$ with a mean value of 0.86. The $`\alpha `$-values obtained from the inversion of this model are not well constrained due to the small number of events available and to possible biases of the inversion method. Indeed, these studies do not take into account the influence of earthquakes below the detection threshold, which may significantly bias the estimation of $`\alpha `$. This method may also be biased by the incompleteness of the catalog just after the mainshock, and by possible trade-offs between the ETAS parameters. The regime $`\alpha b`$ of the ETAS model is probably not relevant for real seismicity. If we do not assume a roll-off of the magnitude distribution $`P(M)`$ for large $`M`$, this regime gives a finite time singularity of the seismicity rate which goes to infinity in finite time $`t_c`$ as $`1/(t_ct)^m`$ SH02 . Such a power-law increase of seismic activity can describe the acceleration of the deformation preceding material failure as well as a starquake sequence SH02 , but cannot describe a stationary seismic activity.

In this study, we use a stacking method to estimate the average rate of earthquakes triggered (directly or indirectly) by a previous earthquake as a function of the magnitude of the triggering earthquake. We use the seismicity catalog of Southern California provided by the Southern California Data Center SCEC , which covers the time period 1975-2003, and which is complete above $`M=3`$ for this time period. The magnitude distribution shown in Fig. 1b follows the Gutenberg-Richter law for $`M3`$ and attests for the completeness of the catalog for this time period above magnitude 3. We do not use the usual distinction between “foreshocks”, “mainshocks” and “aftershocks”, and the constraint that the aftershocks must be smaller than the mainshock because this classification is not based on physical differences. Indeed, recent studies have shown that a simple model that assumes that each earthquake can trigger earthquakes of any magnitude, without any distinction between “foreshocks”, “mainshocks” and “aftershocks” can reproduce many properties of real seismicity including realistic foreshock sequences KK87 ; K91 ; O88 ; HSominv ; F02 ; HSG03 . Constraining “triggered earthquakes” to be smaller than the mainshock would obviously underestimate the number of earthquakes triggered by small mainshocks and thus overestimate $`\alpha `$.

We define a “triggered earthquake” as any event occurring in a space-time window $`R\times T`$ after a preceding “mainshock” above the background level, whatever the relative magnitude of the triggered and triggering earthquakes. We consider as a potential “mainshock” each earthquake that has not been preceded by a previous larger earthquake in a space-time window $`d\times T`$ in order to estimate the rate of seismicity triggered by this “mainshock” removing the influence of previous earthquakes.

This definition of “triggered earthquakes” and “mainshocks” contains unavoidably a degree of arbitrariness in the choice of the space-time windows but the estimation of $`\alpha `$ is found to be robust when changing $`T`$, $`R`$ and $`d`$. We have tested different methods for the choice of $`R`$, either fixed or increasing with the mainshock magnitude. We use a distance $`R`$ increasing with the mainshock magnitude because the aftershock zone is usually found to scale with the rupture length U61 ; K02 . We use $`R`$ equal to 1 rupture length of the mainshock. For small mainshock magnitudes, this choice would lead to unacceptable values of $`R`$ smaller than the location error, and thus to underestimate the number of triggered earthquakes of small mainshocks. Therefore, we impose $`R>5`$ km, larger than the location error. Taking $`R`$ fixed has the advantage of not introducing by hand any scaling between the aftershock zone and the mainshock magnitude. However, it may overestimate the number of earthquakes triggered by the smallest mainshocks if $`R`$ is too large, or underestimate the number of triggered events of the largest mainshock if $`R`$ is too small.

The results obtained for $`T=1`$ year, $`R=0.01\times 10^{0.5M}`$ km and $`d=50`$ km are presented in Fig. 1. The rate of triggered earthquakes is found to decay according to Omori’s law $`K(M)/t^p`$, with the same exponent $`p0.9`$ for all mainshock magnitudes $`M`$ (Fig. 1a). The amplitude $`K(M)`$ increases exponentially $`10^{\alpha M}`$ as a function of $`M`$ with $`\alpha =0.81`$ (Fig. 1b). This confirms that the scaling of the rate of triggered earthquakes with $`M`$ follows (1). Our method is more accurate than previous studies SS88 ; YS90 ; MD92 ; DL01 . Indeed, these studies SS88 ; YS90 ; MD92 ; DL01  determine the scaling of $`n(M)`$ with $`M`$ using the total number of aftershocks SS88 ; YS90 ; MD92 ; DL01  in a time window $`[0T]`$ after a mainshock, and can thus be biased by the incompleteness of the catalogs just after a large mainshock or by the background seismicity at large times after a mainshock. In contrast, in order to deal with these problems, we use the seismicity rate in the time window where we observe the Omori law decay characteristic of triggered seismicity.

The value of $`\alpha `$ is robust when increasing or decreasing the distance $`R`$ used for the selection of triggered earthquakes between 1 to 5 rupture lengths, or when increasing the minimum value of $`R`$ from 2 to 10 km. Selecting earthquakes within a disk of fixed radius $`R=50`$ km for all mainshock magnitudes yields a slightly smaller value $`\alpha =0.72`$. Decreasing $`R`$ leads to a smaller value of $`\alpha `$ because it underestimates the number of events triggered by the largest mainshocks, which have a rupture size larger than $`R`$. When increasing $`R`$ from 10 to 100 km, the value of $`\alpha `$ first increases with $`R`$ and then saturates around $`\alpha =0.72`$ for $`R`$ 50 km. We have also checked that $`\alpha `$ is not sensitive to the parameter $`d`$ used for the selection of mainshocks if $`d50`$ km. All values of $`\alpha `$, for reasonable values of the parameters $`R`$ in the range $`30100`$ km, $`T`$ between 0.1 and 2 yrs and $`d>50`$ km, and for different time periods of the catalog, are in the range $`[0.70.9]`$. We have also tested the method on synthetic catalogs generated with the ETAS model. We recover the $`\alpha `$ parameter with an error smaller than 0.05.

For the same catalog of seismicity, we measure using a maximum likelihood method the $`b`$-value of the Gutenberg-Richter law (2) equal to $`b=1.08\pm 0.10`$. We have also tested that the magnitude distribution $`P(M)`$ of triggered events is independent of the mainshock magnitude. Fig. 2 shows the magnitude distribution of triggered events for different ranges of the mainshock magnitude, using the same data as in Fig. 1. This figure shows that a large earthquake can be triggered by a smaller earthquake. Our results suggest that $`\alpha `$ is significantly smaller than the $`b`$ exponent of the magnitude distribution. Whether the exponent $`\alpha `$ varies with region and maybe even with time is an interesting question that is outside the scope of this Letter but we urge further studies in that direction.

We now propose a simple explanation for this non-trivial scaling of the number of triggered earthquakes with the mainshock magnitude, and we suggest that $`\alpha `$ can be related to the fractal structure of the spatial distribution of seismicity. It is widely accepted that the aftershock zone scales with the rupture length U61 ; K02 . While the area affected by the stress variation induced by an earthquake increases with the rupture length, the stress drop is independent of the mainshock magnitude KA75 ; IB01 . The stress variation at a distance from the mainshock proportional to the fault length $`L`$ is thus independent of the mainshock magnitude, neglecting the effect of the finite width of the crust and the visco-elastic deformation in the lower-crust. Therefore, assuming that earthquakes triggered by the stress change induced by the mainshock, the density of earthquakes triggered at a distance up to $`RL`$ from the mainshock is independent of the mainshock magnitude. The increase of the number of triggered events with the mainshock magnitude results only from the increase in the aftershock zone size with the rupture length.

The rupture length is usually related to the magnitude by KA75

$$L10^{0.5M}.$$
(4)
The same relation thus holds between the aftershock zone size $`R`$ and the mainshock magnitude.

In order to estimate the scaling of the number of triggered events with the rupture length, we need to make an assumption about the spatial distribution of earthquakes around the mainshock. Assuming that triggered earthquakes are uniformly distributed on the fault plane, and using (4), the number of earthquakes triggered by a mainshock of magnitude $`M`$ is given by $`n(M)L^210^M`$ YS90  and thus leads to $`\alpha =1`$. The value $`\alpha =0.5`$ obtained for a numerical model of seismicity HZK00  suggests that in this model earthquakes are triggered mostly on the edge of the fracture area of the mainshock HZK00 . Our result $`\alpha =0.8`$ for Southern California seismicity implies that triggered earthquakes are distributed neither uniformly on the rupture plane nor on the edge of the rupture, but rather on a fractal structure of dimension $`D<2`$. Using the definition of the capacity fractal dimension, the number of aftershocks is

$$n(M)R^D,$$
(5)
where $`R`$ is the characteristic length of the aftershock zone. Using (4) and (5), we obtain the scaling of the number of triggered earthquakes with the mainshock magnitude

$$n(M)10^{0.5DM}$$
(6)
which gives $`\alpha =0.5D`$. Our estimation $`\alpha =0.8`$ for Southern California seismicity thus suggests $`D=1.6`$. This value of the fractal dimension of aftershocks hypocenters has never been measured for Southern California seismicity. Our estimate of $`D`$ is significantly smaller than the value measured in the range $`[22.8]`$ for aftershock sequences in Japan GO97 . This fractal dimension of the spatial distribution of triggered earthquakes results in part from the fractal structure of the fault system B01 , but it may also reflect the non-uniformity of the distribution of the earthquakes on the fault due to the heterogeneity of stress or strength on the fault. The fractal dimension of the aftershock distribution may thus be smaller than the fractal distribution of the fault system.

While the energy release and the total slip on faults is controlled by the largest earthquakes, the suggestion that $`\alpha <b`$ implies that small earthquakes may be more important than large earthquakes in triggering earthquakes. We have also checked that the magnitude distribution of triggered earthquakes is independent of the mainshock magnitude (Fig. 2). This implies that earthquake triggering is driven by the smallest earthquakes at all scales, even for the largest earthquakes. Other observations HSominv  support the conclusion that the same mechanisms can explain the triggering of a large earthquake by a smaller one and the triggering of a small earthquake by a previous larger event.

Recent studies F02  have proposed that secondary aftershocks dominate an aftershock sequence, so that subsequent large aftershocks are more likely to be triggered indirectly by a previous aftershock of the mainshock. Our study further suggests that the smallest earthquakes will dominate the triggering of following earthquakes. The importance of small earthquakes casts doubts on the relevance of calculations of direct stress transfer functions to predict seismicity stein , because large earthquakes are likely to be triggered by the smallest earthquakes below the detection threshold of the seismic network. Small earthquakes taken individually have a very low probability of triggering a large earthquake. But because they are much more numerous than larger earthquakes, collectively, they trigger more earthquakes. This result requires the existence of a small magnitude cut-off $`m_0`$, below which earthquakes may occur but cannot trigger earthquakes larger than $`m_0`$, or a change of the scaling of $`N(M)`$ given by (3) for small earthquakes, otherwise the seismicity at all scales would be controlled by infinitely small earthquakes.

I am very grateful to D. Sornette, J.-R. Grasso, A. McGarr, L. Margerin, C. Voisin and T. Gilbert for useful suggestions and discussions. I am also grateful for the earthquake catalog data made available by the Southern California Earthquake Data Center.