text:
# A non-QE signature of 𝜈_𝑒 appearance in a water Cherenkov detector

## Abstract

We argue that analyzing the $`\nu _e`$-induced CC reaction $`\nu _ene^{}\pi ^0p`$ along with the quasielastic reaction $`\nu _ene^{}p`$ may significantly enhance the sensitivity of a water Cherenkov detector to the subleading oscillation $`\nu _\mu \nu _e`$ at neutrino energies $``$ 2–3 GeV, as projected for the off-axis neutrino beam of NuMI. At the level of standard selections, the multi-ring signal of $`\nu _\mu \nu _e`$ yielded by this $`\pi ^0`$-producing reaction is comparable to the 1-ring (quasielastic) signal in statistical significance. The neutral-current background to $`\nu _ene^{}\pi ^0p`$ can be further suppressed by analyzing spatial separation between the reconstructed primary vertex and the vertices of individual rings.

Detecting the ”subleading” oscillation $`\nu _\mu \nu _e`$ in an off-axis beam with peak energy near 2–3 GeV has emerged as one of the major goals of the NuMI program , addressed at a Fermilab workshop in May 2002  and in a recent Letter of Intent on the subject . A 2.5–3 times bigger baseline than in the proposed JHF2K experiment , that will operate at a lower beam energy of $`E_\nu `$ 0.7–0.8 GeV, will allow to probe the matter effect by comparing the probabilities of the $`\nu _\mu \nu _e`$ and $`\overline{\nu }_\mu \overline{\nu }_e`$ transitions111Without antineutrino running, this can be done by comparing the $`\nu _\mu \nu _e`$ probabilities for a longer (NuMI) and shorter (JHF2K) baselines.. Detector options discussed include a fine-grained low-density calorimeter, a liquid-Argon TPC, and a water Cherenkov spectrometer. Of these, the latter option is based on proven techniques, has an excellent record in neutrino physics , and offers best opportunities in terms of maximum target mass at reasonable cost. It appears however that, in the quasielastic mode only discussed thus far, a water Cherenkov detector will perform worse in NuMI than in JHF2K because of more background to 1-ring electronlike events from single-$`\pi ^0`$ production in NC collisions . It has been estimated  that, due to considerable NC background to quasielastics, a water Cherenkov detector needs to be $`3`$ times as massive as a low-density calorimeter in order to reach similar sensitivity to $`\nu _\mu \nu _e`$ in NuMI conditions. But can the performance of a water Cherenkov detector at $`E_\nu `$ 2–3 GeV be boosted by going beyond quasielastics ?

We believe that at neutrino energies $`2`$ GeV or higher, the sensitivity to $`\nu _\mu \nu _e`$ can be enhanced by also detecting the CC reactions producing a $`\pi ^0`$, $`\nu _ene^{}\pi ^0p`$ and $`\overline{\nu }_epe^+\pi ^0n`$, that largely proceed through excitation of the $`\mathrm{\Delta }(1232)`$ and other baryon resonances222The first observation in a water Cherenkov detector of the corresponding $`\nu _\mu `$-induced reaction, $`\nu _\mu n\mu ^{}\pi ^0p`$, has been reported in .. Compared to $`\nu _ene^{}p`$ and $`\overline{\nu }_epe^+n`$, the cross sections of these reactions333The cross sections of CC and NC reactions are quoted according to NEUGEN . are small for $`E_\nu <1`$ GeV but significant at $`E_\nu `$ 2–3 GeV (see Fig. 1), so that these processes may be relevant to NuMI rather than JHF2K. Depending on whether or not the $`\pi ^0`$ has been fully reconstructed, two different signatures are possible for a water Cherenkov detector:

* Three $`e`$-like rings, of which two fit to $`\pi ^0\gamma \gamma `$;
* Two $`e`$-like rings that would not fit to a $`\pi ^0`$.

Despite a smaller cross section, the $`\pi ^0`$-producing reaction may be competitive with quasielastics because of less neutral-current background: two $`\pi ^0`$ mesons, and not just one, have to be produced in order to mimic the aforementioned 2- and 3-ring signatures. At neutrino energies of a few GeV in particular, the cross section of the NC reaction $`\nu N\nu \pi ^0\pi ^0N`$ should be kinematically suppressed compared to $`\nu N\nu \pi ^0N`$. This simple conjecture is supported by NEUGEN predictions444The uncertainties of these predictions for the cross sections of $`\nu _ene^{}\pi ^0p`$ and $`\nu N\nu \pi ^0\pi ^0N`$ are hard to estimate, as the data are scarce for the former reaction and totally lacking—for the latter reaction., see Fig. 1.

The values of oscillation parameters assumed in the simulation are $`\mathrm{\Delta }m^2`$ = 0.003 eV2, $`\mathrm{sin}^22\mathrm{\Theta }_{23}=1`$, and $`\mathrm{sin}^22\mathrm{\Theta }_{13}=0.1`$ (that is, at the CHOOZ limit ). Matter effects are not accounted for, as here we only wish to compare the oscillation signals in the single-ring and multi-ring channels. Apart from enhancing the matter effect, increasing the baseline shifts the oscillation maximum to higher values of $`E_\nu `$ where the cross sections of the reactions $`\nu _ene^{}\pi ^0p`$ and $`\overline{\nu }_epe^+\pi ^0n`$ are relatively big. Therefore, we select a baseline close to the maximum value for NuMI: $`L=900`$ km. The displacement from the axis of the NuMI medium-energy beam , $`R`$, is varied in the simulation. The $`E_\nu `$ distribution of all $`\nu _\mu `$-induced CC events in the absence of oscillations, illustrated in Fig. 2 for the neutrino mode and $`R=10`$ km, peaks at $`E_\nu 2.6`$ GeV. In the peak region, the intrinsic $`\nu _e`$ component of the beam is some 0.3% of the $`\nu _\mu `$ component (see Fig. 2). Running in the antineutrino mode will yield some 3 times less CC events for the same number of delivered protons.

Taking into account the experimental conditions of a water Cherenkov detector, actually simulated are the ”quasi-inclusive” CC reactions $`\nu _eNe^{}X`$ and $`\nu _eNe^{}\pi ^0X`$ 555Respective antineutrino reactions are implicitly included. and flavor-blind NC reactions $`\nu N\nu \pi ^0X`$ and $`\nu N\nu \pi ^0\pi ^0X`$ in neutrino collisions with water. Here, $`X`$ denotes a system of hadrons other than the $`\pi ^0`$, in which the momenta of all charged particles are below the Cherenkov threshold in water. These reactions are analyzed in terms of visible energy $`E_{\mathrm{vis}}`$, defined as a sum over the energies of all detectable particles: the $`\pi ^0`$ mesons(s) and the charged lepton for CC reactions.

In a water Cherenkov detector, the two photons from $`\pi ^0\gamma \gamma `$ may show up as a single $`e`$-like ring because of a small opening angle (this largely occurs at high $`\pi ^0`$ momenta), or because one of the photons from an ”asymmetric” $`\pi ^0`$ decay is too soft to be detected . The efficiency of $`\pi ^0`$ reconstruction as a function of its momentum will depend on the geometry and instrumentation of a Cherenkov detector; the estimates quoted below are based on the results for the 1-kiloton detector of K2K , as reported in . The momenta of $`\pi ^0`$ mesons emitted in $`\nu _eNe^{}\pi ^0X`$ are plotted in Fig. 3, that also shows the distribution of reconstructed $`\pi ^0`$ mesons (lower histogram). We assume that at least one photon from $`\pi ^0\gamma \gamma `$ is always detected, so that all 1-ring CC events arise from $`\nu _eNe^{}X`$ and all 1-ring NC events—from $`\nu N\nu \pi ^0X`$ with unresolved photon showers. The probability for two photons to form a fake $`\pi ^0`$ candidate is neglected (in SuperK, the r.m.s. width of the $`\pi ^0`$ peak is only $`40`$ MeV ). Depending on whether or not the $`\pi ^0`$ is reconstructible, a CC collision $`\nu _eNe^{}\pi ^0X`$ will produce 3 or 2 rings in the detector. NC events showing 3 (2) rings arise from failing to reconstruct one $`\pi ^0`$ (both $`\pi ^0`$s) in the reaction $`\nu N\nu \pi ^0\pi ^0X`$.

The $`E_{\mathrm{vis}}`$ distributions666Failing to reconstruct a $`\pi ^0`$ will but weakly affect the value of visible energy: in this case, either the two photons from $`\pi ^0\gamma \gamma `$ have merged into a single shower sampled as a whole, or one of them is very soft. of events featuring 1, 2, and 3 $`e`$-like rings are shown in Figs. 46 for incident neutrinos and different values of $`R`$. The three components of the $`E_{\mathrm{vis}}`$ distribution for either channel are: the $`\nu _\mu \nu _e`$ signal (yellow area), the NC background (green area), and the intrinsic-$`\nu _e`$ background (red area). The $`E_{\mathrm{vis}}`$ interval for estimating the effect is selected so as to maximize the ”Figure of Merit” $`S/\sqrt{B}`$, where $`S`$ is the number of signal events and $`B`$ is the total (NC plus intrinsic-CC) background. For either the $`\nu `$ and $`\overline{\nu }`$ settings of the beam, Table 1 compares the 1-ring and multi-ring samples in terms of total $`\nu _\mu \nu _e`$ signals, numbers of signal and background events in the selected $`E_{\mathrm{vis}}`$ windows, and statistical significance. Predictably, the ratio between the multi-ring and 1-ring signals decreases with increasing $`R`$ (or off-axis angle). For incident neutrinos, the multi-ring signal is $``$2 times less than the 1-ring signal in absolute value, but has comparable significance due to less NC background. For incident antineutrinos, the multi-ring signal is substantially less significant than the 1-ring signal.

In a realistic Cherenkov detector, recoil protons often escape detection even for momenta above the Cherenkov threshold . On average, recoil protons have higher momenta in $`\nu _ene^{}\pi ^0p`$ than in $`\nu _ene^{}p`$ due to a broader $`Q^2`$ distribution, so that the multi-ring signal is expected to benefit most from keeping (some) energetic protons. That lifting the upper cut on proton momentum effectively increases the ratio between the multi-ring and 1-ring signals is illustrated by Table 2, to be compared with Table 1.

As indicated in , fast PMT’s and good photocathode coverage may help discriminate between the electron- and $`\pi ^0`$-induced showers by detecting the spatial separation between the conversion points of the two photons from $`\pi ^0\gamma \gamma `$. If shown to be realistic, this will equally apply to 1-ring and multi-ring signatures of $`\nu _\mu \nu _e`$. Yet another geometric handle may be possible for multi-ring topologies only, provided that spatial resolution of the detector is better than photon conversion length $`\lambda _c`$. An important advantage of having more than one ring is that constraining the axes of all rings to a common point in space will yield the position of the primary vertex. Within errors, this should coincide with the reconstructed vertex of a $`e^{}`$-induced shower, whereas the vertex of an unresolved $`\pi ^0`$ shower will be displaced by $`\lambda _c`$ along the shower direction. The spatial resolution of SuperK has been estimated as 18 cm for the vertex of proton decay $`pe^+\pi ^0`$ whose signature is very similar to that of $`\nu _ene^{}\pi ^0p`$, and as 34 cm for the vertex of a single $`e`$-like ring . We have $`\lambda _c40`$ cm for water, so that even a modest improvement in resolution over SuperK will allow to efficiently discriminate between CC and NC multi-ring events and to measure the NC background (this of course needs to be checked by a detailed simulation of detector response).

To conclude, our preliminary results indicate that analyzing the reaction $`\nu _ene^{}\pi ^0p`$ along with the quasielastic reaction $`\nu _ene^{}p`$ may significantly enhance the sensitivity of a water Cherenkov detector to the ”subleading” oscillation $`\nu _\mu \nu _e`$ at neutrino energies $``$ 2–3 GeV. At the level of standard selections, the statistical significance of the multi-ring signal of $`\nu _\mu \nu _e`$ is comparable to that of the 1-ring (quasielastic) signal. The antineutrino reaction $`\overline{\nu }_epe^+\pi ^0n`$ is a less efficient probe of $`\overline{\nu }_\mu \overline{\nu }_e`$ because its cross section is small. The NC background to $`\nu _ene^{}\pi ^0p`$ can be suppressed by reconstructing the vertex of neutrino collision and analyzing spatial separation between the primary and secondary vertices.


format: text

text:
# Unveiling the behavior of the ionized gas in irregular galaxies

## 1. Results and Discussion

We present the results of a detailed kinematical and dynamical analysis performed in two highly interesting gas–rich Irregular galaxies (Irrs): IC 1613 and NGC 4449. The analysis has been accomplished by means of optical Fabry–Perot interferometry mapping the H$`\alpha `$ and \[SII\] lines. Our interest was centered on global and local scales in both galaxies. The results focused on the global radial velocity field in both galaxies show that they are different from each other. IC 1613 displays a highly perturbed velocity field due to local motions caused by superbubbles. NGC 4449 shows a velocity field which is relatively well behaved with a weak decreasing range in radial velocities ($``$60 km s-1) from NE to SW. On local scales, we found that the HII region population for both galaxies displays a velocity dispersion mean value of $``$20 km s-1, which does not reflect the environment of the host galaxy in each case. The HII region population in these galaxies display standard values according to their diameter distribution and luminosity functions (D0= 30–65 pc; $`\alpha `$= 1.5–1.9). The superbubble population in IC 1613 displays expansion velocities (20–40 km s-1), sizes (100–320 pc), dynamical ages (0.7–2 Myr) and mechanical energies (0.6–8$`\times 10^{50}`$ ergs) which can be explained –in the framework of the standard model– in terms of the activity derived from the massive stars (stellar winds and supernova explosions) at those locations. Based on our results (Valdez–Gutiérrez et al. 2001, 2002, Rosado et al. 2001), we conclude that the impressive distribution and kinematical/dynamical behavior of the ionized gas in these two gas–rich Irrs has been moulded –on global and local scales– by the activity of their massive stars; i.e. such activity has “etched” globally and locally the galaxies as they look today. Due to the nature of the present work –and as far as we know–, this is the first time that the kinematics and dynamics (on global and local scales) of the ionized gas in nearby gas–rich Irrs is unveiled in such detail.

## References

Valdez–Gutiérrez, M. et al. 2001, A&A, 366, 35

Valdez–Gutiérrez, M. et al. 2002, AJ, submitted

Rosado, M. et al. 2001, AJ, 122, 194


text:
# Weakly mass-loaded accretion disks

## 1 Introduction

Dynamics of the multi-phase interstellar medium is dominated by mass and momentum exchange between different phases. Ejection of mass into diffuse component from a condensed phase changes velocity of the bulk flow, increases pressure in diffuse phase, and heats diffuse gas due to conversion of kinetic energy of evaporated/destroied clouds into thermal energy of bulk medium (Cowie, McKee & Ostriker 1981, Hartquist et al. 1986). As a result, loaded flows (i.e. the flows dominated by mass input from a condensed phase) can differ qualitatively from the ones in a homogeneous (monophase) medium. In a simplest approach mass loaded flows are described by dynamical equations with source terms. For accretion disks mass loading corresponds physically to a continuing mass infall on the disk surface from the residual gas envelope.

Extremely inhomogeneous structure of the diffuse interstellar gas and molecular clouds shows self-similarity on the scales from 0.02 (the resolution threshold) to 100 pc (see, Falgarone, Phillips & Walker 1991). It seems obvious that such self-similarity can extend to much smaller scales corresponding to protoplanetary condensations in accretion disks. \[Note that the resolution limit, 0.02 pc, corresponds to planetary masses of clumps, $`10^4M_{}`$\]. One can expect that in such conditions clumpy gas from debries of the parent molecular cloud will provide a continuing matter infall on the accretion protoplanetary disk being already formed. In theoretical treatment of protoplanetary disks this circumstance is described by the source terms in mass and momentum equations (see, e.g. Safronov & Vityazev 1983, Papaloizou & Lin 1995). Multiple CI and CO velocity components (Roberge et al. 2000), CaII absorptions (Barnes, Tobin & Pollard 2000) and emission in FeII, SI, SiII, NiII and CI lines (Lecavelier des Etangs et al. 2000), observed in $`\beta `$ Pictoris, are seemingly connected with the infall of clumpy gas on the disk. Similar detections of the infalling circumstellar gas have been reported for Herbig Ae/Be stars (Grinin, Natta & Tambovtseva 1996). One can assume that this phenomenon represents late stages of matter infall from a highly inhomogeneous primordial circumdisk cloud onto the accretion disk. If this is the case, and protoplanetary disks are loaded by infalling and passing through them dense clumps, both vertical and radial structure of such disks, the accretion rate and the temperature distribution can differ considerably from those in homegeneous (monophase) disks. In this paper we analyze properties of steady disks, concentrating mainly on radial temperature distribution, in a simple model of weak mass loading, i.e. when mass loading does not violate the keplerian flow of the bulk motion. In Section 2.1 we describe the model, in Section 2.2 we formulate basic equations describing radial structure of thin accretion disks for weak mass loading limit, in Section 2.3 we describe the solution of energy equation and find radial temperature distribution, and discuss possible observational manifestations, the effects of angular momentum input associated with mass loading on temperature distribution are discussed in Sec. 2.4, Section 3 summarizes the results.

## 2 Model and governing equations

### 2.1 Model of weakly loaded disks

As a model of a mass-loaded disk which admits simple analytical description we will accept the following: an accretion disk is surrounded by a spherical “subsystem” of dense condensations formed of the debries of a parent cloud and orbiting around the central star (Fig. 1). This subsystem can be in a quasistationary state either because the condensations are nearly dissipationless, or an infinite (in practice, very massive) reservoir of gas replenishes dissipation of mass and energy. In the first case, the characteristic life time of the quasistationary orbits determined by interaction of the condensations with the diffuse component, with other condensations and with the disk itself, must be much larger than characteristic rotation period of the disk. For this to fulfill the following relations must hold, respectively

$$N_iN_c,N_cN_c,N_dN_c,$$
(1)
where $`N_i`$, $`N_c`$, and $`N_d`$ are the column densities of the diffuse intercloud gas over the whole circumdisk cloud, of a single clump, and of the disk, respectively, $`N_c=fN_c`$, $`f=𝒩_cR_c^2/R^2`$ is the covering factor of clumps, $`R_c`$ is the clump radius, $`R`$, the radius of the circumdisk cloud. These inequalities definitely fulfill for planetesimals or comet-like bodies. For gaseous clumps it is not so obvious, and one has to assume a large enough density contrast in order that the clumps were long-living and orbiting quasistationary. In the latter case when the clumps are not sufficiently dense, so that for instance $`N_cN_d`$ and they are absorbed by the disk in a single encounter, the overall picture can be kept quasistationary due to replenishment of the clumps from the surrounding molecular cloud as shown schematically on Fig. 1. A gas reservoir with radius $`R310`$ pc can support continuous mass infall during the time $`t13`$ Myr with the rate $`10^710^6M_{}`$ yr-1 if the mean density is only $`n_cf_v0.030.3`$ cm-3, where $`n_c`$ is gas density in clumps, $`f_v`$, the volume filling factor of the clumps; optical depth corresponding to the whole ensemble of such clumps is as small as $`\tau _c10^{21}n_cf_vR0.03`$. As we will see below, mass loading even with one order of magnitude lower mass rate can produce observable effects. We will assume that the quasistationarity takes place in either way. With these assumptions the subsystem of clumps can be treated as stationary. When the clumps intersect the disk they lose a fraction of their mass due to stripping under the action of Rayleigh-Taylor instability (see detailed discussion in Hartquist et al. 1986, Klein, McKee & Colella 1994). The subsystem of clumps is assumed in rest as a whole, i.e. $`𝐯_𝐜(𝐫)=0`$, where the averaging is over the ensemble of clumps at given $`𝐫`$, such that the angular momentum input into the disk from destroying clumps is zero. More general case with $`𝐯_c=𝐯_c(𝐫)`$ will be discussed in Sec. 2.4.

### 2.2 Equations of motions

In this framework the continuity eqation for thin disk is described by the equation

$$\frac{1}{R}\frac{d}{dR}(R\mathrm{\Sigma }u_R)=\mathrm{\Psi },$$
(2)
where $`\mathrm{\Sigma }=_{\mathrm{}}^+\mathrm{}\rho 𝑑z`$, $`u_R=_{\mathrm{}}^+\mathrm{}v_r𝑑z`$, (the interrelation $`\mathrm{\Sigma }u_R=_{\mathrm{}}^+\mathrm{}\rho v_r𝑑z`$ is assumed as well), $`\mathrm{\Psi }=_{\mathrm{}}^+\mathrm{}q𝑑z`$, $`q`$ \[g cm-3 s-1\] is the mass ejection rate from destroying clumps per unit volume. \[In general, $`\mathrm{\Psi }`$ may represent not only mass loading due to destroying cloudletts, but also mass infall onto the disk from the envelope (Safronov & Vityazev 1983), the effects of the mass transfer stream in cataclismic variable stars, and the mass loss through a wind (Papaloizou & Lin 1995)\].

Naviet-Stokes equation for azimuthal component of a stationary disk with zero angular momentum input is written in a standard form (see, Shakura & Sunyaev 1973, Pringle 1981)

$$\frac{d}{dR}(\mathrm{\Sigma }R^3\mathrm{\Omega }u_R)=\frac{d}{dR}\left(\nu \mathrm{\Sigma }R^3\frac{d}{dR}\mathrm{\Omega }\right),$$
(3)
here $`u_\varphi =\mathrm{\Omega }R`$ is explicitly taken into account, $`\mathrm{\Omega }=\mathrm{\Omega }(R)`$ is the local angular velocity, $`\nu `$ is the kinematical viscosity. The change of the angular velocity $`\mathrm{\Omega }`$ in one rotational period $`T_R=2\pi /\mathrm{\Omega }`$ due to mass input is $`|\mathrm{\Delta }\mathrm{\Omega }|\mathrm{\Psi }/\mathrm{\Sigma }`$. From requirement that $`|\mathrm{\Delta }\mathrm{\Omega }|/\mathrm{\Omega }1`$ one obtains restriction on the mass loading $`\mathrm{\Psi }\mathrm{\Sigma }\mathrm{\Omega }/2\pi `$, so that the ratio of the mass input rate $`\dot{M}_p`$ associated with the infalling clumps to the standard mass accretion rate $`\dot{M}_02\pi \nu \mathrm{\Sigma }`$ is

$$\frac{\dot{M}_p}{\dot{M}_0}\frac{1}{4\pi }\frac{R^2\mathrm{\Omega }}{\nu }\frac{1}{4\pi }\frac{u_\varphi }{|u_R|}1.$$
(4)
Obviously, the restrictions on the muss input $`\mathrm{\Psi }`$ are rather weak. Formally, the change in $`\mathrm{\Omega }`$ is larger on a much longer accretion time, $`T_D=R/|u_R|R^2/\nu `$ (here $`u_R\nu /R`$ is explicitly assumed, see Pringle 1981)

$$|\mathrm{\Delta }\mathrm{\Omega }|\mathrm{\Omega }_0\frac{\dot{M}_p}{\dot{M}_0},$$
(5)
where $`\mathrm{\Omega }_0`$ is characteristic angular velocity in the absence of mass loading.

For the radial velocity component we have

$$u_R\frac{du_R}{dR}\frac{u_\varphi ^2}{R}+\frac{1}{\rho }\frac{dp}{dR}+\frac{GM}{R^2}=\frac{\mathrm{\Psi }u_R}{\mathrm{\Sigma }}.$$
(6)
Assuming that for a weak loading the radial component $`u_R`$ is of the order $`\nu /R`$ (see, Pringle 1981) we arrive at

$$u_\varphi ^2=\frac{GM}{R}\left[1+O(^2)+\beta (R)\right],$$
(7)
where $`\beta =\mathrm{\Psi }u_RR^2/GM\mathrm{\Sigma }`$, $`=u_\varphi /c_s`$ is the local azimuthal Mach number, $`c_s`$, the local sound speed. Equations (2)-(6) differ from the standard equations of thin accretion disk (Shakura & Sunyaev 1973, Lyndel-Bell & Pringle 1974, see also reviews by Pringle 1981, and Papaloizou & Lin 1995) by the sources due to mass loading in the r.h.s. of continuity equation (2) and in equation for radial velocity (6). The parameter $`\beta `$ can be estimated as

$$\beta \frac{\alpha }{\pi }\frac{T_R}{T_D}\frac{\dot{M}_p}{\dot{M}_0}1,$$
(8)
with the assumption that $`\dot{M}_p/\dot{M}_0`$ can have low or moderately high values, here $`\alpha `$ is the viscosity parameter (Shakura & Sunyaev 1973), and thus its contribution to the r.h.s. of (7) is negligible, and in this sense the loading will be treated as weak. It is readily seen that approximately $`\beta 2\pi \mathrm{\Psi }/\mathrm{\Omega }\mathrm{\Sigma }`$, and (8) is equivalent to the condition that on one rotational period $`|\mathrm{\Delta }\mathrm{\Omega }|\mathrm{\Omega }`$. This means that the infalling settles onto the generic keplerian motion quickly, and as long as the mass input from clumpuscules is small enough on one rotational period, i.e. (4) fulfills, it can be treated as a weak loading.

Defining the accretion rate as $`\dot{M}=2\pi R\mathrm{\Sigma }u_R`$ (see Pringle 1981) one can obtain the following solutions of equations (2) and (3) for a weakly loaded (i.e. $`\beta 1`$) disk

$$\dot{M}=\dot{M}_{}2\pi \underset{R_{}}{\overset{R}{}}\mathrm{\Psi }R𝑑R,$$
(9)
$$\nu \mathrm{\Sigma }=\frac{\dot{M}_{}}{3\pi }\left[1\left(\frac{R_{}}{R}\right)^{1/2}\right]\frac{2}{3}\underset{R_{}}{\overset{R}{}}\mathrm{\Psi }R𝑑R,$$
(10)
where $`\dot{M}_{}=2\pi (R\mathrm{\Sigma }u_R)_{}`$ is the total mass accretion rate (i.e. the rate on the stellar surface); note, that the mass rate provided by the accretion disk itself in the absence of mass loading is equal

$$\dot{M}_0=\dot{M}_{}2\pi \underset{R_{}}{\overset{\mathrm{}}{}}\mathrm{\Psi }R𝑑R.$$
(11)
The quantity $`\nu \mathrm{\Sigma }`$ determines dissipation at a rate

$$D(R)=\frac{1}{2}\nu \mathrm{\Sigma }(R\mathrm{\Omega }^{})^2=\frac{3GM\dot{M}_{}}{8\pi R^3}[1(R_{}/R)^{1/2}]\frac{3GM}{4R^3}\underset{R_{}}{\overset{R}{}}\mathrm{\Psi }(R)R𝑑R,$$
(12)
and the corresponding luminosity integrated over the disk

$$L_D^\nu =4\pi \underset{R_{}}{\overset{\mathrm{}}{}}D(R)R𝑑R=\frac{1}{2}\frac{GM\dot{M}_{}}{R_{}}3\pi GM\underset{R_{}}{\overset{\mathrm{}}{}}\frac{dR}{R^2}\underset{R_{}}{\overset{R}{}}\mathrm{\Psi }(R^{})R^{}𝑑R^{}.$$
(13)
### 2.3 Energy equation

One of the most interesting and important consequence of the mass loading is connected with possible heating of the disk due to conversion of kinetic energy of the falling gas into heat. The energy equation in the presence of mass loading (see, Cowie, McKee & Ostriker 1981, White & Long 1991), integrated over $`z`$, is written as

$$\mathrm{\Sigma }\dot{\epsilon }=\frac{3}{2}\frac{}{\mu }T\mathrm{\Psi }\left[\frac{1}{2}(\gamma 1)\frac{\rho u^2}{p}\gamma \right]+\frac{1}{2}\nu \mathrm{\Sigma }(R\mathrm{\Omega }^{})^2(T,\rho ),$$
(14)
where $`\epsilon `$ is specific (per unit mass) thermal energy, the first term in the r.h.s. describes the energy input due to inelastic interaction of gas stripped from clumps and involved in keplerian motion of the disk, $`u^2`$ is the mean square relative velocity of clumps and disk (the averaging is over the ensemble of clumps), $``$ is gas constant, $``$ is the radiative cooling rate, we will assume $`=\sigma T^4`$, $`\sigma `$ is Stefan-Bolzmann constant, the second term is the dissipative energy rate (12). In steady case, disk temperature is determined

$$\sigma T^4=\frac{3}{2}\frac{}{\mu }T\mathrm{\Psi }\left[\frac{1}{2}(\gamma 1)\frac{\rho u^2}{p}\gamma \right]+\frac{1}{2}\nu \mathrm{\Sigma }(R\mathrm{\Omega }^{})^2.$$
(15)
The mean square velocity is $`u^2=u_\varphi ^2+u_c^2`$, where $`u_c`$ is the velocity dispersion of clupms encountering the disk, for clumps orbiting around the central star or falling from external molecular cloud $`u_cu_\varphi `$, so we will assume $`u^22u_\varphi ^2`$; the contribution from radial disk velocity $`u_R`$ is neglected. The corresponding contribution to disk luminosity from the first term in the r.h.s. of (14) is

$$L_D^c=6\pi \frac{}{\mu }\underset{R_{}}{\overset{\mathrm{}}{}}T\mathrm{\Psi }\left[\frac{1}{2}(\gamma 1)\frac{\rho u^2}{p}\gamma \right]R𝑑R,$$
(16)
and when compared to the viscous luminosity in the absence of mass loading it is

$$\frac{L_D^c}{L_{D0}^\nu }\frac{R^2\mathrm{\Psi }}{\nu \mathrm{\Sigma }}\frac{\dot{M}_p}{\dot{M}_0}.$$
(17)
In dimensionless units Eq. (15) reads as

$$\theta ^4+6\delta \frac{\mathrm{\Lambda }1}{I}\psi \theta =2(\gamma 1)\frac{\mathrm{\Lambda }1}{I}\frac{\psi }{x}+\frac{\mathrm{\Lambda }}{x^3}\left(1\frac{1}{\sqrt{x}}\right)\frac{\mathrm{\Lambda }1}{I}\frac{1}{x^3}\underset{1}{\overset{x}{}}\psi x𝑑x,$$
(18)
where $`x=R/R_{}`$, $`\mathrm{\Psi }`$ assumed in the form $`\mathrm{\Psi }=\mathrm{\Psi }_0\psi (x)`$, $`\theta =T/T_{}`$, $`\delta =(c_{}/c_M)^2`$

$`T_{}=\left({\displaystyle \frac{3GM\dot{M}_0}{8\pi \sigma R_{}^3}}\right)^{1/4},\mathrm{\Lambda }={\displaystyle \frac{\dot{M}_{}}{\dot{M}_0}},\dot{M}_{}=2\pi \mathrm{\Psi }_0R_{}^2I+\dot{M}_0,`$
$`I={\displaystyle \underset{1}{\overset{\mathrm{}}{}}}\psi x𝑑x,c_{}^2=\gamma {\displaystyle \frac{}{\mu }}T_{},c_M^2={\displaystyle \frac{GM}{R_{}}}.`$ (19)
The problem contains two dimensionless parameters, $`\delta `$, and $`\mathrm{\Lambda }`$. It is readily seen that $`\delta 1`$: for $`R_{}=R_{}`$, $`M=M_{}`$ and $`\dot{M}_010^6M_{}`$ yr-1 this ratio is $`\delta 10^3`$. For weak loading, i.e. $`\dot{M}_{}\stackrel{>}{}\dot{M}_0`$, $`\mathrm{\Lambda }1`$ is also small, however even in this case mass loading can change significantly temperature distribution in the disk and its spectrum. In Fig. 2 we show the temperature distribution $`\theta (x)`$ for the mass-loading function

$$\psi (R)=\left[\frac{RR_{}}{R_L}\right]\mathrm{exp}\left(\frac{RR_{}}{R_L}+1\right),$$
(20)
with the peak $`\psi (R)=1`$ at $`R=R_{}+R_L`$, where $`R_L`$ is the scale of the radial distribution of matter infall, $`R_L=LR_{}`$, $`L=30`$ and a set of $`\mathrm{\Lambda }`$.

It is seen that the most important contribution from additional mass input comes from regions $`x\stackrel{>}{}L`$ where temperature can increase by factor of 2 even for a moderate mass loading $`\mathrm{\Lambda }=2`$. At larger distances, $`x1020L`$, the temperature profiles for $`\mathrm{\Lambda }>1`$ approach the standard profile ($`\mathrm{\Lambda }=1`$), and then go lower because at such distances the heat input from kinetic energy of clumps ($`u^2x^1`$, first term in the r.h.s. of Eq. 18) is smaller than the energy required to heat the infalling (presumably cold) gas (second term in the l.h.s. of Eq. 18). Similar changes are seen for more extended in radius mass input, $`L=100`$ (Fig. 3). In this case, the surface mass input rate for equal $`\mathrm{\Lambda }`$ is $`3`$ times lower than for $`L=30`$, as $`\mathrm{\Lambda }1\mathrm{\Psi }_0L`$, however the increase in total luminosity of loaded disks from the area $`L^2`$ must be larger. For mass input in the form (20) and $`\mathrm{\Lambda }11`$ at $`x<L/3`$ the second term in the r.h.s. of (18) dominates and the first and third terms contribute only around 0.01 % to 1%, at $`x\stackrel{>}{}L/3`$ the situation changes and main contribution (up to 97%) is from the first tirm in (18). At $`x>10L`$ the dominant contribution is again from the second term; the third term at $`x>L/3`$ remains always about 30-50 % of the second. For $`\mathrm{\Lambda }10.1`$ the second term dominates up to $`x2L/3`$ and then at $`x>8L`$, the last term always remains as small as 5-10 % of the second term.

In order to illustrate how the mass loading can change the spectra of accretion disks, we show here the spectral energy distribution (SED) for optically thick disks seen face-on

$$\mathrm{\Sigma }_\nu =2\pi \underset{R_{}}{\overset{R_{out}}{}}B_\nu [T(R)]R𝑑R,$$
(21)
with $`T(R)=T_{}\theta (x)`$, $`B_\nu [T]`$ being the Planck function, $`R_{out}`$, the outer radius of the disk, in the models shown in Figs (2)–(5) $`R_{out}=500`$, contribution from regions $`R>500`$ is small. Defining dimensionless frequency $`y`$ as

$$y=\frac{h\nu }{k_BT_{}},$$
(22)
one can reduce (21) to the form

$$\mathrm{\Sigma }_\nu =\mathrm{\Sigma }_0S(y),$$
(23)
where

$$\mathrm{\Sigma }_0=\frac{16\pi ^2k_B^3T_{}^3R_{}^2}{h^2c^2}.$$
(24)
Fig. 4 shows the spectra from disks with temperature distribution shown in Fig. 2. It is seen that even for relatively small loading, $`\mathrm{\Lambda }=2`$, the spectrum changes substantially: in maximum it increases by factor of 5, the maximum itself is more sharp showing an excess of energy in this range (for $`R=R_{}`$, $`M=M_{}`$ and $`\dot{M}_0=10^6`$ $`M_{}`$ yr-1 the frequency at maximum is $`\nu 11.5\times 10^{13}`$ Hz, or $`\lambda 1020\mu `$m). For a more extended disk ($`L=100`$) the changes in their spectra are more pronounced as shown in Fig. 5: at $`\mathrm{\Lambda }=2`$ in maximum it reaches an order of magnitude increase in comparison with mass-unloaded disks, the maximum itself shifts to a lower frequencies, and the energy excess in this frequency range is more obvious. These two examples show that additional heat input due to mass loading can be responsible (or may contribute) for FIR excess observed from circumstellar disks around young stars (see for a review Natta, Meyer, Beckwith 1997). In order to illustrate such a possibility we compare in Fig. 6 the spectrum expected from an optically thick mass-loaded disk seen face-on with the one observed in the Chameleon cloud (Natta et al. 1997). For the modelled spectrum we assumed the distribution of infalling material of $`L=9.2`$ AU in radius (integration was over $`R=3L`$), its accretion mass rate $`\dot{M}_0=2.2\times 10^7M_{}`$ yr-1, and the mass loading only $`\mathrm{\Lambda }1=0.08`$, the central star is of $`M=1M_{}`$, $`R_{}=R_{}`$. Filled squares show the FIR excess observed from CS Cha and which cannot be attributed to a standard reprocessing accretion disk; open diamonds in optics show presumably the photospheric emission. \[The lack, or negligibly low level, of emission in optical and NIR range from the accretion disk is considered as an indication of a hole in the inner disk (about 0.3 AU in size).\]

For a rough estimate of the effects of mass loading on SEDs from accretion disks, one can use the following approximate expression for temperature in the region where heating due to dissipation of kinetic energy of falling clumps (first term in r.h.s. of Eq. 18) reaches the maximum (at $`x_m\sqrt{L}`$ with $`(s/x)_me/L`$)

$$\theta _m\left[2(\gamma 1)\frac{\mathrm{\Lambda }1}{I}\frac{e}{L}\right]^{1/4},$$
(25)
which correctly by order of magnitude describes the asymptotics at distances $`x`$ where accretion heating (last two terms in Eq. 18) and cooling due to mixing of cold infalling gas (second term in the l.h.s. in Eq. 18) are unimportant – this region corresponds to the plateaus in $`\theta (x)`$ profiles in Figs. 2 and 3. In dimensional units $`T_m0.8(\mathrm{\Lambda }1)^{1/4}L^{1/2}T_{}`$. For the above model it gives $`T_m100`$ K, and the corresponding frequency $`\nu _m6\times 10^{12}`$ Hz.

### 2.4 Effects of angular momentum input

When infalling gas has non-zero angular momentum equation, so that $`v_{c,\varphi }0`$ equation (3) is written in the form (Safronov & Vityazev 1983, Papaloizou & Lin 1995)

$$\frac{d}{dR}(\mathrm{\Sigma }R^3\mathrm{\Omega }u_R)=\frac{d}{dR}\left(\nu \mathrm{\Sigma }R^3\frac{d}{dR}\mathrm{\Omega }\right)+\mathrm{\Psi }(R)R^2v_{c,\varphi }(R),$$
(26)
and the contribution to energy equation is determined by

$$\frac{2\pi }{\dot{M}_0R^2\mathrm{\Omega }}\underset{R_{}}{\overset{R}{}}\mathrm{\Psi }(R)R^2v_{c,\varphi }(R)𝑑R.$$
(27)
In an extreme case when the angular velocity of the infalling gas approaches the keplerian value $`v_{c,\varphi }(R)=\sqrt{GM/R}`$ energy equation reads as

$`\theta ^4+6\delta {\displaystyle \frac{\mathrm{\Lambda }1}{I}}\psi \theta =2(\gamma 1){\displaystyle \frac{\mathrm{\Lambda }1}{I}}{\displaystyle \frac{\psi }{x}}+{\displaystyle \frac{\mathrm{\Lambda }}{x^3}}\left(1{\displaystyle \frac{1}{\sqrt{x}}}\right)`$
$`{\displaystyle \frac{\mathrm{\Lambda }1}{I}}{\displaystyle \frac{1}{x^3}}{\displaystyle \underset{1}{\overset{x}{}}}\psi x𝑑x+{\displaystyle \frac{\mathrm{\Lambda }1}{I}}{\displaystyle \frac{1}{x^{7/2}}}{\displaystyle \underset{1}{\overset{x}{}}}\psi x^{3/2}𝑑x.`$ (28)
It is readily seen that for the mass input rate in the form (20) the contribution from the last term is of the same order of magnitude (about 2/3) as of the third term, and thus never exceeds 30% for $`\mathrm{\Lambda }1`$ and 3–6% for $`\mathrm{\Lambda }11`$.

## 3 Conclusions

Dissipation of kinetic energy of clumpy gas falling onto accretion disks from the debries of a parent protostellar molecular cloud significantly alters the radial temperature distribution. It can be important in outer parts of a reprocessing disk, where temperature decreases sufficiently $`TR^3`$, while the heat input rate associated with mass loading is approximately proportional to the square of the azimuthal velocity $`\dot{E}R^1`$. Thus, mass loading can naturally produce excess of emission in the long-wavelength range of SEDs of the accretion disks, and with an appropriate radial distribution of mass loading can reproduce observed FIR excess. For a one-peak radial distribution of the infall rate localized at $`RL`$, the corresponding temperature in the hot ring for $`R_{}=1R_{}`$, $`M_{}=M_{}`$ is

$$T_m800\left(\frac{\dot{M}_p}{\dot{M}_0}\right)^{1/4}\left(\frac{L}{1\mathrm{AU}}\right)^{1/2}\left(\frac{\dot{M}_0}{10^6M_{}\mathrm{yr}^1}\right)^{1/4}\mathrm{K},$$
(29)
and the wavelength

$$\lambda _m6\left(\frac{L}{1\mathrm{AU}}\right)^{1/2}\left(\frac{\dot{M}_p}{\dot{M}_0}\right)^{1/4}\left(\frac{\dot{M}_0}{10^6M_{}\mathrm{yr}^1}\right)^{1/4}\mu \mathrm{m}.$$
(30)
It is readily seen that the extra luminosity produced by the infalling gas is proportional to the mass loading rate, $`L(\dot{M}_p)\sigma T_m^4L^2\dot{M}_p`$, and for weak mass loading, $`\dot{M}_p\dot{M}_0`$, is small as compared to the total luminosity of an unloaded disk. However, it can result in an observable effect if the loading mass is deposited at large radii, so that the wavelength in the peak of emission (30) falls into FIR range where the reprocessing disk does not contribute singnificantly. This case is shown in Fig. 6 where for mass loading $`\dot{M}_p0.1\dot{M}_0`$ the FIR luminosity of the infalling gas is approximately ten times smaller than the optical luminosity of the reprocessing disk.

The results presented in this paper have only an illustrative meaning, and can be directly applied for explanation of spectral features observed in protoplanetary disks only when clear evidences of mass loading will be found. Nonetheless, they definitely show that mass loading can have an important influence on temperature distribution in accretion disks and their spectra.

I thank S. A. Lamzin, A. V. Lapinov, A. M. Sobolev, I. I. Zinchenko for valuable discussions. This work was supported by RFBR (project No 99-02-16938) and INTAS (project No 1667). This research has made use of NASA’s Astrophysics Data System Abstract Service.

Note added in manuscript 2002 August 8.– After acceptance of this paper I have been informed about the papers by Falcke H. & Melia F., 1997, Astrophys. J., 479, 740, and Coker R. F., Melia F., & Falcke H., 1999, Astrophys. J., 523, 642 which deal with accretion disks around massive black holes loaded by wind infall with application to Sgr A, in the framework of hydrodynamical equations similar to that used in this paper. Their emission spectrum has an excess in the long-wavelength limit due to additional energy input connected with the mass-loading similar to that obtained in our paper.


text:
# On the Temperature Structure of M87

## 1 Introduction

M87 is the giant elliptical galaxy at the center of the Virgo cluster. It is amongst the most studied objects in the sky. Observations at radio and X-ray wavelengths, over the past few years, have allowed us to improve substantially our understanding of the physical processes in this source. The ROSAT analysis of M87 (Böhringer et al. 1995) has shown excess emission correlated with the radio structure and that this emission originates from gas cooler than the surrounding intracluster medium. A detailed radio map of M87 obtained with the VLA at 90 cm (Owen et al. 2000) provided evidence for a very complex radio structure. The authors identified two “bubbles” of synchrotron emission located respectively E and SW some 20-30 kpc from the nucleus. The bubbles are most likely supplied with energy and relativistic plasma which come ultimately from the active nucleus and the inner jet of M87. Largely on the basis of the above X-ray and radio observations Churazov et al. (2001) have developed a model for M87. These authors propose that the bubbles discovered by Owen et al. (2000) buoyantly rise through the cluster atmosphere. During their rise the bubbles uplift relatively cool X-ray emitting gas from the central regions to larger distances. Since the entrained gas is cooler than the ambient gas, the radio bubbles are trailed by cool elongated X-ray features. From the analysis of the XMM-Newton observation of M87, Belsole et al. (2001) have shown that cospatially with the radio bubbles the X-ray emission can be modeled as the sum of two single temperature components, with temperatures of about 2 and 1 keV respectively. In a recent paper based on a Chandra observation of M87, Young et al. (2002) have confirmed the results reported by Belsole et al. (2001).

In this Letter we revisit the XMM-Newton observation, and also make some use of a Chandra observation of M87. We perform a detailed spectroscopic analysis with the aim of better constraining the temperature structure of M87 and consequently gain a deeper understanding of the physical processes at work in this source. The remainder of the paper is organized as follows. In Section 2 we briefly review the XMM-Newton observation and data preparation. In Section 3 we describe the spectral modeling of the M87 data. In Section 4 we focus on the cool component emission. Finally in Section 5 we summarize our main results.

## 2 Observations and Data Preparation

We use XMM-Newton EPIC data from the PV observation of M87/Virgo. Details on the observation, as well as results from the analysis of this object, have already been published in various papers including Böhringer et al. (2001), Belsole et al. (2001), Molendi & Pizzolato (2001), Molendi & Gastaldello (2001), Gastaldello & Molendi (2002) and Matsushita et al (2002).

We have obtained calibrated event files for the MOS1, MOS2 and pn cameras with SASv5.2. In this Letter we use pn data only for the imaging analysis described in Section 4, spectral analysis has been performed on MOS1 and MOS2 data. At the time of writing, the calibration of the MOS detector in the relevant energy range (0.4-4.0 keV) is better than that of the pn detector, moreover the superior spectral resolution of the MOS with respect to pn allows a better modeling of the line rich spectra observed in M87.

Data were manually screened to remove any remaining bright pixels or hot column. There have been no soft-proton flares during the M87 observation and consequently no rejection for time periods with high background has been applied. We have divided the M87 observation in 8 concentric annuli centered on the emission peak. Each ring has been further divided into sectors for a total of 107 regions.

Spectra have been accumulated for MOS1 and MOS2 independently. The Lockman hole observations have been used for the background. Periods in which the background is increased by soft-proton flares have been excluded using an intensity filter; we rejected all events accumulated when the count rate exceeds 15 cts/100s in the $`[1012]`$ keV band. Background spectra have been accumulated from the same detector regions as the source spectra. The vignetting correction has been applied to the spectra rather than to the effective areas, similarly to what has been done by other authors who have analyzed EPIC data (e.g. Arnaud et al. 2001). Spectral fits were performed in the 0.4-4.0 keV band. Data below 0.4 keV were excluded to avoid residual calibration problems in the MOS response matrices at soft energies. Data above 4 keV were excluded because above this energy the spectra show a substantial contamination from hotter gas emitted farther out in the cluster, on the same line of sight.

## 3 Spectral Modeling

All spectral fitting has been performed using version 11.0.1 of the XSPEC package. All models discussed below include a multiplicative component to account for the galactic absorption on the line of sight of M87. The column density is always fixed at a value of 1.8$`\times 10^{20}`$ cm-2, which has been determined from a detailed 21 cm measurement in the direction of M87 (Lieu et al. 1996). The above value is somewhat smaller than the one derived by Stark et al. (1992), 2.5$`\times 10^{20}`$ cm-2, obtained at a lower spatial resolution. We note that such small variations in column density do not affect in any substantial way our results.

We have first compared our data with a single temperature model (vmekal model in XSPEC), which allows to fit separately individual metal abundances. This model has 14 free parameters: the temperature $`T`$, the normalization and the abundance of $`N`$, $`O`$, $`Ne`$, $`Na`$, $`Mg`$, $`Al`$, $`Si`$, $`S`$, $`Ar`$, $`Ca`$, $`Fe`$ and $`Ni`$, which are all expressed in solar units.

As shown in Gastaldello & Molendi (2002) using the apec rather than the mekal model yields very similar results as far as the temperature and the normalization are concerned. Since temperatures and normalizations are the quantities we are interested in we have not run the apec model on our data.

In Figure 1 we show the temperature map obtained from the single temperature model run on the data. The map shows that in general the temperature is smaller at smaller radii. However the most interesting result is the presence of cooler regions which extend from the core in the E and SW directions. We will refer to these regions from now on as the E and SW “arms”. Our findings are in agreement with previous analysis of the XMM M87 observation by Belsole et al. (2001) and Matsushita et al. (2002).

We have then compared our data with a two temperature model (vmekal + vmekal model in XSPEC). This model has two free parameters more than the single temperature model: the temperature and the normalization of the second component. Since we are interested in the temperature structure we have chosen to link the metal abundance of the second component to that of the first component. From tests on sample spectra we have found that, with the exception of Fe, allowing the metal abundances of the two components to vary independently results in little difference in the normalizations and temperatures of the best fitting model. In the case of Fe there is a correlation between the normalization and the Fe abundance for the cooler component, because the spectral range where this component contributes most significantly to the observed spectra is around 0.8-1.0 keV where the cooler component emission is dominated by the Fe-L shell line blend. Tests on sample spectra show that the correlation between normalization and Fe abundance of the cooler component is such that, within the 90% confidence contours, the former can be reduced up to 30-50% if the latter is allowed to increase by a similar amount.

To establish if spectra are better represented by a 2 temperature model than by a 1 temperature model we have made use of the F-test. In Figure 2 we show the value of the probability, associated to the F-test, for the improvement to be statistically significant. Quite interestingly it is only for the two arms extending from the core in the E and SW directions that we have evidence of more than 1 temperature. For all other regions spectra are consistent with being emitted by a single-phase gas.

One possible way to test the robustness of our characterization of the cooler component is to take the spectrum from a region where we have evidence of two temperatures and subtract from this the spectrum of a nearby region at the same angular distance from the center of M87 showing no evidence of the cool component. Since, as detailed later on in this section, the hotter component does not feature strong azimuthal variations, we expect that the spectrum resulting from the above subtraction should be well represented by a single temperature model with temperature and normalization similar to those found for the cool component. We have conducted the above test on 4 sectors, some in the SW and others in the E arm, in all cases we find that the resulting spectra can be well fitted by a single temperature model with temperature and normalization differing no more than 10% from those previously derived from the two temperature fits of the same regions.

We note that the presence of a second temperature component in the E and SW arms was already discussed in the Matsushita et al. (2002) and Belsole et al. (2002) paper; the statistical test we conduct here firmly establishes that this second thermal component is found in the SW and E arms and nowhere else at a comparable statistical significance level.

In general the reduced $`\chi ^2`$ we obtain from our 2 temperature fits indicate that this model provides a fair description of the data. There are three regions for which the probability associated to the $`\chi ^2`$ is smaller than 0.01, visual inspection of the residuals shows that in all three cases an important contribution to the total $`\chi ^2`$ comes from two spectral ranges: the first around 1.1-1.3 keV and the second around 1.7-1.9 keV range. Inspection of other spectra, where the fits are formally acceptable, indicates that residuals in the above energy ranges are rather common. For the 1.1-1.3 keV range the residuals may reflect an incomplete modeling of the Fe-L lines in the mekal code while for the 1.7-1.9 keV they are probably due to an insufficient modeling of the Si edge fine structure in the detector quantum efficiency.

We have fitted spectra showing evidence of a second thermal component with a 3 temperature model (vmekal + vmekal + vmekal model in XSPEC) and with a 2 temperature plus cooling-flow model (vmekal + vmekal + vmcflow in XSPEC). By making use of the F-test we find that the introduction of a third temperature or of a cooling-flow component does not provide a statistically significant improvement to the fit. Thus, from an observational point of view, our results on the temperature structure may be summarized as follows. For most regions the gas can be well described by a single temperature model. For the two arms extending respectively E and SW the spectra can be adequately described by a two temperature model. The lack of improvement in the fit when going to models more complicated than 2 temperatures does not mean that the gas is strictly composed of two phases and there is no emission at all at any other temperature. The correct way of interpreting our result is that the bulk of the emission is coming from gas with a temperature close to either one of the two temperatures found in our best fits. Emission from any other component is relatively small. To quantify the latter statement we have performed the following exercise on the sector with bounding radii 2 and 3 and bounding angles 150o and 180o (angles are counted starting from West and going counter clockwise), which is characterized by a particularly intense cool component. We have added a third temperature component to our two temperature model and frozen the temperatures of the first 2 components as well as the normalization of the hotter component to the best fitting values found from the 2T fit. We have then run fits stepping the temperature of the third component from 0.5 to 5 keV. For each fit we have derived the 90% confidence upper limit for the normalization of the third component. In Figure 3 we show the ratio of the 90% upper limit for the normalization of the third component over the cool component normalization as a function of the temperature of the third component. The plot nicely illustrates that a third component can contribute substantially with respect to the cool component measured in the two temperature fit, only if its temperature is larger than $`0.8`$ keV, emission from gas at temperatures smaller than 0.6 keV has an associated normalization that is no more than 2% of that of the cool component. Some emission from gas with temperatures in the range between 1.7 and 5 keV cannot be ruled out, the intensity of this emission can be relatively large if compared to that of the cool component, however it is much smaller (a few percent) when compared to that of the hot component.

In Figure 4 we report the “normalized” emission integrals (EI), i.e. the emission integrals, as derived by spectral fitting, divided by the area of the region over which they have been measured, as a function of radius (for more details see the Figure caption). For those sectors where we have evidence of two temperatures, we report the normalized EI of the hot and cool components while for all other sectors we plot the normalized EI obtained from the one temperature fit. Similarly in Figure 5 we report the temperatures as a function of radius, showing both hot and cool component temperatures for two temperature regions and the temperature obtained from the 1T fits for all other regions. We note that neither the normalized EIs nor the temperatures for the hot or single phase component show evidence of a large scatter at any given radius, indicating that no strong azimuthal gradients are present. Moreover for those annuli where we have 1 and 2 temperature measurements we find that the temperature and normalized EI of the hot component are not clearly separated from the temperature and normalized EI of the single temperature fits, as is to be expected if the hot component observed in the 2 temperature regions is essentially the same we measure in the 1 temperature regions. Thus the picture emerging from our analysis is that of a 2 temperature structure, with the hot component that is distributed in a regular and symmetric fashion and a cool component that is observed only in the SW and E arms. It is to this latter component that we now turn our attention.

## 4 The cool component

In a number of papers based on ASCA observations various authors (e.g. Ikebe et al. 1996, Makishima et al. 2001) have discussed 2 temperature models as an alternative to the multi-phase cooling-flow model preferred by other authors. In the two temperature model presented by the above authors the higher and lower temperature components are associated respectively to gas in the potential well of the cluster and of the cD galaxy. The potential would thus be characterized by a sort of hierarchical structure. From Figure 2 it is evident that the lower temperature component has a highly asymmetric distribution which is reminiscent of the radio emission (cfr. Owen et al. 2000) and is not associated to the cD which has a much more regular structure. Thus the EPIC observation provides evidence against a model in which the lower temperature component is associated to gas in the potential well of the cD galaxy.

We now wish to further investigate the nature of the cooler component, and more specifically to determine the volume occupied by this component with respect to that occupied by the hot component. For a given region this can be expressed as $`V_cEI_c/n_c^2`$, where $`V_c`$ is the volume occupied by the cool component and $`EI_c`$ and $`n_c`$ are respectively the emission integral and gas density of the cool component. Note that the above estimate is valid under the assumption that the cool component is confined to a scale in projection comparable to the observed extent of the arms within the plane of the sky. To compute $`V_c`$ we need an estimate of $`n_c`$, which can be obtained by assuming that the hot and cool components are in pressure equilibrium with one another.

There are a number of reasons that favor pressure equilibrium between the various components: 1) the X-ray data does not show evidence of shocks anywhere in M87, which would be expected in the presence of substantial pressure jumps and/or supersonic motions; 2) radio observations (Owen et al. 2000) indicate that the radio bubbles are most likely not overpressurized with respect to the X-ray gas; 3) if strong pressure gradients were present between the hotter and cooler X-ray emitting gas and or the radio plasma they would be reduced on a few sound crossing timescales which, assuming a typical size of a few kpc would be smaller than $`10^7`$ yr. Finally we note that, as far as our calculations are concerned, we do not require exact pressure equilibrium. Indeed as we shall see below our conclusions would not change if the pressure of the hot and cool components were within a factor of 2 of each other.

Under the assumption of pressure equilibrium the volume occupied by the cool component can be expressed as $`V_cEI_c/n_h^2kT_c^2/kT_h^2`$, where $`kT_c`$ and $`kT_h`$ are respectively the temperature of the cool and hot components and $`n_h`$ is the gas density for the hot component. The ratio between the volumes occupied by the cool and hot components may then be estimated in two different ways. In the first case we assume that, as for the cool component, $`V_hEI_h/n_h^2`$, where $`V_h`$ and $`EI_h`$ are respectively the volume and the emission integral of the hot component. In this case the ratio between the volumes occupied by the cool and hot components is simply

$$V_c/V_h=EI_c/EI_h\left(\frac{kT_c}{kT_h}\right)^2.$$
$`(1)`$
Alternatively the hot component density can be recovered by deprojecting the M87 data and the volume occupied by the hot component can be estimated by assuming that the linear extent of a given region in projection is similar to that in the plane of the sky, i.e. $`V_h=A^{3/2}`$, where $`A`$ is the area from which the spectrum has been extracted. We have estimated the $`V_c/V_h`$ ratio for all regions where we have evidence of a statistically significant cooler component applying both methods. From direct application of eq. (1) we find that $`V_c/V_h`$ varies between a few times $`10^3`$ and a few times $`10^2`$. For the second method the adopted deprojected $`n_h`$ profile has been derived from Matsushita et al. (2002) and from our own deprojection analysis (Pizzolato et al. in prep.), the two profiles are in good agreement with one another and so are the derived $`V_c/V_h`$ ratios. Estimates of $`V_c/V_h`$ based on the deprojected $`n_h`$ method are about a factor of two larger than those obtained from direct application of eq. (1). The difference between the two estimates is due to the fact that in the first case the emission integral for the hot component, $`EI_h`$, overestimates $`V_hn_h^2`$ because it contains contributions from gas farther out in the cluster observed on the same line of sight. This is borne out by estimates of the relative contribution to $`EI_h`$ from gas in a given spherical shell and from the overlaying gas observed on the same line of sight. Given the approximations we have made to compute $`V_c/V_h`$, and the difference in values returned by our two methods, it is fair to say that our measurements are good within a factor of $``$ 2.

As shown in Figure 4 and already discussed in Section 3, for those annuli where we have 1 and 2 temperature sectors, the normalized EIs for the hot component in 2 temperature sectors are similar to the normalized EIs for the single phase gas in one temperature sectors. This implies that the SW and E arms do not contain substantially less hot phase gas than regions which do not feature strong radio emission and are located at similar radial distances from the cluster core. If, as we have argued above, the radio emitting plasma is in pressure equilibrium with the thermal gas, the above result implies then the filling factor of the radio plasma must be small. Owen et al. (2000), from the analysis of an image of M87 at 90 cm, find that the radio halo is highly filamentary and speculate that the interfilament region could be dominated by the thermal plasma, our findings indicate that this may well be the case.

Our analysis shows that the volume occupied by the cooler phase is, for all regions, much smaller than the volume occupied by the hot phase. Moreover the size of the individual structures, possibly blobs, producing the cool emission, must be considerably smaller than the size of the regions for which we have accumulated spectra. To place tighter constraints on the size of these structures we have produced an emission integral ratio map according to the following procedure. We have produced two images: a “soft” image, S, extracting events from the 0.8-1.0 keV band, where the ratio of the cool component to the hot component is largest and a “hard” image, H, extracting events from the 0.5-0.7 keV and 1.1-4.0 keV bands where, on the contrary, the hot dominates over the cool emission. Assuming a cool component temperature of 0.85 keV and a hot component temperature of 1.7 keV, we have expressed the ratio of the soft component to the hot component in the soft band as a function of the counts in the soft and hard bands. We have then converted this in a cool over hot component emission integral ratio and reported the resulting map in Figure 6. As can be seen, the cool component emission is concentrated in the SW and E arms. Moreover the emission integral of the cool component is always between 0.01 and 0.2 of that of the hot component, roughly corresponding to a $`V_c/V_h`$ ratio in the range 0.0025 to 0.07 if we use eq. (1) and 0.005 to 0.14 if we correct for the factor of 2 found when comparing results obtained from eq. (1) to those obtained using deprojected values for $`n_h`$.

The obvious implication is that, even when pushing the resolution of the EPIC image to its limits and perhaps a little beyond (the pixels in our images are $`10\mathrm{}\mathrm{}\times 10\mathrm{}\mathrm{}`$ ) we still do not resolve the cool component blobs. To further constrain the size of the cool component structures we have produced a new emission integral ratio map using a Chandra observation of M87. The observation was carried out with the S3 chip in the aim point between the 29th and 30th of July 2000, the total effective exposure time, after minor cleaning for flares, is $``$35.3 ks, data was reprocessed using version 2.2.1 of CIAO and 2.10 of the calibration database. Thanks to the better spatial resolution of Chandra and high surface brightness of M87 we can employ $`4\mathrm{}\mathrm{}\times 4\mathrm{}\mathrm{}`$ pixels corresponding to 300 pc $`\times `$ 300 pc. The map has been produced following the same procedure adopted for the EPIC data. The Chandra map, see Figure 7, yields emission integral ratios similar to those found from the EPIC map, indicating that in most regions the size of the cool blobs is smaller than the resolution of the Chandra image. There is however one region 5$`\times `$2 pixels in size, localized $`10\mathrm{}\mathrm{}`$ N of the nucleus, where the ratio has an unusually high value of $`0.7`$. In this region we may be observing a cool component blob with a linear size of $``$ 1 kpc. Our assessments of the typical sizes of blobs are ofcourse rather crude, and are only meant as order of magnitude estimates. More detailed estimates would probably require comparisons with simulations in three dimensions.

The survival of small cool blobs embedded in a hotter medium brings us to the role of thermal conduction. In the absence of any inhibition, the timescale over which the blobs will evaporate from contact with the hotter gas is very short, assuming a typical size of 300 pc and density of $`3\times 10^2`$ cm-3, the conduction timescale is a few times $`10^5`$ yr (note that, even though we consider relatively small sizes, we do not expect conduction to be saturated as the mean free path for electrons is roughly ten times smaller). Thus, either blobs are continuously generated by an unknown mechanism or they survive thanks to substantial suppression of thermal conduction. It is hard to envisage a mechanism capable of continuously generating cool blobs for at least two reasons, firstly because the thermal conduction timescale is shorter than other timescales, such as the sonic or the cooling timescales on which such a mechanism may operate; secondly because, as we shall discuss in more detail later, the entropy of the cool blobs is smaller than that of the hotter gas from which they would presumably be formed. Thus the most likely scenario is one in which thermal conduction is suppressed. Assuming that the age of the blobs is equal or larger than the timescale over which the radio bubbles rise through the cluster atmosphere (a few times 107 yr) we estimate that conduction must be suppressed by a factor of 100 or more. For many years it has been recognized that heat conduction might be suppressed in the ICM (e.g. Binney and Cowie 1981), and many authors (e.g. Chandran et al. 1999) have invoked magnetic fields as a possible means to achieve suppression. It is therefore quite remarkable that we find evidence of suppression only in those regions of M87 where we know magnetic field are present because we see radio halos.

In Figure 8 we show the entropy $`S`$, defined as $`ST/n^{2/3}`$, for the cooler and hotter components as a function of the radius. The density is estimated from a spectral deprojection of the data (Pizzolato et al. in prep.), while the temperature is the emission weighted temperature measured directly from spectral fits, using the density profile presented in Matsushita et al. (2002) produces very similar results. As can be seen, the entropy of the cool phase is virtually everywhere smaller than the entropy of the hot phase. This brings us to the rather important question of how the cool blobs where formed.

In a recent paper Churazov et al. (2001) have developed a detailed model to explain the radio and X-ray structure observed in M87. The authors envisage a scenario in which radio bubbles rise subsonically through the cluster atmosphere. As the bubbles rise, they capture ambient gas, dragging it upwards, during its upward motion the entrained gas is expected to evolve adiabatically. The above model was constructed to satisfy, amongst others, constraints coming from X-ray observations carried out with the ROSAT satellite. It is therefore important to test whether the model can live up to the better quality XMM-Newton EPIC and Chandra ACIS data. The model fails to explain two fundamental observational facts. Firstly, the model predicts that as the blobs rise their temperature drops adiabatically, while from Figure 5 we observe that the blob temperatures are remarkably contained everywhere between 0.8 and 1.0 keV, with no indication of a temperature drop with increasing radius. Secondly, since the cooler gas entropy is virtually everywhere smaller than the hotter gas entropy (see Fig. 8), the cooler gas cannot result from adiabatic expansion of the hotter gas entrained by the radio bubbles. Even assuming a conservative indetermination of a factor 2 on our relative entropy estimates we can allow for an adiabatic evolution only of those blobs currently located at radii larger than $`2^{}`$, blobs at smaller distances from the core would remain unexplained. Note also that, since the timescale on which the bubbles rise, a few times $`10^7`$ years, is smaller than the cooling timescale of the gas, a few times $`10^8`$ years, the entropy of the entrained gas cannot have been reduced through cooling.

The lack of a temperature gradient for the cool gas argues against a long adiabatic rise of the blobs through the cluster atmosphere. As suggested by Nulsen et al. (2001), the entrainment of the cool blobs from the radio bubbles may be unstable, thus the blobs could be captured and then released after having risen only a short way. If conduction is suppressed by magnetic fields associated to the bubble than the release of the blob from the bubble could be rapidly followed by heating of the blob to the ambient gas temperature, this could explain the clear separation between the cool and the hot gas temperatures.

Perhaps the most puzzling issue is the low entropy of the cool blobs. As already pointed out, this argues against the cool blobs originating from the ambient gas. One possible solution is to assume that, as envisaged by Fabian et al. (2001), metals are distributed in a highly inhomogeneous way in the ICM. If, for example, about 1% percent of the gas had a metallicity of 50, in solar units, and the rest of the gas had a metallicity of 1/2, again in solar units, than the metal rich gas would cool at a rate more than 20 times faster than the metal poor gas. Outside the region dominated by radio bubbles thermal conduction would maintain the metal rich blobs at the same temperature of the metal poor gas. In the radio bubbles, where magnetic fields suppress conduction, the metal rich gas would cool on timescales comparable or even shorter than those on which the radio bubbles rise through the cluster atmosphere. In this scenario the cool blobs are metal rich lumps which are kept at a roughly constant temperature due to a balance between radiative cooling and suppressed thermal conduction. Indeed, even if conduction is suppressed, cooling of the blobs increases the gradient between the cool and the hot phase making conduction more efficient. The entropy problem would also be solved because, as the metal rich lumps cool, they radiate away entropy. Unfortunately this rather attractive scenario does not do well when tested against the EPIC data. When fitting some of the spectra of the 2 temperature regions with a 2 temperature model allowing the Fe abundance of the two components to vary independently, we find that the Fe abundance of the cool component cannot be larger than 3-4 in solar units, even when using 99% confidence intervals.

Another alternative to explain the low entropy of the cool blobs is to assume that they have been generated by heating of cooler gas, perhaps through interaction with the radio lobes. However detailed studies at longer wavelengths have shown that the mass of cold gas present in M87 is orders of magnitudes smaller than the mass of the cool phase gas which we observe in X-rays. For instance Sparks et al. (1993) estimate a total mass of the order of $`10^510^7M_{}`$ in the form of ionized gas in dusty filaments, while Edge (2001) derives an upper limit of $`1.5\times 10^8M_{}`$ for any molecular gas, this is to be compared with a total mass of about $`4\times 10^9M_{}`$ for the cool component, which we estimate from our measurements. A possible way out would be to assume that the gas has all been heated up, as might be expected if the core of M87 had recently undergone a starburst phase. If this were to be the case, we would expect to observe an enhancement in the abundance of heavy element. As discussed in Gastaldello and Molendi (2002) and in Canizares et al. (1982) we do indeed have evidence of an O overabundance in the core of M87, however a simple estimate of the mass of the starburst required to explain the O excess yields $`45\times 10^8`$ M, which is an order of magnitude smaller than the mass of the cool gas we observe at X-ray wavelengths.

## 5 Summary

We have performed an analysis of M87 using XMM-Newton EPIC and Chandra ACIS data. Our main findings may be summarized as follows.

* Spectra for most regions of M87 can be adequately fit by single temperature models. It is only for a few regions, which are cospatial with the E and SW radio arms, that we find evidence of a second temperature. Fitting more complicated spectral models to these regions (i.e. including a third temperature or a cooling-flow component) does not improve the quality of the fits.
* The lack of X-ray “holes” at the location of the SW and E radio arms indicates that the radio plasma has a small filling factor. As already suggested by Owen et al. (2000) the radio structures are most likely highly filamented with the interfilamentary regions dominated by the hot phase thermal plasma.
* The cool component fills a small volume compared to the hot component, it is probably structured in blobs with typical linear sizes smaller than a few 100 pc.
* Thermal conduction must be suppressed for these blobs to survive in the hotter ambient gas. Since the cool gas is observed only in those regions of M87 where we have evidence of radio halos our results provide important evidence that magnetic fields could be responsible for the suppression of heat conduction.
* The model proposed by Churazov et al. (2001) is found to be in contradiction with our observational results in two major aspects. Firstly, it predicts that as the blobs rise their temperature drops adiabatically, while we find no indication of a temperature drop with increasing radius. Secondly, the entropy of the cool blobs is in general smaller than that of the hot phase gas, thus cool blobs cannot originate from adiabatic evolution of hot phase gas entrained by buoyant radio bubbles, as suggested by Churazov et al. (2001).
* We have investigated two alternative solutions to solve the puzzle of the low entropy of the cool blobs. The first solution requires that the blobs, which originate from higher entropy ambient gas, rapidly radiate away their excess entropy thanks to their high metal content. However detailed spectral modeling of the EPIC data does not allow for a very high metal content of the cool blobs. The second assumes the cool gas originates from heating of cold gas. However, current estimates of the mass of cold gas are orders of magnitude smaller than the mass of cool gas we observe $`4\times 10^9M_{}`$. Thus, at the present time, we are not capable of providing a convincing origin for the cool gas we observe in M87.

It is a pleasure to acknowledge fruitful discussions with Fabio Pizzolato, Fabio Gastaldello, Simona Ghizzardi, Stefano Ettori, Gianfranco Brunetti and Sabrina De Grandi. S. De Grandi is also thanked for reduction of the Chandra data. An anonymous Referee is thanked for providing helpful comments.