text:
# Pulsar bow-shock nebulae

## 1 Introduction

Pulsars are known to be sources of ultrarelativistic magnetized winds (see e.g. Michel & Li michel99 (1999)). If during the supernova explosion the newly formed neutron star acquires a high velocity, it may escape from the supernova remnant before its pulsar activity has ceased. In such a case the pulsar wind interacts with the interstellar medium (thereafter ISM), possibly forming a steady flow, with a bow shock in the direction of the stellar motion. So far only four such nebulae are known: PSR 1957+20 (Kulkarni & Hester kulkarni88 (1988)), PSR 2224+65 (Cordes et al. cordes93 (1993), Chatterjee & Cordes chatterjee02 (2002)), PSR J0437+4715 (Bell bell95 (1995)), and PSR 0740-28 (Jones et al. jones01 (2002)). Probably even the nebula associated with the neutron star RX J185635-3754 (van Kerkwijk & Kulkarni vankerkwijk01 (2001)) may be included in the list, even if there is no agreement whether the nebula is a bow shock or the result of photoionization.

All the nebulae discovered have been detected in Balmer lines (mostly H$`\alpha `$). A pure Balmer spectrum is the signature of a non-radiative shock moving through a partially neutral medium (Chevalier & Raymond chevalier80 (1980), Chatterjee & Cordes chatterjee02 (2002)). When recombination times are long the relevant processes for optical emission are exciting charge-exchange and collisional excitation of neutral H atoms.

As mentioned in a previous paper of this series (Bucciantini bucciantini02 (2002), hereafter Paper II), the ISM is seen by the pulsar as a plane parallel flow with a constant density (at least on the length scale of the bow shock). If the ISM is partially neutral, and if the neutral atoms have mean free paths comparable with, or longer than the size of the bow shock a pure fluid treatment is invalid. The interaction of the relativistic magnetized pulsar wind with the ionized component of the ISM is mediated by the magnetic field advected by the pulsar wind and compressed in the head of the nebula. The neutral atoms interact with the ions (mainly protons) thermalized in the bow shock. The penetration thickness for the H atoms changes from pulsar to pulsar depending on their velocity, their luminosity as well as the density and the neutral fraction of the ISM; and the same holds for the diffusion and ionization (Paper I). Thus to have a more realistic picture of bow-shock nebulae in this regime, it is necessary to include the effect of mass loading from a neutral component as well as the diffusion and ionization of dragged atoms.

As a starting point we have used a two thin layers bow shock model, for which an analytic solution may be derived (Comeron & Kaper comeron98 (1998)). In Paper II we have verified, using a numerical hydrodynamic code, that this analytic model gives reasonable estimates for the values of surface density and tangential velocity in the external layer. These estimates describe the fluid structure in the external layer and have been used to model the interaction and the evolution of the neutral atoms.

Once the density of the neutral H atoms in the external layer has been derived, it is used to evaluate the H$`\alpha `$ luminosity of the nebula, and its variation along the bow-shock. The comparison of the observed luminosity with that obtained from a thin layer model may be used to determine the properties of the local ISM, as well as those of the pulsar wind. In Section 2 we present the equations we have used to derive the distribution of neutral atoms, and the solutions derived integrating these equations; in Section 3 we derive the H$`\alpha `$ luminosity of the nebula and discuss what informations may be obtained from observations.

## 2 The model

### 2.1 The zeroth order two thin layer solution, and the effect of dragging charge-exchange

As said in Paper I, the parameters defining the scenario of a runaway pulsar moving supersonically in the ISM are: the pulsar energy losses, $`L`$; its velocity in the rest frame of the ambient medium, $`v_\mathrm{o}`$; the density respectively of the ionized and neutral component of the ambient medium, $`n_\mathrm{o}`$, $`\mathrm{\Xi }_\mathrm{o}n_\mathrm{o}`$. A supersonic motion implies that the ram pressure of the ISM is greater than its thermal pressure, so that the latter may be neglected in the model. For typical pulsar velocities the related Mach number $`M_\mathrm{o}`$ is $`\stackrel{>}{}4`$ (for a diffuse warm ionized medium).

In Paper II we stressed the fact that a simple bow-shock model (Wilkin wilkin96 (1996)), may not be adequate for the case of a pulsar. The pulsar wind is magnetized; thus the basic assumption of the Wilkin solution (efficient mixing between stellar wind and ambient medium) is not satisfied. For a better description we need to follow separately the two layers, one containing the shocked stellar wind and the other the shocked external medium, allowing two different tangential velocities. In this work we are mainly interested in the behaviour of the external layer of thermalized ISM.

If we retain the thin layer assumption, equations describing the shape, tangential velocity and surface density, can be derived from the momentum and mass conservation. We shall start from the freely crossing case (the neutrals do not interact), then adding on it the effect of mass loading from a neutral component. In fact in typical cases the density of neutral hydrogen may be of the same order as that of protons.

Let us first summarize the freely crossing case. With reference to Fig. 1, let us consider a point (P) on the layer with distance $`r`$ from the star, where the incident fluxes of momentum are $`J_\mathrm{o}`$ and $`J_{}`$. The conservation of perpendicular momentum gives:

$$J_\mathrm{o}\mathrm{sin}^2\alpha _1J_{}\mathrm{sin}^2\alpha _2=\frac{\sigma _{fi}v_f^2+\mathrm{\Phi }_r}{R},$$
(1)
where $`\sigma _{fi}`$ is the surface density of the ionized component in the external layer, $`v_f`$ is the tangential velocity in the external layer, $`R`$ is the curvature radius of the front, and $`\mathrm{\Phi }_r`$ is the momentum flux in the internal layer. For the tangential component of the fluxes in the external and internal layer we have respectively:

$$J_\mathrm{o}\mathrm{sin}\alpha _1\mathrm{cos}\alpha _1=\frac{d(\sigma _{fi}v_f^2)}{ds},$$
(2)
$$J_{}\mathrm{sin}\alpha _2\mathrm{cos}\alpha _2=\frac{d\mathrm{\Phi }_r}{ds},$$
(3)
where $`ds`$ is the line element. Finally the mass conservation equation for the external layer, where surface density is required for the evaluation of penetration thickness, reads:

$$\sigma _{fi}v_f=\frac{J_\mathrm{o}}{2v_\mathrm{o}}y.$$
(4)
These equations may be written in a simpler form using dimensionless variables. The lengths $`y`$ and $`x`$ scale with the distance to the stagnation point (Cordes cordes96 (1996)):

$$d_\mathrm{o}=\sqrt{J_{}/J_\mathrm{o}}=\sqrt{L/4\pi cn_\mathrm{o}v_\mathrm{o}^2},$$
(5)
while the velocity and surface density in the external layer are respectively scaled with $`v_\mathrm{o}`$ and $`n_\mathrm{o}d_\mathrm{o}`$, and the momentum flux in the internal layer with $`J_\mathrm{o}d_\mathrm{o}`$. Equations 1 to 4 may be written:

$`(\sigma _{fi}v_f^2+\mathrm{\Phi }_r){\displaystyle \frac{dz}{dy}}n+{\displaystyle \frac{n(yzx)^2}{r^4}}=0,`$
$`{\displaystyle \frac{d(y\sigma _{fi}v_f^2)}{ydy}}{\displaystyle \frac{z}{n}}=0,`$ (6)
$`{\displaystyle \frac{d(y\mathrm{\Phi }_r)}{ydy}}{\displaystyle \frac{(xz+y)(yzx)}{nr^4}}=0,`$
$`2\sigma _{fi}v_fy=0.`$
Where $`z=dx/dy`$ and $`n=(1+z^2)^{1/2}`$. This set of equations can be integrated numerically, starting from the stagnation point, and using a Taylor expansion near the apex. It is remarkable to see that these equations do not contain any free parameter, so there is only one solution in these dimensionless coordinates (solid lines of Figs. 2 and 3).

Once the freely crossing solution has been obtained from the integration of Eqs. 6, it is possible to introduce the effects of a neutral component of ISM (that we assume to be mainly composed by neutral hydrogen). As said in Paper I, for typical pulsar velocities the main process responsible for the dragging of neutral H in the external layer is charge-exchange. Consequences of the presence of a neutral component are: 1- a new source of perpendicular and tangential momentum, 2- the need for a further mass conservation equation for the neutrals. We assume that, once a neutral atom has been captured in the external layer, it can neither escape via diffusion nor be ionized (an assumption whose implications will be analyzed in section 2.2), and that there is no shear between the neutral and ionized material. The penetration thickness of the external layer is not constant but tends to increase moving away from the stagnation point (Paper II). The fraction of neutrals stopped in this layer is:

$$F=1\mathrm{exp}(\sigma _{cx}(v_{rel})\sigma _{fi}v_{rel}/v_\mathrm{o}\mathrm{sin}\alpha _1),$$
(7)
where $`\sigma _{cx}`$ is the cross section of charge-exchange, $`v_{rel}`$ is the relative velocity between the incoming H atoms and the hot shocked plasma and we have assumed a constant state (i.e. no variations of density, temperature, velocity) across the external layer. We want to stress that the constant state assumption is valid only near the head of the bow shock where the geometrical thickness of the layer is small compared to the size of the nebula. In order to evaluate the relative velocity we assume that all the ionized material in the external layer moves tangentially with a speed $`v_f`$: this assumption is good (Paper I) if the penetration thickness is not too high, otherwise it would be better to use the post shock velocity (according to the Rankine-Hougoniot jump conditions). However the use of the post shock velocity is justified only in case $`F1`$, anyway giving results not too different from the other case. The set of Eqs. 6, becomes:

$`((\sigma _{fi}+\sigma _{fn})v_f^2+\mathrm{\Phi }_r){\displaystyle \frac{dz}{dy}}n(1+\mathrm{\Xi }_\mathrm{o}F)+{\displaystyle \frac{n(yzx)^2}{r^4}}`$ $`=`$ $`0,`$
$`{\displaystyle \frac{d(y(\sigma _{fi}+\sigma _{fn})v_f^2)}{ydy}}{\displaystyle \frac{z}{n}}(1+\mathrm{\Xi }_\mathrm{o}F)=0,`$ (8)
$`{\displaystyle \frac{d(y\mathrm{\Phi }_r)}{ydy}}{\displaystyle \frac{(xz+y)(yzx)}{nr^4}}=0,`$
$`{\displaystyle \frac{d(y\sigma _{fi}v_f)}{dy}}y=0,`$
$`{\displaystyle \frac{d(y\sigma _{fn}v_f)}{dy}}y\mathrm{\Xi }_\mathrm{o}F=0,`$
where $`\sigma _{fn}`$ is the surface density of the neutral component in the external layer. It is convenient to parametrize $`F`$ using a power law for the cross section in the typical range of pulsar velocities (Newman & al. newman82 (1982), McClure mcclure66 (1966), Peek peek66 (1966)):

$$\sigma _{cx}(v_{rel})=\sigma _{cx}(v_{rel,0})\left(\frac{v_{rel}}{v_{rel,0}}\right)^{0.37},$$
(9)
(where we have taken $`v_{rel,0}=1`$) so that:

$$F=1F_\mathrm{o}^{(v_{rel}(y)^{0.63}\sigma _{fi}n/1.25)}.$$
(10)
Therefor $`1F_\mathrm{o}`$ corresponds to the fraction of neutrals captured at the stagnation point, in the fully ionized case, which defines the properties of the true physical nebula (the factor $`1.25`$ corresponds to the surface density of the external layer on the axis, in the freely crossing case).

Contrary to the freely crossing case, Eqs. 8 cannot be integrated from the head, because they are not closed (i.e. in the Taylor expansion of the equations the terms of order $`𝒪`$($`y^n`$) are function of those of order $`𝒪`$($`y^{n+1}`$)). Instead we have used a perturbative method: from the freely crossing solution, we evaluate a first approximation of $`F`$; this approximation is then used in Eqs. 8 instead of the self-consistent expression, both to obtain the Taylor expansion and to integrate numerically the equations; we repeat iteratively till $`F`$ converges.

In Fig. 2 we present the solution of Eqs. 8 for the case $`\mathrm{\Xi }_\mathrm{o}=1`$ (an intermediate case between the highly ionized and highly neutral ISM) for different penetration thicknesses in the head. The solid line represents the freely crossing case (it is a good approximation if the thickness remains $`1`$ all along the nebula), while the dotted line refers to the completely dragging case (which is simply a bow shock scaled for a greater external density). The stagnation point distance from the star decreases with increasing values of $`F_\mathrm{o}`$, and its value is:

$$dd_\mathrm{o}/\sqrt{(1+\mathrm{\Xi }_\mathrm{o}(1F_\mathrm{o}))}.$$
(11)
The self-consistency level of our solution decreases moving to the tail of the nebula, away from the star, because the thin layer approximation is no longer acceptable, an so the constant (perpendicularly to the bow-shock) internal state assumption lose its validity.

Fig. 3 shows the surface density of the ionized component in the$`\mathrm{\Xi }_\mathrm{o}=1`$ case: also here the solid line represents the freely crossing case while the dotted line is the completely dragging solution. It is evident that the solution with highest $`F`$ tends to the dotted line in the tail of the nebula. Fig. 4 gives the surface density of the neutral component in the various cases obtained from Eqs. 8. Once the neutral and the ionized component density has been obtained it is quite easy to derive the H$`\alpha `$ luminosity along the nebula. However, before doing this it is necessary to consider the effect of other processes like diffusion and ionization (see Paper I) which modify the neutral distribution.

### 2.2 Importance of diffusion and ionization and their inclusion in the model

As previously mentioned the neutral H atoms, once stopped by charge-exchange with thermalized protons of the external layer, may diffuse or be ionized as a consequence of processes of charge-exchange and ionization in the internal hot plasma. Before introducing these effects in the equations it is necessary to discuss the physical regimes of the pulsar bow shock nebula, and their variations along the nebula. Comparing the time scale of the flow $`\tau _{flux}d_\mathrm{o}/v_\mathrm{o}`$ with, respectively, the time scale of diffusion and ionization, tells what are the importance of these two processes in the various pulsar cases.

The time scale for diffusion is:

$$\tau _{diff}\frac{\mathrm{\Delta }\sigma _{fi}\sigma _{cx}(v_{th})}{v_{th}},$$
(12)
where $`\mathrm{\Delta }`$ is the geometrical thickness of the external layer, the cross section of diffusion $`\sigma _{cx}`$ is that of charge-exchange, and the thermal velocity $`v_{th}`$ is that obtained from the Rankine-Hugoniot condition on the bow shock $`0.8v_\mathrm{o}\mathrm{sin}\alpha _1`$. Moving away from the head this time scale tends to increase, so that if it is greater than the flow time on the head, diffusion may not be taken into account. The time scale for ionization is:

$$\tau _{ion}\frac{\mathrm{\Delta }}{\sigma _{fi}\sigma _{ion}(v_{th})v_{th}},$$
(13)
where $`\sigma _{ion}`$ is the cross section of ionization, that, in the range of pulsar velocities, can be approximated by the relation (Ptak & Stoner ptak73 (1973)):

$$\sigma _{ion}(v_{th})=\sigma _{ion}(v_{th,0})\left(\frac{v_{th}}{v_{th,0}}\right)^{2.9}.$$
(14)
Also this time scale increases in the tail of the nebula, mainly as a consequence of the decreasing of the thermal velocity.

In the case of very luminous pulsars it is also necessary to determine the relative importance of diffusion and ionization, because atoms are ionized before the may diffuse (Paper I). The ratio $`\tau _{diff}/\tau _{ion}`$ tends to decrease away from the head because the cross section of ionization has a steeper velocity dependence than that of charge-exchange.

The effect of diffusion and ionization may be included in the equations to give a more realistic model, but some simplifications are necessary. We have assumed that atoms in the external layer have the same probability to diffuse in the internal stellar wind region as to diffuse in the external region of unperturbed ISM. In reality there may be an asymmetry in the diffusion not taken into account in our model, however the results of its inclusion are not significatively different, both in terms of surface densities and shape of the nebula. Diffusion acts in two ways: 1- reducing the surface density of neutral atoms, 2- bringing away from the layer an amount of tangential momentum, proportional to the number of diffused atoms. To treat properly diffusion we need to know the distribution of neutrals in the layer, but thin-layer models lack this piece of information. Contrary to diffusion, ionization does not change the total surface density in the external layer, but simply changes neutrals into ions. This may modify the geometry of the nebula because the increased ionized surface density makes the layer thicker to the penetration of neutrals. Using Eqs. 9 and 14 for the parametrization of the cross sections we find that the rate of diffusion $`C`$ and ionization $`B`$ may be written:

$`C`$ $`=`$ $`C_\mathrm{o}{\displaystyle \frac{\sigma _{fn}\mathrm{\Delta }_\mathrm{o}}{n^{0.37}\sigma _{fi}\mathrm{\Delta }}},`$ (15)
$`B`$ $`=`$ $`B_\mathrm{o}{\displaystyle \frac{\sigma _{fn}\sigma _{fi}\mathrm{sin}\alpha _{1}^{}{}_{}{}^{3.9}\mathrm{\Delta }_\mathrm{o}}{\mathrm{\Delta }}},`$ (16)
where all the unknown parameters are enclosed in the quantities $`C_\mathrm{o}`$ and $`B_\mathrm{o}`$, and $`\mathrm{\Delta }_\mathrm{o}`$ is the thickness of the layer along the axis in the freely crossing case. From a numerical simulation (Paper II) we find the parametrization $`\mathrm{\Delta }=\mathrm{\Delta }_\mathrm{o}d(1+y^2/5)/d_\mathrm{o}`$. Eqs. 8 then become:

$`((\sigma _{fi}+\sigma _{fn})v_f^2+\mathrm{\Phi }_r){\displaystyle \frac{dz}{dy}}n(1+\mathrm{\Xi }_\mathrm{o}F)+{\displaystyle \frac{n(yzx)^2}{r^4}}`$ $`=`$ $`0,`$
$`{\displaystyle \frac{d(y(\sigma _{fi}+\sigma _{fn})v_f^2)}{ydy}}{\displaystyle \frac{z}{n}}(1+\mathrm{\Xi }_\mathrm{o}F)+Cv_f=0,`$ (17)
$`{\displaystyle \frac{d(y\mathrm{\Phi }_r)}{ydy}}{\displaystyle \frac{(xz+y)(yzx)}{nr^4}}=0,`$
$`{\displaystyle \frac{d(y\sigma _{fi}v_f)}{dy}}y(1+B)=0,`$
$`{\displaystyle \frac{d(y\sigma _{fn}v_f)}{dy}}y(\mathrm{\Xi }_\mathrm{o}FCB)=0.`$
They can be integrated in the same way as discussed in the previous subsection, in fact the quantities $`C`$ and $`B`$ may be evaluated self-consistently.

Figs. 5 and 6 show respectively the effect of diffusion and of ionization on the surface density of the neutral component in the external layer, for various values of the coefficient $`C_\mathrm{o}`$ and $`B_\mathrm{o}`$. Both diffusion and ionization tend to reduce the density of H atoms. However this effect is reduced in the tail: diffusion decreases both as a consequence of the greater thickness of the layer and as a conseguence of the reduction of the mean free path of neutral H atoms; ionization simply vanishes because the cross section decreases rapidly with the thermal velocity. Fig. 7 shows the surface density of the ionized component: in case of hight ionization rates it tend to decrease in the head with respect to the values on the axis. This influences the H$`\alpha `$ luminosity, making it less that what may be derived from a simple case with only dragging charge-exchange.

## 3 H$`\alpha `$ luminosity

Once a neutral component has been introduced in the model, it may be used to evaluate the H$`\alpha `$ luminosity of the bow shock nebula and its variation moving away from the head. In Paper I, we have given an approximate estimate of such a luminosity, using some simplifications whose validity were limited to the stagnation point. The knowledge of the global distribution of neutrals and ions in the nebula makes it possible to evaluate the variation of luminosity along the bow shock.

Before investigating the effect of diffusion and charge-exchange we need to discuss the distribution of neutrals in the external layer. In those cases in which the penetration thickness for H atoms if not too high, we assume that the neutrals are homogeneously distributed in the external layer so that the luminosity per unit volume is constant in the direction normal to the layer. But in the highly dragging case it is possible that neutrals cannot reach the contact discontinuity, and remain confined in a thinner shell near the bow shock (al least in the head of the nebula). Anyway only a very high resolution mapping may discriminate among the above cases. In fact the instrumental resolution may spread the pattern. This is a problem if the signal to noise ratio is low so that only the more luminous regions of the nebula are detected.

To investigate properly the emission from such nebulae we have to consider the various processes responsible of the excitation of H atoms. They may be divided in two types: primary processes, which involve neutral atoms coming directly from the ambient medium; secondary processes, related to atoms dragged and thermalized in the external layer. In a non-radiative shock moving through a partially neutral medium the dominant processes for optical emission are exciting charge-exchange and collisional excitation of neutral H atoms, which give the observed Balmer spectrum. However primary processes usually have a very low efficiency for incoming neutrals (Paper I) so, the nebula can easily be detected only if the emission is dominated by the secondaries, otherwise a large amount of neutrals in the ISM is required to give enough flux. In what follows only secondary charge-exchange and ionization are considered.

Fig. 8 shows the behaviour of peak luminosity along the nebula for secondary charge-exchange. The various curves refer to the same cases of panel A in Fig. 5. There are various effects that contribute to the shape of the H$`\alpha `$ luminosity. The solid line represents an upper limit while the others show the effect of diffusion. The luminosity tends to increase as a consequence of the increasing column length of integration, but this effect is reduced by the diminishing of the ionized and neutral density, and this gives the behaviour of the various lines. While the case without diffusion gives a quite constant emission along the nebula, the presence of diffusion reduces the H$`\alpha `$ luminosity in the head (where the efficiency of diffusion is grater), and this effect could be used as a tool to determine the physical conditions of such nebulae.

We have not taken into account the presence of a threshold for the excitation of H atoms: we considered a pulsar so fast with respect to the ambient medium that even in the tail the thermal velocity is above the excitation energy ($`10`$ eV). In the case of a slower pulsar, the temperature in the tail may not be enough to allow exciting charge-exchange, so the H$`\alpha `$ flux may be lower. This is the case of PSR J0437+4715 (Bell & al. bell95 (1995)): only the head of the bow shock nebula have been detected so far despite the fact that this is the nearest object. In fact its velocity is estimated to be $`80`$ km/s, very near the threshold, while the other three nebulae, associated with faster pulsars, are detected even in the tail.

Fig. 9 shows the behaviour of peak luminosity along the nebula for secondary collisional excitation. The various lines refer to the same cases of panel A Fig. 5. Only models with ionization are shown. The H$`\alpha `$ emission rapidly decreases in the tail because the rate of interaction is a steep function of the temperature, and vanishes in the tail. So it is possible to discriminate between a flux dominated by charge-exchange or by collisional processes, analyzing the H$`\alpha `$ luminosity along the nebula. However even if the efficiency of collisional excitation is one-two orders of magnitude higher than that of charge-exchange, its cross section is much less at typical pulsar velocities, so that only in the faster cases ionization may give a non neglectable contribution, at least in the head of the nebula.

## 4 Conclusion

We have modified the two thin layer solution (Comeron & Kaper comeron98 (1998)) for bow shock produced by runaway stars (with the aim of deriving a model in the case of pulsar wind bow shock nebulae) with the inclusion of a neutral component. We find that the ISM neutrals lead only to minor effects to the shape of the nebulae (see Fig. 2), but plays an essential role in their luminosity. We have also introduced in our model the possibility that the H atoms, once stopped by the thermalized protons, may diffuse out of the layer or be ionized, as a consequence of charge-exchange and collisional ionization. Both processes act depleting the neutral surface density in the external layer, but their efficiency is higher in the head of the nebula both as a consequence of the decrease of the rate of charge-exchange and ionization in the tail, and the increase of the geometrical thickness of the layer.

The distributions obtained have been used to evaluate the Balmer emission along the nebula. The study of such emission can be a good tool to understand the properties of local ISM as well as those of the pulsar. Even if this is not a simple task, and requires very good measurements, which only for two objects are actually available, an accurate study may show which hypothesis can be applied (like, for example, the absence of preionization both from the pulsar surface radiation, or from the emission of shocked ISM). In Paper I, we stressed the point that the total H$`\alpha `$ flux depends on the velocity of the pulsar.

Here we evaluate its variations along the nebula. The behaviour is different if the dominant process is exciting charge-exchange or collisional excitation, so not only the total flux but also its variation moving away from the stagnation point can be used. Concerning the emission we find that: 1- charge-excange and collisional excitation have different trends; 2- processes like diffusion or ionization may play an essential role changing the emission in the head of the nebula. This not only allows one to put stronger constraints on the regime of such nebulae, but it may be used to verify the validity of the model, by comparison with the observed cases.

###### Acknowledgements.

I want to acknowledge Rino Bandiera for the fruitful discussion and the reading of this article. This work has been supported by the Italian Ministry for University and Research (MIUR) under the grant Cofin 2001-02-10.


format: text

text:
# LA-UR-02-3761 The Software Anatomy of a Flexible Accelerator Simulation Engine

## 1 Background

### 1.1 Discovering a Simulation Architecture

Our group has designed and implemented a unified accelerator Application Programming Interface (API) called XAL. XAL is designed to aid in the development of science control applications for beam physics. Accordingly, the XAL API is a physics-centric software programming interface. The physics applications interact with a model of an accelerator that resides in computer memory. XAL also contains the software infrastructure that creates the accelerator model. XAL loads a text-based (XML) description of an accelerator and assembles software objects such that an accurate model of the accelerator exists in computer memory. XAL is based on UAL , the Unified Accelerator Library.

The original motivation for XAL was to provide an accelerator independent interface for applications to interact with I/O from a live accelerator. This allows physicists to write beam physics control applications (Orbit Correctors, Beam Profile Monitors, RF Tuners, etc.) to the XAL API so that they can run on any accelerator. Some pseudo-code illustrating the principles of an XAL-based Orbit Correction Application may illustrate the essence of the concept.

```
Accelerator  theAccel = XALFactory.newAccel(‘‘sns.xml’’)
BPM[]        theBPMs       = theAccel.getNodesOfType(BPM)
HorzDipole[] theCorrectors = theAccel.getNodesOfType(DCH)
for each BPM in theBPMs
      read BPM.avgPos() and set a corrector magnet accordingly

```

To aid in writing applications that take into account design values, the accelerator description file contains all design information for the accelerator. This condition allows, for example, a physics application to compare the design field of a quadrupole with its read-back (runtime) field.

With all design information incorporated into a software model of an accelerator, we have discovered an excellent simulation engine. As long as the software accelerator has a convenient means for traversing beam-line devices in a spatially sequential manner, we can use design values along the way to simulate beam-dynamics.

This scenario allows for a drastic departure from traditional accelerator simulation codes. Traditionally simulators have been isolated software products. They load some type of lattice description of an accelerator and apply predefined beam-dynamics to an initial beam. Ultimately this design has led to huge codes (to account for various beam-line element types). Further, these codes typically operate with only one type of simulation (multi-particle or rms envelope, but not both).

The architecture presented here contains a novel approach to the simulation domain. It is our conjecture that the method presented here better captures reality in that there is some sort of software beam actually traversing a software model of a real accelerator.

### 1.2 The Architecture

Our approach is based upon the Element-Algorithm-Probe Design Pattern . The core concept of this design pattern is the separation of beam-dynamics code from the actual beam-line elements. It is desirable to keep the code that corresponds to beam-line elements as simple as possible so that the application writer has a clean interface to a beam-line element. The Element-Algorithm-Probe pattern enforces this concept by requiring beam-dynamics code to exist in a separate entity, called an IAlgorithm .

Deferred until runtime is the binding of beam-dynamics to actual beam-line elements. This deployment strategy allows for conceptually correct simulations. First it is truly modular. The three concepts, beam-line elements, beam-dynamics, and the beam are compartmentalized into separate code. Second it is truly maintainable. To support a new beam type or new beam-line element type does not cause code bloat. Finally it is truly extensible. Via the mechanism of a Java interface, various beam-dynamics algorithms can be written for the same type of beam-line element and switched at will at runtime. Modularity, maintainability, and extensibility provide true power and flexibility to our architecture.

## 2 IProbe, IAlgorithm, and IElement

### 2.1 Technology Introduction

It may help to understand the facets of Java that we exploit in order to implement the Element-Algorithm-Probe Design Pattern.

At the center of the Element-Algorithm-Probe pattern is the concept of a Java interface. Essentially, an interface is a contract between a user and an implementor. The contract says that the implementor of an interface is required to provide an implementation of the methods defined in the interface.

For example, consider the interface

```
public interface Thermometer {
   public double getTemperature();
}

```

Using this interface, a programmer can assume being able to perform operations on a thermometer no matter how the thermometer actually obtains the temperature. This is desirable because a thermometer implementor can change how the temperature is actually obtained (if, say, a new sensor system was installed) without requiring all thermometer users to recompile their code.

We use the same idea with beam-dynamics code. Beam-dynamics reside in files that implement (the computer science term for acknowledging involvement in the contract from the implementors point of view) the IAlgorithm interface. Since the simulation engine knows how to do beam-dynamics calculations solely by interacting with IAlgorithms, it is trivial to swap beam-dynamics algorithms at will.

The IAlgorithm interface looks like this.

```
public interface IAlgorithm {
   public void propagate(IElement, IProbe);
   public Class legalElementType();
   public Class legalProbeType();
}

```

Conceptually, an IAlgorithm implementor is required provide an implementation of the method propagate() to modify the the beam (IProbe) according to the beam dynamics of the the beam-line element (IElement).

In essence, all that the simulation engine knows about are the three data types (all defined in interfaces) IAlgorithm, IProbe, and IElement.The beauty of this design is that there are separate code locations for beam-line elements, beam-dynamics, and the beam itself.

### 2.2 Probes

The IProbe interface should ideally contain the bare minimum information to fully represent a beam. Such beam information consists of beam current, beam charge, particle charge, particle rest energy, particle kinetic energy, etc. Further, since a probe represents the state of the beam at a position in the beam-line, a probe also contains a beam-line-position attribute. The current IProbe specification serves the purpose of representing a beam for a single particle , particle ensemble, and envelope simulations (in both two and three dimensions). Figure 1 represents a suitable inheritance hierarchy of probe types to handle the aforementioned simulation types.

### 2.3 Propagation

It is important to note that there are various approaches toward simulating beam-dynamics. For example, accurate approaches may involve slicing nodes up into small pieces. An aggregate of approximations done on sufficiently small elements is typically more accurate than one overall approximation. However, normally this is only practical in elements that have special behavior. So the question arises: How are probes propagated through elements?

It is the responsibility of the IAlgorithm implementor to handle all beam-dynamics, including the propagation mechanism. Sample propagation mechanisms will be presented later in this paper. However, keep in mind the most appropriate propagation mechanism when implementing algorithms for the particular problem at hand.

### 2.4 Algorithms

We have already introduced the concept of the IAlgorithm interface. Now let us pursue a few more details regarding its implementation.

The IAlgorithm interface provides a generic way of assembling algorithms in a simulation engine. In practice any particular IAlgorithm implementation only makes sense in the context of a particular beam-line element type and probe type. For example, a hypothetical IAlgorithm called QuadParticleMapper would expect a Quadrupole as its IElement and a Particle as its IProbe . Providing such specificity is the job of the legalElementType() and legalProbeType() methods. An implementation of the QuadParticleMapper could look like this.

```

public class QuadParticleMapper(IElement p_elem,
                                IProbe p_probe){
  public Class legalElementType(){return Quadrupole.class;}
  public Class legalProbeType(){return Particle.class;}
  public void propagate(){Quadrupole/Particle beam dynamics}
}

```

By providing these methods, the simulation engine can do type checking upon algorithm binding. It would not make sense to bind this algorithm to a WireScanner. Providing these methods helps to avoid that condition.

## 3 Design of a Single Particle Simulation

Designing an actual simulation merely involves putting together Elements, Algorithms, and Probes in a semantically meaningful way. It turns out that, to the first order, the beam dynamics through a particular node type can be captured by a transfer matrix. This property allows for a straight-forward means of simulating a particle traveling down a beam-line. An object-oriented approach would be to create a ParticleMapper class that transforms the Particle probe by the simple vector-matrix multiplication

$$\stackrel{}{𝐳}_{n+1}=\mathrm{\Phi }_n\stackrel{}{𝐳}_n$$

where $`\stackrel{}{𝐳}_n`$ is the coordinate vector of the particle $`\left(\begin{array}{cccccc}x& x^{}& y& y^{}& z& z^{}\end{array}\right)^T`$ at the start of the node, and $`\mathrm{\Phi }_n`$ is the transfer matrix of the node. Further, $`\mathrm{\Phi }_n`$ can be obtained by the ParticleMapper via the use of an abstract method that is implemented by beam-dynamics algorithms for individual nodes (QuadrupoleParticleMapper, RFCavityParticleMapper, etc.). A suitable class design can be seen in Figure 2.

To further illustrate some of these concepts, the basic layout of the ParticleMapper class looks like this.

```
abstract public Matrix computeTransferMatrix();

public void propagate(IElement pElem, IProbe pProbe)
{
        //type-cast the probe and element to what we expect
        Particle theProbe = ((Particle)pProbe);
        AcceleratorNode theNode = ((AcceleratorNode)pElem);Ψ

        //do the vector-matrix multiplication
        theProbe.setCoords(computeTransferMatrix
                             .times(theProbe.getCoords()));

        //advance the probe the length of the node
        theProbe.advancePosition(theNode.getLength());
        return;
}

```

Once the computeTransferMatrix() operations are implemented for the node-specific dynamics, all that remains is writing a driver program. A driver program binds algorithms to nodes and injects the probe. Here is a pseudo-code driver111It is anticipated as a logical extension to these ideas that an AlgorithmManager GUI will be written so that simulations can be bound at run-time to the XAL accelerator model and easily hooked into analysis software..

```

//instantiate the XAL accelerator model
Accelerator  theAccel = XALFactory.newAccel(‘‘sns.xml’’)

//bind the algorithms
Quadrupole[] theQuads      = theAccel.getNodesOfType(QUAD)
RFCavity[]   theCavities   = theAccel.getNodesOfType(RFC)
for each QUAD in theQuads
      bind a QuadParticleMapper instance to QUAD
for each RFCav in theCavities
      bind a RFCavityParticleMapper instance to RFCav

//instantiate a probe
Particle p1 = new Particle(initial conditions...)

//run the probe down the beam-line
AcceleratorNode[] theNodes = theAccel.getAllNodes();
for each NODE in theNodes
      NODE.propagate(p1)

```

And that is it!111Notice that beam-dynamics between nodes are left out of the demonstration case. In the actual simulation engine, the propagate() method accounts for space between the position of the probe and the start of the node by calculating the beam-dynamics through a drift space. However future implementations of XAL may include drift-spaces as an actual node type which would warrant writing a specific DriftSpaceMapper class.

The particle probe will be transformed by each beam-line element according the the bound algorithm. Note that the pseudo-code is a basic proof of concept and does not contain the code necessary to broadcast probe increment intermediate data to produce, for example, a plot.

The single particle simulation can be applied to a two-dimensional case by only considering the first four elements of $`\stackrel{}{𝐳}`$. Further, the single particle simulation can be extended to a multi-particle simulation (in two and three dimensions) by constructing a container of particle probes and writing beam-dynamics algorithms that properly transform the collection. The only matter that complicates (and complicate it does!) a multi-particle simulation is the concept of space-charge. Before biting off this task, however, a presentation of another type of simulation that accounts for space-charge is warranted.

## 4 Design of an RMS Envelope Simulation

### 4.1 The Concept

The RMS qualities of a beam can be represented by the 6x6 symmetric matrix $`\sigma `$ that statistically expresses the boundaries of a beam in transverse, longitudinal, and phase space by using moments of the beam distribution. RMS Envelopes are convenient because applying beam-dynamics involves a simple matrix operation. Namely, the same transfer matrix $`\mathrm{\Phi }`$ used in single particle simulations can propagate rms envelopes according to the conjugation

$$\sigma _{n+1}=\mathrm{\Phi }\sigma _n\mathrm{\Phi }^T$$

The other important concept in this simulation is space charge. A $`\sigma `$ matrix is a statistical representation of a beam, which is a multi-particle entity. Therefore, each particle in the beam is aware (electromagnetically) of all other particles in the beam. It turns out that to the first order the effects of space charge can be captured in a $`\mathrm{\Phi }`$ matrix. While it may not be mathematically trivial to calculate the matrix, having the calculation in such a form makes the integration into our simulation engine simple. However it should not be overlooked that this quantity is very important to the correctness of simulation.

### 4.2 Propagation

The envelope simulation is more complex than a single-particle simulation in that we will propagate envelopes through elements using more than one propagation mechanism. Specifically, we may be able to compute a better approximation of behavior through quadrupoles than RF Cavities.

This condition is due to the fact that the $`\mathrm{\Phi }`$ matrix for a quadrupole adheres to the semi-group property.

$$\mathrm{\Phi }(\mathrm{\Delta }s_1+\mathrm{\Delta }s_2)=\mathrm{\Phi }(\mathrm{\Delta }s_1)\mathrm{\Phi }(\mathrm{\Delta }s_2)$$

or

$$\mathrm{\Phi }(n\mathrm{\Delta }s)=\mathrm{\Phi }^n(\mathrm{\Delta }s)$$

where $`\mathrm{\Delta }s`$ is the length of the quadrupole being considered.

To more accurately consider space charge, we take advantage of the semi-group property of the transfer matrices. In the propagate() method of the SemiGroupEnvelopeMapper (see Figure 3) we subsection the node (e.g., a quadrupole) into $`n`$ slices of length $`\mathrm{\Delta }s=\frac{l}{n}`$ where $`l`$ is the length of the quadrupole. Then we run the probe through these $`n`$ slices, applying space charge kicks after every subsection (See Figure 4).

Since RF Cavity transfer ($`\mathrm{\Phi }`$) matrices do not in general adhere to a semi-group property, we are forced to take a more simplistic approach toward transforming the envelope. We will slice the node in two, treating each half as a drift-space (to account for space charge) and hit the envelope in the middle of the node with the numerically approximated $`\mathrm{\Phi }`$ matrix (see Figure 5).

## 5 Design of a Particle Ensemble Simulation

As a final exercise it will be useful to consider the design of a multi-particle simulation. The true complication of designing a multiple-particle (ensemble) simulation is the computation of space-charge effects. Unfortunately, to model multiple particles, space-charge effects cannot be accurately captured by a transfer matrix. On the other hand, the architecture outlined in this paper keeps the details of the space-charge calculations from interfering with code cleanliness.

```

```

The two core concepts of a multi-particle (ensemble) simulation are

* Calculation of the electric self-fields of the ensemble
* Using the calculated fields to update the particle coordinates.

There are various approaches that can be taken for both tasks. All that we attempt to show here is that by correctly isolating these concepts, a clean software architecture can be maintained.

Namely, an Ensemble probe should encapsulate the logic necessary to obtain the electric self fields of the ensemble. That being the case, various Ensemble probe implementations could be swapped at will to employ different field calculation techniques. For example, many electric field calculation techniques involve solving Poisson’s equation to obtain the electric potential of the ensemble. By hiding this code from the simulation engine (it is contained within the Ensemble probe implementation), the implementor could exploit parallel processing facilities (See Figure 6).

By moving the calculation of electric fields out of the beam-dynamics code, the beam-dynamics algorithm developer is free to choose space-charge consideration techniques with minimal impact to code clarity. One may decide to take the “thin lens kick” approach that has been used previously in this paper. One may alternatively decide to apply a “trajectory integration” based approach. The key point here is that by separating codes into their logical components allows for a high degree of flexibility in simulation technique.

## 6 Conclusions and Future Directions

We are enthusiastic to report that the results obtained in the RMS Envelope Simulation have been validated against Trace3D  Figure 7111The slight inconsistencies arise from the fact that Trace3D only outputs at the end of nodes whereas our simulation engine produced a higher resolution of intermediate data points. shows agreement between simulation results of the SNS Medium Energy Beam Transport (MEBT) using both Trace3D and the XAL simulation engine. It is encouraging that a problem domain with so many interdependencies(particle physics) can be simulated with a clean architecture.

As we move toward the future, we are anticipating the ability to implement model reference control techniques. That is, within the XAL model there is access to a live accelerator and a simulated accelerator. Having both at hand allows the comparison of live behavior with simulated behavior to develop control strategies.

The key to effectively implementing an environment conducive to model reference control is architectural discipline when designing both the I/O and simulation aspects of XAL. As long as the interface to the two are respectively clean, hybridization of the two will be a straightforward extension.


text:
# Convergence and stability in numerical relativity

## Abstract

It is often the case in numerical relativity that schemes that are known to be convergent for well posed systems are used in evolutions of weakly hyperbolic (WH) formulations of Einstein’s equations. Here we explicitly show that with several of the discretizations that have been used through out the years, this procedure leads to non-convergent schemes. That is, arbitrarily small initial errors are amplified without bound when resolution is increased, independently of the amount of numerical dissipation introduced. The lack of convergence introduced by this instability can be particularly subtle, in the sense that it can be missed by several convergence tests, especially in 3+1 dimensional codes. We propose tests and methods to analyze convergence that may help detect these situations.

Convergence is a central element of any numerical simulation. It refers to the property that if one refines the simulation by adding more points to the grid, numerical errors should diminish. In the limit of zero spacing, they should go to zero and one should get the exact solution. Having a convergent code is a key element for numerical simulations to have predictive power: although one will always be limited in practice to a finite number of points in the grid, one can extrapolate the results for more and more refined simulations and have a very good approximation to the true (continuum) results. Codes that do not converge produce answers that, even if they remain finite (at least for a while), have little predictive power: there generically is no way to know if the solutions found approximate the desired continuum solution.

In this paper we want to emphasize that discretization schemes that yield convergent code for strongly hyperbolic (SH) systems of equations do not necessarily do so for weakly hyperbolic (WH) or completely ill posed (CIP) systems. The relevance of this observation is that several formulations of the Einstein evolution equations commonly used in numerical relativity, including, for example, the Arnowitt–Deser–Misner (ADM) ADM  and Baumgarte–Shapiro–Shibata–Nakamura (BSSN) BSSN  formulations with fixed lapse and shift, are not SH.

We concretely prove that a simple system of linear WH equations with constant coefficients, when discretized with either the iterated Crank–Nicholson (ICN) with fixed number of iterations or second order Runge–Kutta (2RK) methods, leads to unconditionally unstable codes, even if numerical dissipation is explicitly added. This simple system of equations is related to the equations one encounters in the WH formulations of general relativity. It therefore strongly suggests that these formulations should produce code that is not convergent. The lack of convergence is of a particularly pernicious nature, since in the WH case it might grow slower than the “usual” von-Neumann numerical instability. It is such that if one tests the code with stationary solutions (for instance simulating a single black hole) and performs convergence tests, one could easily be confused into believing one has convergent code.

Consider a linear system of $`m`$ partial differential equations in $`m`$ variables of the form

$$_tv=A_xv,$$
(1)
where $`A`$ is a constant $`m\times m`$ matrix and $`v=(v_1,\mathrm{},v_m)^T`$. The system is SH if $`A`$ has real eigenvalues and is diagonalizable, WH if it has real eigenvalues but is not diagonalizable, and CIP if it has complex eigenvalues.

The iterated Crank–Nicholson (ICN) method is an explicit approximation to the (implicit) Crank–Nicholson method first proposed by Choptuik. In this proposal the number of iterations is not fixed but, instead, might depend on the spatial gridpoint and time step 111We thank Matt Choptuik for discussions on this point.. A different notion of ICN is to a priori fix the number of iterations ICN . Teukolsky ICN  showed that for the advective equation, $`_tv=_xv`$ (which is SH), such scheme was stable if one used $`2,3,6,7,10,11,\mathrm{}`$ iterations (one iteration corresponds to the 2RK method). This paper seeks to extend this proof to more general systems of equations, including explicit artificial viscosity.

A difference scheme is said to be stable with respect to a given norm if there exist positive constants $`\mathrm{\Delta }x_0`$ and $`\mathrm{\Delta }t_0`$ and $`f(t)`$ such that $`u^nf(t)u^0,`$ for $`0<\mathrm{\Delta }x\mathrm{\Delta }x_0`$, $`0<\mathrm{\Delta }t\mathrm{\Delta }t_0`$, and $`t=n\mathrm{\Delta }t`$ where $`u^n`$ is the numerical solution at time $`t`$ and $`u^0`$ the initial data. That is, the growth of the solution is bounded by some function of time $`f(t)`$ that is independent of the initial data and resolution. Some schemes are only conditionally stable. This means that a relationship between $`\mathrm{\Delta }t`$ and $`\mathrm{\Delta }x`$ is needed for stability. Through Lax’s theorem, numerical stability in this sense is equivalent to convergence, provided the scheme is consistent.

Working in Fourier space and writing explicitly the solution at time $`t`$ in terms of the amplification matrix $`Q`$, $`\widehat{u}^n=Q^n\widehat{u}^0`$ yields a necessary condition for stability: the von-Neumann condition (VNC). It states that, for the cases considered in this paper, the spectral radius $`\rho (Q)`$ of $`Q`$ (the maximum eigenvalue in norm) is not greater than one. If this condition is not satisfied, a numerical instability that grows like $`\mathrm{exp}(n\mu )`$, with $`\mu >0`$ a constant, is present. A well known example of this instability is the forward in time, centered in space scheme for the advective equation. It is sometimes thought that the VNC is not only necessary but also sufficient for stability. This often leads to (wrongly) concluding that a discretization for a WH system is stable based on an discrete analysis for the advection or wave equation. The VNC is sufficient only if $`Q`$ can be diagonalized by a similarity transformation whose norm, and that of its inverse, has an upper bound that is independent of resolution kreiss . An observation that will be important for the stability analysis below is that, for the cases we are discussing, $`Q`$ is diagonalizable if and only if $`A`$ is. Therefore, in the WH case the VNC will not be sufficient for stability.

On the other hand, if the system is SH $`A`$ is diagonalizable (with real eigenvalues) and so is $`Q`$. Working in the diagonal basis we then have an uncoupled system of $`m`$ equations for the $`m`$ grid functions. Repeating Teukolsky’s analysis for the ICN scheme, this system is stable if $`\rho (A)\lambda 2`$ ($`\lambda =\mathrm{\Delta }t/\mathrm{\Delta }x`$ is the Courant factor) and the number of iterations is $`p=2,3,6,7,10,11,\mathrm{}`$, if no explicit dissipation is added. In evolutions of SH equations, the addition of dissipation may help stabilize schemes that are otherwise unstable. For example, the VNC for the 2RK case with third order dissipation is

$$(\rho (A)\lambda )^48\lambda \stackrel{~}{\sigma }1,$$
(2)
where $`\stackrel{~}{\sigma }`$ is the dissipation parameter. Since $`Q`$ is diagonal, this condition is necessary and sufficient for stability.

Below we will explicitly prove that for CIP or WH cases the ICN and 2RK methods are unstable, independently of the amount of dissipation, number of iterations and Courant factor. Although our results are only proven for a system of linear equations with constant coefficients, they are significant in several ways for numerical relativity. First of all, if one considers situations that are small perturbations off Minkowski spacetime, then a system of equations identical to the one we considered does appear in the linearized Einstein equations written in first order form in space and time if the system of evolution equations used is WH or CIP. Since far away from binary black holes the spacetime is approximately flat, this situation does arise in realistic simulations. Could non-linearities change our conclusions? It is possible, but unlikely. Usually when instabilities are established for linear systems the addition of non-linear terms only makes instabilities worse kreiss .

The type of lack of convergence differs in the CIP and WH case. In the CIP case the VNC is easily violated, codes crash very fast and non-convergence is easy to spot. In the WH case the VNC may be satisfied and crashes may occur very late (if they occur at all). If one performs routine convergence analyses in the region before the crash one might appear to see convergence, especially if the initial data has few (and low) frequency components, unless one increases enough the resolution or runs for long enough time.

We now show that these different types of lack of convergence do appear in numerical relativity simulations, even in very simple 1D ones with periodic boundary conditions. We have performed non-linear evolutions of a plane gauge wave spacetime (that is, the spacetime is flat, though in a time-dependent slicing),

$$ds^2=e^{A\mathrm{sin}(t+x)}(dt^2+dx^2)+dy^2+dz^2.$$
(3)
Figure 1 shows results of an evolution using 2RK with $`\sigma =1/2,0,1/2`$ in the Einstein–Christoffel (EC) system ec  of evolution equations. The latter is a symmetric-hyperbolic reformulation of Einstein’s evolution equations that includes a parameter $`\sigma `$ that densitizes the lapse. In the EC formulation, $`\sigma =1/2`$, but by tuning this parameter we can make the system WH ($`\sigma =0`$) or CIP ($`\sigma =1/2`$). The whole $`30`$ equations of the EC system are evolved, but all quantities are assumed to depend only on $`t`$ and $`x[\pi ,\pi ]`$. The Courant factor is set to $`1/2`$, the dissipation parameter to $`\stackrel{~}{\sigma }=0.02`$, and the evolution is followed for $`1,000`$ crossing times. In the SH case the code is convergent for all resolutions tested, see stability  for a detailed discussion. On the other hand, in the same code with $`\sigma =1/2`$ a lack of convergence becomes apparent immediately (before one crossing time), the errors become bigger and the code crashes earlier as resolution is increased. In the WH case the code is also not convergent but the effect is less noticeable, since with the chosen values of dissipation and Courant factor the VNC, Eq.(2), is satisfied.

If one performed only a few runs, with, say $`120`$ and $`240`$ gridpoints as is typical in 3D convergence studies, one would have to wait for around $`150`$ crossing times in order for the lack of convergence in this WH example to become obvious. To put these numbers in context, suppose one had a similar situation in a 3D black hole evolution. Suppose the singularity is excised, with the inner boundary at $`r=M`$ and the outer boundary at $`20M`$. In this case $`120`$ and $`240`$ gridpoints correspond to grid spacings of, approximately, $`M/5`$ and $`M/10`$, respectively. If one had to evolve up to $`150`$ crossing times in order to notice the lack of convergence, that would correspond to $`t3000M`$, which is more than what present 3D evolutions typically last. One could therefore be misled to think that the code is convergent. Repeating the same runs with an initial data with more frequencies (say, a non-stationary black hole) would make the instability manifest in a shorter timescale. By making a spatial Fourier decomposition (Figure 2) of the numerical solution, we have found that in the WH and CIP cases there are always non-zero frequency modes growing exponentially from the very beginning, though sometimes starting at truncation error. By performing such decomposition one can detect that the code is not convergent much before this becomes obvious in the overall errors.

We have done simulations with different spectral distributions on the initial data, different number of iterations in the ICN method, different Courant factors, and different values of dissipation. The time and resolution at which the numerical instability becomes obvious depends on all these factors, but lack of convergence is always present in the CIP and WH cases, while the SH runs do converge. Too much dissipation violates the VNC, resulting in a more severe instability and this is immediately seen in numerical experiments, see stability  for details.

We have also found similar results performing simulations with the same initial data but with Kidder-Scheel-Teukolsky’s many parameter family of formulations of Einstein’s equations kst , leaving all parameters fixed and changing only one of them at a time, achieving different levels of hyperbolicity. In particular, we have found lack of convergence in the ADM equations as well (which is a subsytem of this family).

The results presented in this paper do not represent a special feature of the ICN method but instead, only reflect the fact that the definition of numerical stability is just a discrete version of well posedness. Therefore, difference schemes approximating ill posed problems can be expected to be non-convergent. In the absence of boundaries this is the case for CIP or, generically, WH formlations. If there are boundaries, strong hyperbolicity is not enough and extra care has to be taken in order to guarantee well posedness; wrong boundary conditions not only lead to inconsistencies but also to lack of convergence (see, e.g., tadmor ).

Although it is not possible to prove in a definitive way that a code is convergent only by numerical experiments, the previous examples suggest some of the pathologies that one should look for. Namely, Fourier modes that are not convergent but that are very small for a while, since they are initially suppressed. Also, notice from the WH simulations above shown that a code does not need to crash nor exhibit violent growth in order not to converge. The main lesson learned from this paper is that one should exercise significant care in numerical simulations before empirically concluding that the simulation is convergent, especially if the formulation of Einstein’s equations used is not SH or its level of hyperbolicity is unknown.

We wish to thank D. Arnold, M. Choptuik, P. Laguna, L. Lehner, M. Miller, R. Price, O. Reula, and S. Teukolsky for comments. This work was supported in part by NSF grant PHY9800973, the Horace Hearne Jr. Institute for Theoretical Physics, the Swiss National Science Foundation, and Fundación Antorchas.

Proof of non-convergence for CIP and WH cases: A usual way of discretizing the right hand side of Eq.(1) is using centered differences plus third order explicit dissipation. By this one means solving $`_tv=Cv/\mathrm{\Delta }t`$, with $`C=A\lambda \delta _0/2\stackrel{~}{\sigma }\lambda \delta ^4`$, where $`I`$ is the identity matrix, $`\delta _0v_k=v_{k+1}v_{k1}`$, and $`\delta ^4v_k=v_{k+2}4v_{k+1}+6v_k4v_{k1}+v_{k2}`$. This is equivalent to discretizing $`_tv=A_xvI\stackrel{~}{\sigma }(\mathrm{\Delta }x)^3_x^4v`$ using centered differences (using first order dissipation in our proof, i.e. discretizing $`_tv=A_xv+\stackrel{~}{\sigma }\mathrm{\Delta }x_x^2v`$, gives similar results).

The ICN method with $`p`$ iterations for the partial differential equation (1) can be written as

$$v_k^{n+1}=\left(1+2\underset{j=1}{\overset{p+1}{}}\frac{C^j}{2^j}\right)v_k^n,$$
(4)
where $`v_k^n=(v_{1k}^n,\mathrm{},v_{mk}^n)^T`$. The index $`n`$ corresponds to the time step and $`k`$ to the spatial mesh point. Since we are considering an initial value problem on all of $`𝖱\text{ }`$, $`\mathrm{}<k<+\mathrm{}`$ (our proof can be easily modified for periodic boundary conditions, and the results are the same).

In order to analyze stability we work in the basis in which $`A`$ takes its canonical form. That is, one multiplies both sides of Eq.(4) by $`T`$, where $`TAT^1=J`$ has the canonical Jordan form, and analyzes the equation

$$u_k^{n+1}=\left(1+2\underset{j=1}{\overset{p+1}{}}\frac{L^j}{2^j}\right)u_k^n,$$
(5)
where $`L=TCT^1=J\lambda \delta _0/2I\stackrel{~}{\sigma }\lambda \delta ^4`$. Any conclusions regarding stability will hold as well for the original variable $`v`$. Working in Fourier space, $`u_k=_\pi ^\pi \mathrm{exp}(ik\xi )\widehat{u}(\xi )𝑑\xi `$, $`Q=1+2_{j=1}^{p+1}\widehat{L}^j/2^j`$ and $`\widehat{L}=J\lambda i\mathrm{sin}(\xi )16\stackrel{~}{\sigma }\lambda \mathrm{sin}^4(\xi /2)`$.

If the system is CIP, $`J`$ has at least one complex eigenvalue $`c`$, say $`J_{11}=c`$. In this case, $`\widehat{L}_{11}=c\lambda i\mathrm{sin}(\xi )16\stackrel{~}{\sigma }\lambda \mathrm{sin}^4(\xi /2)`$, and the eigenvalue $`Q_{11}`$ has norm $`1\lambda Im(c)\xi +O(\xi ^2)`$. Therefore, the VNC is violated for sufficiently small $`\xi `$ with the appropriate sign and the scheme is unconditionally unstable, as expected.

If the system is WH then $`J`$ has at least one Jordan block and so does $`Q`$. In this basis all the Jordan blocks are uncoupled and we can consider one at a time. We assume there is one block of dimension $`2\times 2`$ (the proof actually holds for higher dimensionality as well). The canonical form for such a block, the resulting amplification matrix, and its n-th power are

$$J=\left(\begin{array}{cc}c& 1\\ 0& c\end{array}\right),Q=\left(\begin{array}{cc}a& b\\ 0& a\end{array}\right),Q^n=\left(\begin{array}{cc}a^n& a^{n1}nb\\ 0& a^n\end{array}\right).$$

Suppose the VNC $`|a|1`$ is satisfied (this rules out the $`p=0`$ case), we will show that the scheme is still unstable. Consider as initial data a small perturbation of the form

$$u_0^0=(0,2\pi ϵ)^T,u_k^0=(0,0)^T\text{ otherwise}$$
(6)
with $`ϵ`$ arbitrary small. In Fourier space we have $`\widehat{u}^0(\xi )=ϵ(0,\mathrm{\hspace{0.17em}1})^T`$ for all $`\xi `$. The solution at time $`t=n\mathrm{\Delta }t`$ is, then, $`\widehat{u}^n=ϵ(a^{n1}nb,a^n)^T`$. We will show that its norm grows without bound when the number of gridpoints is increased (while keeping the time and Courant factor fixed). We first notice that ($`n_22n2`$),

$$\frac{u^n^2}{u^0^2}=_\pi ^\pi \left(n^2|a|^{n_2}|b|^2+|a|^{2n}\right)\frac{d\xi }{2\pi }n^2_\pi ^\pi |a|^{n_2}|b|^2\frac{d\xi }{2\pi }.$$
(7)
Expressions $`|a|^2`$ and $`|b|^2`$ are analytic functions of $`\xi `$ with Taylor expansion $`|b|^2=\lambda ^2[\xi ^2+O(\xi ^4)]`$ and $`|a|^2=1[2\stackrel{~}{\sigma }\lambda +d(c\lambda )^4]\xi ^4+O(\xi ^6)`$, with $`d=1`$ for $`p=1`$, $`d=1`$ for $`p=2`$, and $`d=0`$ otherwise. From these expansions one can see that, for all $`p`$, there are positive constants $`\alpha `$ and $`\rho `$ such that

$$|a(\xi )|^21\alpha \xi ^40,|b(\xi )|^2\rho ^2\xi ^20,$$
(8)
for small enough $`\xi `$, say $`|\xi |ϵ_1`$. For $`n>1`$ and small enough $`\xi `$, say $`|\xi |ϵ_2`$,

$$|a(\xi )|^{2n}1n\alpha \xi ^40,$$
(9)
where the last inequality holds for $`|\xi |(n\alpha )^{1/4}`$, provided that $`n`$ is large enough so that $`(n\alpha )^{1/4}\mathrm{min}(ϵ_1,ϵ_2)`$. Using bounds (8,9) and the VNC in inequality (7) gives

$$\frac{u^n^2}{u^0^2}n^2_{(n\alpha )^{1/4}}^{(n\alpha )^{1/4}}(1n\alpha \xi ^4)\rho ^2\xi ^2\frac{d\xi }{2\pi }=\frac{4}{21\pi }\frac{\rho ^2}{\alpha ^{3/4}}n^{5/4}$$

which diverges for $`n\mathrm{}`$. This shows that in the WH case the ICN and 2RK schemes are unstable, independent of the amount of dissipation. In fact, by examining $`|a(\xi =\pi )|`$, one can see that a necessary condition for the VNC (for any number of iterations) is $`0\stackrel{~}{\sigma }\lambda 1/8`$. Adding too much dissipation violates the VNC and the instability becomes worse.

Perturbations like that of (6) are to be expected in a numerical simulation due to truncation or roundoff errors. For high enough resolution, such a perturbation will be amplified without bound, spoiling any convergence. One expects that in the non-constant coefficient or in the non-linear case the rate of growth of the instability with number of gridpoints will be even faster (see example in page 216 of kreiss ), as in figure 2.


text:
# The effect of disk magnetic fields on the truncation of geometrically thin disks in AGN

## 1 Introduction

Accretion onto supermassive black holes can occur via an advection-dominated accretion flow (ADAF) or via a standard geometrically thin accretion disk. The mode depends on the accretion rate. In an ADAF advection rather than radiation removes the locally generated accretion heat. This is only possible in the inner region around the central accretor. At larger distance the mass accretes via a cool geometrically thin standard disk. The situation is the same for supermassive and galactic black holes (for a review on the now already “classical” advection-dominated accretion and applications see Narayan et al. 1998). Confirmation for this picture, an ADAF in the inner region near the black hole and a surrounding cool standard disk with a corona above came from the successful modeling of the spectral transitions of low-mass X-ray binaries (LMXRBs) which arise from the changing mode of accretion (Esin et al. 1997, 1998).

The change from disk accretion to an ADAF was investigated for several low luminosity AGN. The low luminosity in the elliptical galaxy M87 (NGC 4486) was pointed out by Reynolds et al. (1996). Quataert et al. (1999) investigated the low-luminosity galactic nuclei of M81 and NGC 4579 and from the spectral fits found evidence for thin disks truncated at about 100 Schwarzschild radii, similar to the result for NGC 4258 (Gammie et al. 1999). Di Matteo et al. (1999, 2000) studied the nuclear regions of six elliptical galaxies and derived larger truncation radii from the X-ray and high-frequency radio observations. In these fits wind loss was included which seems to be demanded at least in some galaxies. This means a significant fraction of mass, angular momentum and energy is removed from the accretion flow by a wind. Large radii were found also for M 87 and NGC 4696 by Quataert et al. (1999) (see detailed discussion of constraints for the different objects in their paper). That mass loss via wind may be important was also stated by Quataert & Narayan (1999), who pointed that different combinations of parameters characterizing wind loss and micro physics lead to equally good models. That means comparably good fits to the spectra can be found for different mass flow rates and truncation radii. But despite this ambiguity the fits showed, that there is a clear discrepancy of truncation radii in nuclei of elliptical galaxies of order $`10^3`$ to $`10^4r_\mathrm{S}`$ and those in low luminosity AGN (LLAGN) of only $`10^2r_\mathrm{S}`$ as pointed out by Quataert et al. (1999) (truncation radii are measured in Schwarzschild radius $`r_\mathrm{S}=2GM/c^2`$, $`M`$ mass of central black hole).

Since the transition to the ADAF can be understood as caused by evaporation of mass from the thin cool disk to a hot corona above, the truncation of the thin disk depends on the amount of mass flow there (Meyer & Meyer-Hofmeister 1994, Meyer et al. 2000b). If one measures the radii $`r`$ in Schwarzschild radius and the mass accretion rate $`\dot{M}`$ in Eddington mass accretion rates ($`\dot{M}_{\mathrm{Edd}}=40\pi GM/\kappa _{\mathrm{es}}c`$, $`\kappa _{\mathrm{es}}`$ electron scattering opacity) theory predicts a unique relation between truncation radius and mass accretion rate. But in fact the disks in LLAGN reach much closer in, to only 100 Schwarzschild radii, demanded by the observed UV radiation while the accretion rates $`\dot{M}/\dot{M}_{\mathrm{Edd}}`$ are comparable to those in the nuclei of elliptical galaxies.

In the following we suggest that the presence or absence of a disk dynamo explains the observed difference of truncation radii between LLAGN and elliptical galaxies with the same mass accretion rate. Sect. 2 gives a short description of the transition from disk accretion to a coronal flow. In Sect. 3 we investigate the presence or absence of disk dynamos in AGN. In Sect. 4 we present results of our computations including a magnetic field. A discussion follows in Sect. 5, conclusions in Sect. 6.

## 2 The transition from thin disk accretion to an ADAF caused by evaporation

In previous work (Meyer et al. 2000b) we presented a model for the corona above a geometrically thin standard disk around a black hole. This corona is fed by matter continuously evaporating from the cool disk. The accretion flow is thus divided into a hot coronal flow and a part remaining in the cool disk. The strength of the coronal flow increases inward so that at a certain distance in the inner region all matter is transferred to the coronal flow and the disk is truncated. The situation is the same for galactic and supermassive black holes. The theoretically derived relation for the dependence of evaporation rate on black hole mass and distance allows to determine the truncation radius for each accretion rate. This relation was successfully applied to X-ray transients. Liu & Meyer-Hofmeister (2001) discussed the application to AGN and found reasonable agreement, except that the truncation radii for LLAGN clearly were much smaller than predicted.

## 3 Disk dynamos in AGN

Up to now the evaluation of the transition from disk accretion to a coronal flow/ADAF has neglected a possibly important aspect: If the temperature in the underlying cool disk is high enough for dynamo action to occur the magnetic fields generated will penetrate also into the corona and affect the coronal accretion flow. Magnetic dynamos require sufficiently long Ohmic decay times compared to the dynamical time of the dynamo process. This ratio is measured by a magnetic Reynolds number, which is strongly temperature dependent for low temperatures.

Numerical simulations by Hawley et al. (1996) determined a critical value of the magnetic Reynolds number of about $`10^{3.5}`$. Below this value dynamo action becomes significantly suppressed. Gammie & Menou (1998) showed that dwarf nova accretion disks in quiescence are sufficiently cool to be below this critical value. AGN accretion disks have much larger physical dimensions. Menou and Quataert (2001) have demonstrated that the magnetic Reynolds numbers in such disks are above the critical value in the corresponding quiescent state. Here we are concerned with even further out ranges in elliptical galaxies where temperatures become extremely low.

In elliptical galaxies the black hole masses are about $`10^9M_{}`$ or more, in LLAGN only a few $`10^6M_{}`$. For higher black hole masses the midplane temperature is lower when normalized radii $`r/r_\mathrm{S}`$ and accretion rates $`\dot{M}/\dot{M}_{\mathrm{Edd}}`$ are kept the same. As an example we consider the situation for a mass flow rate $`10^{2.5}\dot{M}_{\mathrm{Edd}}`$ as derived for M81. For this rate we consider the disk structure at a radius of $`10^{3.5}r_\mathrm{S}`$ where the truncation would be expected from the standard curve (Fig.1).

Structure and evolution of accretion disks in AGN were investigated mainly in connection with the ionization instability in these disks (Mineshige & Shields 1990, Cannizzo & Reiff 1992, Cannizzo 1992, Siemiginowska et al. 1996). The radii of interest for our evaluation are larger than those considered there. For high central masses self-gravitation becomes important at such radii.

To explore a possible influence of disk dynamo action on the corona we refer to investigations for self-gravitating disks. Shlosman et al. (1990) discussed an attractive model of a “disk of clouds”: when the angular momentum transport is locally mediated gravitational interaction, the disk breaks up and possibly forms a layer of clouds (which should occur if the cooling time is short compared to the Kepler time). Recent numerical computations of Gammie (2001) for a thin horizontally extended layer confirm this picture with the formation of blobs. The clouds then by interaction have to transport the angular momentum, no effective dynamo would be expected. How can in this case the accretion flow via the clouds be transformed to a coronal flow ? According to the standard evaporation model these clouds must be embedded in a corona in equilibrium with the cool surfaces of the blobs. Thermal conduction of the hot corona to the cool surfaces establishes an equilibrium density and a mass flow rate from the blobs to the coronal gas, similar to the case of an underlying cool disk. The conductive flux requires a large enough surface area. Estimates give that a cloud covering factor $`C10^2`$ is required, which seems possible. Then the cloud layer is cut off and closer to the center only a coronal flow (ADAF) exists.

Balbus & Papaloizou (1999) discuss an alternative picture, in which waves propagating over long-range distances transport angular momentum and also energy (note that Gammie (2001) gives arguments against the formation of large scale coherence). The effective removal of angular momentum in this way might allow the existence of a disk marginally stable against break up. Because of non-local energy transport the disk would be cooler than a corresponding “high-$`\alpha `$” disk with the same angular momentum transport. For the latter disk (for a marginally stabilized disk see Paczyński 1978) in our case the formal values of $`\alpha `$ would be $`10`$. The midplane temperature in such a disk can be derived using a technique for solving the full vertically averaged disk equations similar to the accretion disk modeling in Cannizzo & Reiff (1992). We get for the temperature

$$T=10^{2.7}(\frac{M}{10^9M_{}})^{4/7}(\frac{r}{10^{3.5}r_\mathrm{S}})^{9/7}(\frac{\dot{M}}{10^{2.5}\dot{M}_{\mathrm{Edd}}})\mathrm{K}$$
(1)
and for the density

$$\rho =10^{12.6}(\frac{M}{10^9M_{}})^2(\frac{r}{10^{3.5}r_\mathrm{S}})^3\mathrm{gcm}^3.$$
(2)
The value of $`\alpha `$ is then dependent on $`M`$, $`r`$ and $`\dot{M}`$, incorporated in the formulae above. In the derivation of Eqs. (1) and (2) we used dust opacities, $`\kappa =1cm^2/g`$. The heating in the disk is due to internal dissipation.Direct and indirect (by scattering on the clouds) irradiation is negligible for the example M 87 because of the low luminosity.

With these numbers we can evaluate the magnetic Reynolds number

$$Re_\mathrm{M}V_\mathrm{s}H/\eta =\sqrt{2}\frac{\mathrm{}}{\mu }T/\mathrm{\Omega }^2\eta $$
(3)
with $`V_\mathrm{s}`$ isothermal sound velocity, $`H`$ disk scale height, $`\eta `$ magnetic diffusivity, $`\mathrm{}`$ gas constant, $`\mu `$ molecular weight, $`\mathrm{\Omega }`$ Kepler rotation frequency.

With the magnetic diffusivity $`\eta =m_ec^2\nu _{en}/n_ee^2`$ ($`n_e`$ electron number density, $`e`$ electron electric charge, $`m_e`$ electron mass, $`\nu _{en}`$ electron neutral collision frequency proportional to the number of neutral particles $`n_n`$) the Reynolds number becomes proportional to the ratio $`x_e=n_e/n_n`$ (Gammie & Menou 1998). This gives

$$Re_\mathrm{M}=10^{9.7}(\frac{M}{10^9M_{}})(\frac{r}{10^{3.5}r_\mathrm{S}})^3[X(T)T]^{1/2}$$
(4)
with $`X(T)=n_nx_e^2`$ is a function of temperature only determined by the Saha equation. The ionization state of the electron providing alkali metals becomes extremely low for the low temperatures. This yields a Reynolds number below the critical value of $`10^{3.5}`$ for $`T10^{2.9}`$ K. At the low densities of these disks ambipolar diffusion (Cowling 1976) increases the magnetic diffusivity (see also recent simulations by Sano & Stone 2002) and thereby reduces the magnetic Reynolds number. With typical dynamo values $`\beta 20`$ ($`\beta `$ ratio of gas pressure to magnetic pressure) the reduced $`Re_\mathrm{M}`$ is below the critical value for $`T10^{3.1}`$ K. Then no effect of disk magnetic fields is expected.

We conclude that in both these cases, a disk of clouds and a disk marginally stabilized by wave transport, the evaporation of cool gas into the corona occurs without the influence of a disk magnetic field.

For the disks in LLAGN we derived the temperature in the same way (such disks are not self-gravitating and we took $`\alpha `$=0.01). At the same (scaled) distance and mass accretion rate one obtains $`T=10^{3.3}`$ K. This is in the range already considered by Menou & Quataert (2001). In agreement with their analysis the magnetic Reynolds number becomes large enough to allow dynamo action.

## 4 The effect of the disk magnetic field on the truncation radius

We investigated how a magnetic field in the thin disk affects the coronal mass flow. The amount of matter that has evaporated from the disk and flows in the corona at a given distance from the black hole can be determined by solving a set of ordinary differential equations which describe vertical dynamical equilibrium, conservation of mass and energy and heat flux (Meyer et al. 2000b). Our earlier computations take magnetic fields of a coronal dynamo only implicitly, as a source of friction, into account as in standard disk theory. But an analytical estimate was given how $`\beta `$ enters into numerical results (Meyer et al. 2000b). For the present investigation we explicitly model the effect of additional magnetic fields (e.g. from an underlying disk) on coronal dynamical equilibrium, energy release, and diffusive inward mass flow. To a first approximation all effects can be accounted for if in all equations one formally replaces the gas pressure $`P`$ by the product $`P(1+1/\beta )`$. The result is a remarkable reduction of the strength of the coronal mass flow, i.e. of the evaporation efficiency, with increasing magnetic field.

Matter therefore remains in the cool disk down to much smaller radii before the coronal flow has picked up all the accretion flow and the disk is truncated. Fig. 1 shows the truncation radius as a function of the distance from the black hole. The values are scaled to $`\dot{M}_{\mathrm{Edd}}`$ and $`r_\mathrm{S}`$. (Inside of the truncation radius all matter flows in the form of a hot corona). The additional field shifts the truncation radius to smaller distances $`r`$. A moderate magnetic field strength in the corona ($`\beta `$=1) reduces the truncation radius by about a factor 5, a three times higher magnetic pressure results in a reduction by a factor 25. This is an important change.

In addition Fig. 1 shows data derived from observations by various authors. The observational data show a wide range of truncation radii. In particular only radii of $`10^4`$ or 300 $`r_\mathrm{S}`$ were considered in the investigation of nuclei of elliptical galaxies (Di Matteo et al. 2000) dealing with spectral fits based on an extension of the ADAF model including loss of mass and angular momentum by a wind. Note that spectral fits with wind loss yield only lower limits of truncation radii. For M 87 both, models with and without wind loss agree with the observed spectra. (Different results for NGC 4649 are due to different black hole mass.) Despite this spread in acceptable $`r`$ \- $`\dot{M}`$ combinations there is a clear discrepancy between the truncation radii for LLAGN and nuclei of ellipticals with comparable mass flow rates (Quataert et al. 1999).

## 5 Discussion

### 5.1 Dynamos in AGN disks

We expect dynamo action for LLAGN, but not for elliptical galaxies with a very high central mass. Otherwise if accretion rates are low it is possible that the temperature even in disks around $`10^6M_{}`$ is low enough to forbid dynamo action (an example could be Sgr A\* with a central mass of $`2.510^6M_{}`$ indicated by observations (Genzel et al. 1997) and an accretion rate from spectral fits based on ADAF models of $`10^4\dot{M}_{\mathrm{Edd}}`$ (Quataert & Narayan 1999) - if the thin disk really exists, doubts come from new Chandra observations (Narayan 2002)). On the other hand for higher mass accretion rates disks around the high mass black holes in elliptical galaxies are hot enough to allow dynamo action at the standard truncation radius, so that the true truncation radius becomes shifted inward to smaller radii.

### 5.2 The strength of the disk magnetic fields in the corona

To estimate the field strength of disk dynamo fields in the corona is difficult. The data for M 81 (Fig. 1) indicate a value of $`\beta 1/3`$ as appropriate. Near the truncation radius mass flow rates in corona and disk are of comparable size. This implies that the product $`(1+\frac{1}{\beta })HP`$ is about the same. From this one can estimate the required magnetic pressure in the corona $`(\frac{B^2}{8\pi })_{\mathrm{corona}}10^{1.2}(\frac{B^2}{8\pi })_{\mathrm{disk}}`$ where the temperature ratio $`T_{\mathrm{disk}}/T_{\mathrm{corona}}10^{4.6}`$ estimated for M 81 and $`(\beta )_{\mathrm{disk}}=20`$ were used. Such values may be reached if a corresponding fraction of the dynamo magnetic energy is cascaded to scales of the coronal scale height or larger (e.g. Arlt & Brandenburg 2001).

## 6 Conclusion

The investigation of the effect of magnetic fields of the underlying cool disk penetrating into the corona leads to interesting aspects for accretion in AGN. Our work suggests that the difference in truncation radii derived for LLAGN and ellipticals with similar accretion rates is due to the very different black hole masses, $`10^{6.5}M_{}`$ and $`10^9M_{}`$ respectively. For the same mass flow rate and at the same distance from the black hole (when measured in units of $`\dot{M}_{\mathrm{Edd}}`$ and $`r_\mathrm{S}`$) the disks in LLAGN are hot enough, but those in ellipticals are cool and self-gravitating, no magnetic dynamo work. In LLAGN the magnetic fields affect the coronal flow and shift the truncation to much smaller radii. The truncation of disks in X-ray binaries and the spectral transitions from hard to soft state (Meyer et al. 2000a) will also be affected when disk dynamos occur.