text:
# A Magnetohydrodynamic Nonradiative Accretion Flow in Three Dimensions

## 1. Introduction

Observations of underluminous accretion-powered X-ray sources pose a difficulty for accretion theory. The Galactic center is a particularly compelling example: despite the inferred presence of a 2 million solar mass black hole, and the availability of a substantial gas reservoir, only a small fraction of the expected radiation is detected. In a standard Keplerian accretion disk, liberated mechanical energy is promptly radiated, and there is a direct relationship between the accretion rate $`\dot{M}`$ and the luminosity $`L`$. A number of alternatives to the standard disk have been put forward, which are classified by the ultimate repository of the orbital energy. In Advection Dominated Accretion Flows (ADAFs, e.g., Narayan & Yi 1994), the energy is retained by the plasma and advected into the hole. In contrast, Blandford & Begelman (1999) proposed that the bulk of the accreting gas and energy might be carried off by a wind. A more recent variation, referred to as Convective Dominated Accretion Flow (CDAF; e.g., Quataert & Gruzinov 2000), brings enhanced viscous heating to the fore, leading to a strong convective flow that stifles the accretion. Since the common feature of all these models is negligible radiative losses, we shall refer to them generically as Non-Radiative Accretion Flows (NRAFs).

Earlier simulations of NRAFs have used several simplifying approximations. Most (e.g., Igumenshchev & Abramowicz 1999, 2000, hereafter IA99, IA00; Stone, Pringle, & Begelman 1999, hereafter SPB) have assumed the flow evolves due to a fortified kinematic viscosity, $`\nu `$. While illuminating some aspects of NRAFs, there are substantial drawbacks to this approach. Most importantly, the results strongly depend upon the specific recipe adopted for $`\nu `$. The origin of angular momentum transport is not unknown, however. It arises from the magnetorotational instability (MRI; Balbus & Hawley 1991), and hydrodynamic and magnetohydrodynamic (MHD) flows have fundamentally different properties.

Stone & Pringle (2001; hereafter SP) include MHD in their axisymmetric simulations. The stress required to drive the accretion flow emerges self-consistently; it is not a free parameter. But the two-dimensional (2D) restriction is a significant limitation. First, the anti-dynamo theorem (e.g., Moffat 1978) prevents the indefinite maintenance of the poloidal magnetic field in the face of dissipation. Second, axisymmetric simulations tend to over-emphasize the “channel” mode (Hawley & Balbus 1992) which produces coherent internal magnetized flows rather than the more generic MHD turbulence. Consequently, a fully self-consistent NRAF model requires three-dimensional (3D) MHD. In this Letter we present and discuss the results of a prototype 3D NRAF simulation.

## 2. The Simulation

The simulation evolves the equations of ideal MHD, i.e., the continuity equation, the induction equation, the Euler equation,

$$\rho \frac{𝐯}{t}+(\rho 𝐯)𝐯=\left(P+𝒬+\frac{B^2}{8\pi }\right)\rho \mathrm{\Phi }+\left(\frac{𝐁}{4\pi }\right)𝐁,$$
(1)
and an internal energy equation

$$\frac{\rho ϵ}{t}+(\rho ϵ𝐯)=(P+𝒬)𝐯.$$
(2)
The variables have their usual meanings; $`𝒬`$ is an explicit artificial viscosity (Stone & Norman 1992a), and $`\mathrm{\Phi }=GM/(rr_g)`$ is the Paczyński & Wiita (1980) pseudo-Newtonian potential. We use units with $`GM=r_g=1`$. The equation of state is adiabatic, $`P=\rho ϵ(\mathrm{\Gamma }1)`$, with $`\mathrm{\Gamma }=5/3`$. Radiation transport and losses are omitted. Since there is no explicit resistivity or shear viscosity, the gas can heat only by adiabatic compression, or by the artificial viscosity $`𝒬`$. We evolve the equations using time-explicit Eulerian finite differencing with the ZEUS algorithms (Stone & Norman 1992a; 1992b; Hawley & Stone 1995; Hawley 2000).

The major difficulties of a 3D simulation include the large number of grid zones needed to resolve an extended spatial domain, and the number of time steps required to evolve the problem over several orbital periods at large radius. In this simulation, we use a somewhat restricted resolution of $`128\times 32\times 128`$ grid zones in cylindrical coordinates $`(R,\varphi ,Z)`$. The radial grid extends from $`R=1.5`$ to $`R=170`$. There are 36 equally-spaced zones inside $`R=15`$, and 92 zones outside this point that increase logarithmically in size. Similarly, 50 $`Z`$ zones are equally spaced between $`10`$ and $`10`$ with the remainder of the zones logarithmically stretched to the $`Z`$ boundaries at $`\pm 60`$. The azimuthal domain is limited to $`\pi /2`$. Hawley (2001) found that reduced angular coverage preserves qualitative features of a simulation with only a small reduction in magnetic energy and stress. The radial and vertical boundary conditions are simple zero-gradient outflow conditions; no flow into the computational domain is permitted. The $`\varphi `$ boundary is periodic. The magnetic field boundary condition is set by requiring the transverse components of the field to be zero outside the computational domain, while the perpendicular component satisfies the divergence-free constraint.

The prototype simulation starts with of a constant specific angular momentum ($`\mathrm{}`$) torus with a pressure maximum at $`r=100`$ (similar to model F of SP) and an inner edge at $`r=75`$. The initial magnetic field consists of poloidal loops lying along isodensity contours with a volume-averaged energy of $`\beta =P_{gas}/P_{mag}=200`$. The average specific energy of the magnetofluid in the torus is $`4.6\times 10^3`$ (in units where $`GM=r_g=1`$), i.e., it is bound. The simulation is run out to 5 orbits at the initial pressure maximum. This corresponds to 1429 orbits at the marginally stable orbit, $`R=3`$. The simulation requires over 1.7 million timesteps with an average $`\mathrm{\Delta }t=0.0152`$.

## 3. Results

During the first orbit the field in the torus is amplified by the MRI and by shear, causing the torus to expand. The MRI acts most visibly near the equatorial plane, where the field is predominantly vertical, and the long-wavelength, nearly axisymmetric modes of the MRI grow rapidly. The linear growth phase ends shortly after one orbit, and the magnetic energy has increased to $`\beta 2`$–10.

Over the next orbit, low-$`m`$ spiral arms of gas accrete from the inner edge of the torus, forming a vertically thin and very nonuniform disk. Strong magnetic fields surround the gas. Due to the initial field topology, a current sheet forms near the equator which proves unstable to vertical oscillations. As more gas leaves the torus, the disk fills out and thickens. Low density material is driven off the forming disk to create a backflow.

As time advances, the initial constant-$`\mathrm{}`$ torus evolves to a nearly Keplerian angular momentum distribution. There is net accretion inside of $`R=100`$, and net outflow beyond this point. Inside of $`R=100`$ the NRAF forms a modestly thick, nearly Keplerian disk. In the disk, accretion results from MHD turbulence, but the final accretion rate into the black hole need not match the supply rate precisely, and it does not. Because the inflow is hot, geometrically thick, and nearly Keplerian, inflow is difficult past the marginally stable orbit. To accrete the gas must pass through a narrow gap in the centrifugal barrier at the equator. Gas splashes off the centrifugal barrier near the hole, forming a pressure-supported torus and creating a backflow that adds to a growing low density envelope around a higher density Keplerian core.

Between orbits 4 and 5, the disk is well-formed and reasonably steady. Figure 1 is a series of $`(R,z)`$ contour plots of azimuthally averaged values at $`t=5`$ orbits showing that most of the gas is located in a modestly thick disk surrounded by a low density, highly magnetized atmosphere. Near the equator gas pressure is dominant, i.e., $`\beta >1`$, but in the surrounding region $`\beta <1`$. The gas pressure scale height is $`H7`$–10, decreasing rapidly inside $`R=10`$. The total pressure (gas plus magnetic) is much smoother, and has a scale height $`H20`$ at $`R=100`$ that decreases slowly inward. The gas density, vertically averaged over one gas pressure scale height, is nearly constant from $`R=50`$ down to $`R=10`$, where it increases rapidly to a peak at $`R=5.5`$. Inside $`R=10`$, the disk resembles a pressure supported torus.

The specific angular momentum $`l(R)`$ tracks a Keplerian distribution throughout the disk, although in places it is up to 10% below the Keplerian value. The angular velocity $`\mathrm{\Omega }`$ is constant along cylinders through the main portion of the disk, consistent with $`P=P(\rho )`$ there. In hydrodynamic simulations, surfaces of constant entropy $`S`$ coincide with those of specific angular momentum $`\mathrm{}`$. This corresponds to marginal stability by one of the Høiland criteria. As already noted by SP, with MHD the level surfaces of $`S`$ and $`\mathrm{}`$ do not line up, particularly inside the main disk.

Magnetic fields in the disk produce Maxwell stress, MHD turbulence, and angular momentum transport. Between orbit 4 and 5 the ratio of the vertical- and $`\varphi `$-averaged Maxwell stress to the gas pressure ranges from $`0.1`$ to 0.2 inside of $`R=100`$.

The mass inflow and outflow rates through every cylindrical radius are computed as a function of time. An average over the last orbit shows that these two are nearly equal, although inflow exceeds outflow (Fig. 2). Of course the instantaneous inflow and outflow rates merely represent the flow produced by the MHD turbulent velocity fluctuations. The net accretion rate, $`\dot{M}`$, is due to the resulting drift velocity, which is always considerably smaller.

At the end of the run the total mass on the grid has decreased by 5.4%. Roughly 51% of this leaves through the outer radial boundary, 43% through the upper and lower $`Z`$ boundaries, and the remaining 6% is accreted into the black hole. A higher resolution torus simulation (Hawley & Krolik 2001) shows that the magnetic stress can remain large down to and beyond the marginally stable orbit. This effect could increase $`\dot{M}`$, but the present simulation is not sufficiently well-resolved to address this point. The flow through the outer radial boundary is a consequence of the gain in $`\mathrm{}`$ in the outer part of the torus. The flow through the $`Z`$ boundaries is driven from the accretion disk by gas and magnetic pressure. Part of this is an unbound, high temperature hollow conical outflow confined to the axis region by surrounding magnetic pressure. This can be seen in the gas pressure plot of Figure 1. The remainder of the outflow remains bound, although its specific energy is greater than that of the initial torus.

The initial average specific energy for the gas is $`4.6\times 10^3`$. At 5 orbits this has decreased to $`4.9\times 10^3`$ as the gas remaining on the grid becomes more bound. The total thermal, magnetic, poloidal kinetic, and orbital energies of the gas that remains on the grid all increase with time. The magnetic energy has increased the most, receiving 42% of the total net energy increase; 95% of this is in toroidal field. The thermal, orbital and kinetic energies comprise 25%, 21%, and 12% of the increase respectively.

## 4. Discussion

### 4.1. Comparison with $`\alpha `$ models

Full 3D MHD simulations are difficult and expensive; are they necessary? With respect to the hydrodynamical alternative the answer is clearly yes. Such simulations (IA99; SPB; IA00) have shown that the results depend strongly on the magnitude of $`\alpha `$ and the form assumed for the stress. Thus, a self-consistent stress is essential for understanding NRAF dynamics.

One should not lose sight of the fact that this is a high Reynolds number turbulent system, not a low Reynolds number laminar flow. For example, the result that large Shakura-Sunyaev $`\alpha `$ ($`\nu =\alpha c_s^2/\mathrm{\Omega }`$) accretion flows are stable, laminar and accrete fully into the central hole (IA99; IA00) is a literal consequence of using a Navier-Stokes viscosity for the stress. With stress from MHD turbulence $`\alpha 1`$ necessarily implies velocity fluctuations with speeds comparable to the sound speed $`c_s`$ on scales of order the size of the system.

A second hydrodynamic result is that NRAFs are driven to a state of marginally convective instability (SPB; Quataert & Gruzinov 2000; IA00). This is the basis for CDAF solutions. But this cannot hold in MHD. Recall the form of the Høiland criteria that applies to adiabatic perturbations in magnetized flow (Balbus 1995):

$$N^2+\frac{\mathrm{\Omega }^2}{\mathrm{ln}R}>0$$
(3)
$$\left(\frac{P}{Z}\right)\left(\frac{\mathrm{\Omega }^2}{R}\frac{\mathrm{ln}P\rho ^{5/3}}{Z}\frac{\mathrm{\Omega }^2}{Z}\frac{\mathrm{ln}P\rho ^{5/3}}{R}\right)>0$$
(4)
where $`N^2=(3/5\rho )P\mathrm{ln}P\rho ^{5/3}`$ is the Brunt-Väisälä frequency. Note the symmetry between the angular velocity gradients and the entropy gradients, and the absence of any stabilizing epicyclic frequency. At least one of these criteria must always be strongly, not marginally, violated in an accretion flow: this is what is responsible for the existence of the turbulent stress. Convection does not arise in addition to this fundamental instability, it is part of the instability that produces the anomalous stress in the first place. In particular, one cannot argue for the marginal stability of a rotationally modified convective process independently of what is called $`\alpha `$. Whatever the sign of the Brunt-Väisälä frequency, the instability is essentially the same one at work in a standard Keplerian disk. Indeed, in contrast to hydrodynamical treatments, we find that the MHD flows show no tendency at all toward marginal stability.

The 3D MHD simulation does show large scale turbulent rolls and outflows qualitatively similar to those reported in small $`\alpha `$ hydrodynamic simulations. In the hydrodynamic simulations, however, all the rotational energy released by the stress goes immediately into heat. Some of this heat returns to kinetic energy in the form of large-scale convection. This implies that large-scale shear is converted into heat, producing large-scale convective rolls, which then prevent the accretion process that was driving the heating in the first place. In MHD, the orbital energy goes directly into turbulent kinetic and magnetic energies. The magnetic energy constitutes a particularly important component of the total energy budget. Heating results from dissipation at small scales. In a real plasma, this will depend on microscopic resistivity and viscosity. In the simulation, turbulent rolls are not thermal in origin. More generally the kinetic energy of convective motions is small compared with the dynamical activity associated with MHD instabilities and with the large scale flow itself.

We also remark that in any turbulent flow, the net accretion rate will always be small compared to any locally sampled instantaneous value. As Figure 2 shows, inward and outward accretion rates almost exactly cancel in our purely dynamical simulation. This is not a unique property of inwardly-transported convective energy fluxes. It is a hallmark of turbulence. In both viscous hydrodynamical and MHD treatments, turbulent eddy velocities are much larger than the residual inward drift velocity. It is precisely this ordering that allows the $`\alpha `$ formalism for the mean flow dynamics to be derived from the equations of motion themselves (Balbus & Papaloizou 1999).

### 4.2. Comparison with 2D MHD models

Next we compare the present simulation with the 2D models of SP. In the equivalent run (SP’s Run F), the initial infall phase during orbits 2–3 was relatively smooth and dominated by the channel flow of the MRI. Subsequent evolution was turbulent. In 3D, the initial smooth phase does not occur; the evolution is always turbulent. The 3D flow also generates greater outflow from the disk.

After 3 orbits in both 2D and 3D, the inner regions of the flow are dominated by MHD turbulence driven by the MRI. Interestingly, there appear to be few differences between time- and angle-averaged quantities in either case. For example, comparison of the contours in Fig. 1 with Figure 4 in SP show that both produce a nearly barytropic disk near the midplane with $`\mathrm{\Omega }`$ contours parallel to cylindrical radii, but with no correlation between the specific angular momentum and entropy. Moreover, the net $`\dot{M}`$ values are similar: $`3\times 10^3`$ in 2D versus $`2\times 10^3`$ in 3D, in units of the initial mass of the torus per orbit at the initial pressure maximum. Radial profiles of time-averaged data in the main disk show many similarities: in both cases $`P_{gas}10P_{mag}`$ and $`v_\varphi >c_s>v_{Alfven}`$ at all radii (see figure 6 in SP). The most obvious differences between the 2D and 3D models are near the axis, and these are attributable to spherical 2D versus cylindrical 3D inner boundary geometries.

Eventually all axisymmetric models are limited by the antidynamo theorem. Toward the end of the SP simulations the turbulence in the initial torus dies down, and remains only in the flow close to the black hole. Higher resolution is possible in 2D, but this only prolongs the duration of the turbulent phase before eventual decline. In addition, the toroidal field MRI cannot be simulated in axisymmetry, which limits 2D models to initial conditions that contain large-scale poloidal field.

## 5. Conclusion

We have performed a three-dimensional MHD simulation of an NRAF originating from a constant-$`\mathrm{}`$ torus located at 100 gravitational radii. The resulting flow does not resemble a classical ADAF. Most of the infalling mass and released energy does not end in the hole. In this sense the flow might be better described as an inflow/outflow solution of the type outlined by Blandford & Begelman (1999). The general appearance of the flow is similar to the low-$`\alpha `$ hydrodynamic models of SPB, IA99 and IA00. Of course this is what one should expect for any hot, rotating, nearly Keplerian accretion flow. We have argued that the CDAF picture is also inappropriate because the the classic Høilland criteria do not apply to an MHD flow. The outflow is not driven by thermal convection, but by the action of the same instability that accounts for $`\alpha `$, namely the MRI. Marginally stability is not attained.

We conclude by considering some directions for future work. In the present simulation an internal energy equation is used and the only source of dissipative heating is the artificial viscosity $`𝒬`$. Some amount of the total energy is inevitably lost in this formulation. A portion of this can be recaptured through artificial resistivity, although this made no significant difference in flow structure or accretion rates when tested in 2D (SP). In any case, a more detailed study of the energy flow in a 3D simulation would be enlightening. Higher resolution is needed in the inner region to capture the details of the inner disk boundary and the flow past the marginally stable orbit. Models with a greater radial domain would allow the inflow to proceed over several decades in distance and to liberate greater amounts of gravitational energy. High-resolution axisymmetric MHD simulations can play an important role here in mapping out the types of flows that are possible. Longer evolutions would provide greater confidence that the results have lost memory of the initial conditions, but, because of the antidynamo theorem, true steady states can only be obtained in 3D.

JFH and SAB acknowledge support under NSF grant AST-0070979, and NASA grants NAG5-9266, NAG5-7500, and NAG5-106555. JMS acknowledges support under NASA grant NAG5-3836. The simulation was carried out on the IBM Bluehorizon system of the San Diego Supercomputer Center of the National Partnership for Advanced Computational Infrastructure, funded by the NSF.


format: text

text:
# Evaluating the RiskMetrics Methodology in Measuring Volatility and Value-at-Risk in Financial Markets

## 1 Introduction

Risk management is one of the top priorities in the financial industry today. A huge effort is being invested into developing reliable risk measurement methods and sound risk management techniques by academics and practitioners alike. J. P. Morgan was the first to develop a comprehensive market risk management methodology based on the Value-at-Risk (VaR) concept . Their product, RiskMetrics  has become extremely popular and widespread. It has greatly contributed to the dissemination of a number of basic statistical risk measurement methods and to the general acceptance of VaR as an industry standard. Although backtests performed first by J. P. Morgan and later by other market participants lent support to the RiskMetrics model, its basic assumptions were shown to be questionable from several points of view (see e.g. ). Moreover, the existence of fat tails in real market data (see  for a discussion) is in a clear conflict with RiskMetrics’ assumption of normally distributed returns, which can lead to a gross underestimation of risk. Furthermore, serious doubt has recently been raised as to the stability and information content of the empirical covariance matrices used by the model for calculating the risk of portfolios .

Accurate evaluation of risks in financial markets is crucial for the proper assessment and efficient mitigation of risk. Therefore, it is important to see to what extent widespread methodologies like RiskMetrics are reliable and what their possible limitations are. In particular, we try to understand why, despite the evident oversimplifications embedded in the model, it can perform satisfactorily. We will argue that the apparent success of RiskMetrics is due basically to the way risk is quantified in this framework, which does not necessarily mean that this particular risk measure is the most adequate one.

The paper is organized as follows. In Section 2 a sketchy overview of the RiskMetrics methodology will be presented. In Section 3 the reasons for the success of RiskMetrics will be discussed. Finally, Section 4 is a short summary.

## 2 Overview of the RiskMetrics Methodology

It is well-known that daily returns are uncorrelated whereas the squared returns are strongly autocorrelated. As a consequence, periods of persistent high volatility are followed by periods of persistent low volatility, a phenomenon known as “volatility clustering”. These features are incorporated in RiskMetrics by choosing a particular autoregressive moving average process to model the price process (see below). Furthermore, RiskMetrics makes the very strong assumption that returns are conditionally111Conditional here means conditional on the information set at time $`t`$, which usually consists of the past return series available at time $`t`$. normally distributed. Since the standard deviation of returns is usually much higher than the mean, the latter is neglected in the model222The typical yearly mean return in equity markets is 5%, while the typical standard deviation is 15%. For a time horizon of one day, however, the mean becomes 5%/250=0.02%, while the standard deviation becomes $`15\%/\sqrt{250}0.95\%`$, i.e. about 50 times larger than the mean., and, as a consequence, the standard deviation remains the only parameter of the conditional probability distribution function. In order to avoid the usual problems related to the uniformly weighted moving averages, Riskmetrics uses the so called exponentially weighted moving average (EWMA) method \[2, pp. 78–84\] which is meant to represent the finite memory of the market. Accordingly, the estimator for the volatility is chosen to be

$$\sigma _{t+1|t}^2=\frac{_{\tau =0}^{\mathrm{}}\lambda ^\tau r_{t\tau }^2}{_{\tau =0}^{\mathrm{}}\lambda ^\tau }=(1\lambda )\underset{\tau =0}{\overset{\mathrm{}}{}}\lambda ^\tau r_{t\tau }^2,$$
(1)
where $`\lambda `$ is a parameter of the model ($`0<\lambda <1`$). The notation $`\sigma _{t+1|t}`$ emphasizes that the volatility estimated on a given day ($`t`$) is actually used as a predictor for the volatility of the next day ($`t+1`$). The daily VaR at confidence level $`p`$ (e.g. 95%) can then be calculated (under the normality assumption) by multiplying $`\sigma _{t+1|t}`$ with the $`1p`$ quantile of the standard normal distribution. Moreover, this technique can be used to measure the risk of individual assets and portfolios of assets as well. For linear portfolios (i.e. containing no options) the usual method to obtain the volatility is to estimate the covariance matrix of asset returns, element-by-element, using the EWMA technique and then calculate the portfolio volatility as $`\sigma _p^2=_{i,j}w_iw_j\sigma _{ij}`$, where $`w_i`$ is the weight of asset $`i`$ in the portfolio. Equivalently, however, one can calculate the return on the portfolio first, and then apply the EWMA technique directly to the whole portfolio . Finally, the value of the parameter $`\lambda `$ is determined by an optimization procedure. On a widely diversified international portfolio, RiskMetrics found that the value $`\lambda =0.94`$ produces the best backtesting results \[2, pp. 97–101\].

## 3 Why Does RiskMetrics Work?

In order to explain the reasonably successful performance of RiskMetrics, first we recall the work by Nelson  who showed that even misspecified models can estimate volatility rather accurately. More explicitly, in  it is shown that if the return generating process is well approximated by a diffusion, a broad class of even misspecified ARCH models333 See  for a survey on ARCH type models (e.g. ARCH, GARCH, IGARCH, etc.). can provide consistent estimates of the conditional volatility. Since RiskMetrics can be considered as an IGARCH(1,1) model444 In the IGARCH(1,1) model, returns are generated by the following stochastic process: $`r_t=\sigma _t\epsilon _t`$, where $`\epsilon _t`$ i.i.d.(0,1), $`\sigma _t^2=\omega +\beta _1\sigma _{t1}^2+\alpha _1\epsilon _{t1}^2`$ and $`\alpha _1+\beta _1=1`$. Since from Eq. (1) $`\sigma _{t+1|t}^2=\left(1\lambda \right)\left(r_t^2+\lambda \frac{\sigma _{t|t1}^2}{1\lambda }\right)`$, it follows that $`\sigma _{t+1|t}^2=\lambda \sigma _{t|t1}^2+(1\lambda )r_t^2`$ just as in the IGARCH(1,1) model ($`\omega =0`$, $`\beta _1=\lambda `$, $`\alpha _1=1\lambda `$). the results of  offer a natural explanation for the success of RiskMetrics in estimating volatility. Actually, in the RiskMetrics framework, this estimate is used as a one-day ahead volatility forecast, nevertheless it seems that this does not significantly worsen its accuracy. However, if one uses this estimate to calculate (as often required by regulators) a multiperiod forecast using the simple “square-root-of-time” rule555 The rule says that $`\sigma _{t+h|t}=\sqrt{h}\sigma _{t+1|t}`$, which is derived on the basis of the assumption of uncorrelated returns \[2, pp. 84–88\]., the quality of the forecast is bound to decline with the number of periods. A comparative study of the rate of deterioration of these forecasts with time within the RiskMetrics model resp. other, more sophisticated volatility models is an interesting topic for further research.

Let us turn to the apparent success of RiskMetrics in estimating the VaR now. The model sets the confidence level at 95%. The prescription to obtain this 5% quantile is to simply multiply the volatility estimate by 1.65 (as if returns were conditionally normally distributed). Why such a recipe can pass the backtests can be understood if one analyzes numerically the 5% quantiles of a few leptokurtic (fat tailed) distributions. It is very often found that despite the presence of fat tails, for many distributions the 5% quantile is roughly -1.65 times the standard deviation. For example, the 5% quantile of the Student $`t`$ distribution with 7 degrees of freedom (which is leptokurtic and has a kurtosis of 5 similar to the typical kurtosis of returns in financial markets) is -1.60, very close to -1.65, or, conversely, the -1.65 percentile is 4.6%. This is illustrated in Fig. 1, where the PDF’s of the Student and the normal distributions are shown. It can be seen from the figure that the areas under the two PDF’s to the left of the -1.65 line are roughly equal. We also analyzed the empirical frequencies of RiskMetrics 95% VaR violations which correspond to returns $`r_{t+1}<1.65\sigma _{t+1|t}`$, where $`\sigma _{t+1|t}`$ is the RiskMetrics volatility estimate obtained on day $`t`$. For the 30 stocks of the Dow Jones Industrial Average (DJIA), which are among the most liquid stocks traded on the New York Stock Exchange, it was found666 The sample consisted of the daily closing prices of the 30 stocks of the DJIA from August 1996 to August 2000, about 1000 data points for each stock. that the VaR violations frequency was $`4.7\pm 0.6\%`$, i.e. close to 5%. This explains why RiskMetrics is usually found so successful in evaluating risk (which it defines as the VaR at 95% confidence).

It is evident, however, that for higher significance levels (e.g. 99%) the effect of fat tails becomes much stronger, and therefore the VaR will be seriously underestimated if one assumes normality. For example, the 1% quantile of the Student $`t`$ distribution considered above is -2.54, significantly larger than under the normality assumption (-2.33), while the percentile corresponding to -2.33 is 1.43%. This can also be seen from Fig. 1. Furthermore, for the DJIA data considered above, the rejection frequencies were $`1.4\pm 0.3\%`$, significantly larger than 1%. These effects are even more pronounced for a truncated Lévy distribution (TLD). For example, a Koponen-like  TLD with Lévy exponent $`\nu =1.50`$, scale parameter $`c=0.50`$ and truncation parameter $`\lambda =0.17`$, which provides an excellent fit to the Budapest Stock Exchange Index (BUX) data, has a 5% quantile equal to $`1.76\sigma `$ whereas the 1% quantile is already $`3.75\sigma `$!

Therefore, it can be concluded that the satisfactory performance of RiskMetrics in estimating VaR is mainly the artifact of the choice of the significance level of 95%. However, existing capital adequacy regulations require 99% confidence, and at this level RiskMetrics systematically underestimates risk.

## 4 Summary

In this paper we analyzed the performance of RiskMetrics, perhaps the most widely used methodology for measuring market risk. The Riskmetrics model is based on the unrealistic assumption of normally distributed returns, and completely ignores the presence of fat tails in the probability distribution, a most important feature of financial data. For this reason, one would expect the model to seriously underestimate risk. However, it was commonly found by market participants that RiskMetrics performed satisfactorily well and this helped the method to quickly become a standard in risk measurement. Nevertheless, we found that the success of RiskMetrics is actually the artifact of the choice of the risk measure: the effect of fat tails is minor when one calculates Value-at-Risk at 95%, however, for higher significance levels fat tails in the distribution of returns will make the simple RiskMetrics rule of calculating VaR to underestimate risk.

RiskMetrics has played and continues to play an extremely useful role in disseminating risk management ideas and techniques, even if oversimplified. It is available free of charge and, coming with a careful documentation, it is completely transparent and amenable for a study like the present one: its limitations can be explored and, given sufficient resources, overcome. This is far from being the case with the overwhelming majority of the commercially available risk management systems which incorporate at least as strong simplifications as RiskMetrics, but coming in the the form of “black boxes,” are completely impossible to modify. The continuing dominance of the Gaussian paradigm in risk management software packages represents an industry-wide model risk.

## Acknowledgements

It is a pleasure to thank B. Janecskó for useful interactions. This work has been supported by the Hungarian National Science Found OTKA Grant No. T 034835.


text:
# Stability of Accretion Disks in Presence of Nucleosynthesis

## 1 Introduction

Black hole accretion flow can reach a temperature of about $`T_{virial}=\frac{1}{k}\frac{GM_{BH}m_p}{r_g}5.2\times 10^{12}`$K if no cooling processes are included. Here, $`r_g=2GM_{BH}/c^2`$ is the gravitational radius of a black hole, $`k`$ is the Boltzmann constant, $`G`$ is the Gravitational constant, $`M_{BH}`$ is the mass of the black hole, $`m_p`$ is the mass of the proton and $`c`$ is the velocity of light. In a realistic flow, the temperature is smaller. For instance, in a Shakura-Sunyaev (1973) Keplerian disk, the temperature becomes close to a few times $`10^7`$K. In a two-temperature transonic flow, the ions may remain hot ($`T_p10^{811}`$K) while the electrons may be cooler (($`T_e10^{79}`$K) depending on accretion rate of the Keplerian and sub-Keplerian components (Colpi, Maraschi & Treves, 1984; Chakrabarti & Titarchuk, 1995). Since temperature inside a star is much less compared to the ion temperature mentioned above, considerable nucleosynthesis is expected inside an accretion disk. It is true that the free-fall timescale $`t_{ff}rv_{ff}(\frac{r}{r_g})^{3/2}\frac{2GM_{BH}}{c^3}`$ s is very small compared to the stellar age and density of gas inside an accretion disk, $`\rho _d\frac{\dot{M}_{Ed}t_g}{r_g^3}(\frac{r}{r_g})^{3/2}`$gm cm<sup>-3</sup> is also very small compared to the stellar density but the reaction cross-sections are often very temperature sensitive ($`T^{416}`$) and the higher temperature could compensate for the lower density.

Considerable work has been done in the past on nucleosynthesis inside a disk taken into consideration. Paczyński and Jaroszyński (1978) have studied the importance of nucleosynthesis at high temperature thin disks made up of Helium and Carbon. In these disks, electron degeneracy pressure supports gravity. Taam & Fryxell (1985) studied thin accretion disks of hydrogen-rich matter and showed that under certain condition, thermal instability of a Keplerian disk may totally disappear. They consider very high density disks where degenerate pressure in the mid-plane supports the disk structure. Chakrabarti, Jin & Arnett (1987) and Jin, Arnett & Chakrabarti (1989) studied nucleosynthesis in thick accretion disks of realistic density and temperature. These disks have very little radial velocity and a significant amount of nucleosynthesis was possible only if viscosity parameter of Shakura-Sunyaev (1973) is very low $`\alpha _{ss}\stackrel{<}{}10^4`$. They showed that in order that thin disks remain stable against convection, only $`pp`$ reactions are allowed in these disks in the radiation pressure supported regions. However, for hot thick disks, even CNO cycle would be allowed without destabilizing the disk. These results were later repeated by independent groups (Arai & Hashimoto, 1992; Hashimoto et al. 1993). More recently, Chakrabarti & Mukhopadhyay (1999) and Mukhopadhyay & Chakrabarti (2000) revisited the problem with a self-consistent accretion disk solution which includes heating and cooling. They find that nucleosynthesis in some range of the parameter space could be high enough to have an energy release or absorption by nuclear reaction comparable to the viscous heating. These studies are finding attention in recent years since one could detect lines from various isotopes from jets and companions and can, in principle, have some estimate of the conditions around a black hole.

However, work mentioned above was purely numerical in nature and was mainly due to repeated use of two types of analysis: First, the steady state thermodynamic condition of a disk is determined by using fourth order Runge-Kutta method without nucleosynthesis taken into account. Using temperature and density distributions thus obtained at each point, nucleosynthesis code is run which mainly involves solving 255 linear equations simultaneously, one equation for each isotope. For each infinitesimal advancement of matter toward the black hole, nuclear composition was updated and energy release or absorption was computed locally. After one complete run, the energy release is incorporated into the hydrodynamic equations and the integration scheme using the Runge-Kutta method is carried out to have updated thermodynamic conditions. This process is continued till convergence. In the present paper, our approach is completely different. We attack the problem analytically and estimate effects of nucleosynthesis by using sonic point analysis. This way the problem becomes fully self-consistent to begin with and no iterative step is needed. Presence of a saddle type sonic point is a necessary condition for a steady solution, but it is not sufficient. The flow must connect a black hole horizon with infinity. Of course, in order to have the problem under control, we chose only a few reactions which are likely to be very important, namely, the Proton-Proton and photo-dissociation of deuterium and helium. Since we are interested in stability, we considered only those reactions (both exothermic and endothermic) which change specific energy (either way) significantly. We, however, neglect viscous heating or possible radiative transfer effects.

In any successful model of the black hole accretion, a source of the soft photon and a hot, Comptonising cloud must be present. In our model, this Comptonising cloud is the inner part of the disk itself. For concreteness, we concentrate on this model of black hole accretion, namely, one in which the flow close to a black hole is transonic (this is a necessary condition for any model), but far away, it is sub-divided into Keplerian and sub-Keplerian components. According to the properties of the transonic disks (see, Chakrabarti & Titarchuk, 1995; and references therein) higher viscosity flow settles into the equatorial plane, while quasi-spherical flows with lower viscosity surround the Keplerian disk. Soft photons emitted from the Keplerian component are reprocessed by the electrons at the centrifugal pressure supported boundary layer and hard X-rays are emitted. Recent observational results strongly support the presence of Comptonising clouds close to a black hole (Mirabel & Rodriguez, 1999; Homan et al., 2000) which have properties very similar to our centrifugal pressure supported regions in the inner accretion disks.

We had to choose a particular model in order to compute the optical depth of matter through which soft photons are Comptonized. Since at high temperature, photo-dissociation is important one has to compute number of hard photons present in this inner disk. Ordinarily, photo-dissociation at a high temperature is taken for granted, but in the present situation, there are paucity of very hot photons since only few of them are emitted from the Keplerian component and even fewer are intercepted by the sub-Keplerian flow located ahead of the Keplerian disk. What is more, optical depth may be so low that most of the time these photons do not interact adequately with deuterium and helium. Because of this, we compute photon number density self-consistently from the emitted spectrum of a sub-Keplerian flow.

In the next Section, we present the basic hydrodynamic equations which govern the steady state flow and perform the sonic point analysis. We show that effect of the nuclear term is significant in changing topological properties of the flow. In §3, we show that number of sonic points is larger in presence of nucleosynthesis. We present a complete global solution and show how it is modified in presence of nucleosynthesis. Finally, in §4, we draw our conclusions.

## 2 Model Equations

In what follows, we use the unit of distance (x) to be $`r_g=2GM_{BH}/c^2`$, unit of velocity to be $`c`$ and the unit of mass to be $`M_{BH}`$, where $`G`$, $`M_{BH}`$ and $`c`$ are the gravitational constant, mass of the black hole and the velocity of light respectively. We use Paczyński & Wiita (1980) potential $`\varphi (x)=\frac{1}{2(x1)}`$ to describe the flow around a Schwarzschild black hole. The basic dimensionless hydrodynamic equations which govern the infalling matter in the steady state are given by:

(a) The radial momentum equation:

$$\vartheta \frac{d\vartheta }{dx}+\frac{1}{\rho }\frac{dP}{dx}+\frac{\lambda _{Kep}^2\lambda ^2}{x^3}=0,$$
$`(1a)`$
(b) The continuity equation:

$$\frac{d}{dx}(\mathrm{\Sigma }x\vartheta )=0,$$
$`(1b)`$
(c) The entropy equation:

$$\frac{2na\rho \vartheta h(x)}{\gamma }\frac{da}{dx}\frac{a^2\vartheta h(x)}{\gamma }\frac{d\rho }{dx}=Q_{nuc}.$$
$`(1c)`$
We do not have an azimuthal momentum equation, since we assumed an inviscid flow where angular momentum is conserved. The sound speed is defined by, $`a^2=\frac{\gamma P}{\rho }`$, where, $`P`$ and $`\rho `$ are the pressure and the density respectively and $`\gamma `$ is the adiabatic index. Also, $`\lambda `$ is the conserved specific angular momentum (in units of $`2GM_{BH}/c`$) of the infalling matter, $`\lambda _{Kep}`$ is the specific angular momentum in a Keplerian disk, $`\lambda _{Kep}^2=\frac{x^3}{2(x1)^2}`$, $`\mathrm{\Sigma }`$ is the vertically integrated density, $`h(x)`$ is the half thickness of the disk \[$`ax^{1/2}(x1)`$\], $`n=\frac{1}{\gamma 1}`$ is the polytropic index, $`Q_{nuc}`$ is the height averaged heat generation/absorption due to the nuclear effect of the disk. In order to isolate the effect of nucleosynthesis alone, we neglect all the other heating/cooling terms. The factor $`Q_{nuc}`$ may be positive or negative depending on whether net nuclear energy is locally exothermic or endothermic. In future, we plan to include viscous heating and radiative cooling terms.

The energy release rate for $`pp`$ reaction is given by,

$$Q_{pp}=\frac{\rho _p^2}{2}<\sigma v>_{pp}q_{pp},$$
$`(2)`$
where, $`\rho _p=\rho N_{Av}X_p/A_p`$ is the number density of protons. The energy absorption rate due to $`D`$ and $`{}_{}{}^{4}He`$ dissociation is,

$$Q_{disso}=\rho _i<\sigma v>_iq_i,$$
$`(3)`$
where, the number density of the $`i`$th species (proton, deuterium or helium for this case) is $`\rho _i=\rho N_{Av}X_i/A_i`$. $`N_{Av}`$ is the Avogadro Number, $`X_i`$ is the mass abundance and $`A_i`$ is the mass number of the $`i`$th element and $`q_i`$ is the Q-value of the corresponding reaction. The reaction rates have been taken from standard literature (e. g., Clayton, 1983). For simplicity, and in order that this work may be generalized even when a large number of nuclear species are considered, we can express this as,

$$<\sigma v>_i=exp(g_i),$$
$`(4)`$
where,

$$g_i=c_i^1+c_i^2/t9+c_i^3/t9^{1/3}+c_i^4t9^{1/3}+c_i^5t9+c_i^6t9^{5/3}+c_i^7log(t9)$$

is a seven parameter family of functions (Thielemann, 1980). For different reactions, the constant coefficients $`c_i`$s will be different. In our case, $`Q_{nuc}`$ is the dimensionless height integrated nuclear energy generation rate $`Q_{nuc}=(Q_{pp}+Q_d+Q_{He})h(x)`$.

Using Eqs. (1b-1c), we eliminate all $`\frac{d\rho }{dx}`$ and $`\frac{da}{dx}`$ terms and get,

$$\frac{d\vartheta }{dx}=\frac{\frac{(\gamma 1)}{\vartheta }(Q_{pp}+Q_d+Q_{He})+(\gamma +1)[\frac{1}{2(x1)^2}\frac{\lambda ^2}{x^3}]\frac{a^2(5x3)}{x(x1)}}{\frac{2a^2}{\vartheta }(\gamma +1)\vartheta }.$$
$`(5)`$
At the sonic point, numerator and denominator both simultaneously should be zero, so,

$$\frac{2a^2}{\vartheta }(\gamma +1)\vartheta =0$$
$`(6)`$
and

$$\frac{(\gamma 1)}{\vartheta }(Q_{pp}+Q_d+Q_{He})+(\gamma +1)[\frac{1}{2(x1)^2}\frac{\lambda ^2}{x^3}]\frac{a^2(5x3)}{x(x1)}=f(a,\vartheta )=0.$$
$`(7)`$
From (6), we can get the Mach number at the sonic point as,

$$M_c=\sqrt{\frac{2}{\gamma +1}}.$$
$`(8)`$
From (7) and (8) we can obtain the sound speed at the sonic point.

The energy and entropy at the sonic points are expressed as,

$$E_c=\frac{1}{2}v_c^2+na_c^2\frac{1}{2(x_c1)}+\frac{\lambda _c^2}{2x_c^2},$$
$`(9)`$
$$\dot{}=a_c^{(2n+1)}x_{c}^{}{}_{}{}^{3/2}(x_c1)v_c.$$
$`(10)`$
Unlike the case of a flow with conserved energy (Chakrabarti, 1989; hereafter C89) one cannot compute everything about the flow using only two parameters, i.e., energy and specific angular momentum. In the present context of nuclear energy generation, one requires to supply the accretion rate as well which in turn determines the density of the flow. This is explicitly needed in eq. (7) above.

During the infall of matter, rate of change of mass fraction for proton, deuterium and helium can be written as,

$$\frac{dX_p}{dt}=\vartheta \frac{dX_p}{dx}=X_pexp(g_{pp})\rho \frac{X_dexp(g_d)}{2},$$
$`(11)`$
$$\frac{dX_d}{dt}=\vartheta \frac{dX_d}{dx}=X_dexp(g_d)\frac{X_{He}exp(g_{He})}{2}X_pexp(g_{pp})\rho ,$$
$`(12)`$
$$\frac{dX_{He}}{dt}=\vartheta \frac{dX_{He}}{dx}=X_{He}exp(g_{He}).$$
$`(13)`$
During the Runge-Kutta integration, these linear differential equations are also solved simultaneously to obtain their mass fractions at each point.

## 3 Calculation of Photon Number Density

Before we present results of our global analysis, we like to spend some time on how the photon number density is computed inside the transonic flow. Usual treatment of photo-dissociation in the literature assumes the presence of a black body radiation in the flow. However, advective flow close to a black hole is likely to be optically thin, and they would intercept soft photons from Keplerian disks located farther out. These soft photons are Comptonized and in reality, the photo-dissociation would be done by these Comptonized photons and the synchrotron photons rather than black body photons. We need to compute the number of hard photons above $`2`$ MeV in order to check if this is much larger compared to the deuterium nuclei.

If one considers a pure black body spectrum then the photon number density can be written as,

$$N_\nu ^{BB}=\frac{8\pi }{c^3}\left(\frac{kT}{h}\right)^3_{x_m}^{\mathrm{}}\frac{x^2dx}{e^x1}$$
$`(14)`$
where, $`x=\frac{h\nu }{kT}`$. But the spectrum of photons is not purely blackbody in nature. To calculate the number density of photon in our system we use results of Chakrabarti (1997, hereafter C97). We concentrate on Fig. 4 of C97 where variation of photon spectrum other than black body is discussed. In Fig. 4a of that paper, spectral variation is shown with energy for different values of Keplerian accretion rate (in units of Eddington rate) $`\dot{m}_d`$. For a complete dissociation of deuterium, at least $`2`$ MeV energy is needed (although with the energy of slightly less than $`2`$ MeV dissociation may start.). We first concentrate on those cases of $`\dot{m}_d`$ where curves are extended at least to the region greater than $`2`$ MeV. Each of these curves has one peak due to multi-colour black body and the other is a hump due to Comptonization. We concentrate on the hump at higher frequency to compute the number of hard photons. For the dissociation of a helium nuclei one requires a higher energy ($`32`$ MeV). In Fig. 4a of C97, except for $`\dot{m}_d=1.5`$ where there is a large peak due to the blackbody emission, other cases have significant amount of hard photon. Along with these photons, we include hard photons due to synchrotron radiation as obtained from a standard distribution $`I_\nu \nu ^{1.25}`$ normalized in a way such that its number is the same as that of the Comptonized photons at the ‘knee’ where blackbody distribution meets the Comptonized power-law component. This procedure provides a realistic estimate of hard photon. In any case, our main aim is to see if the resulting photon numbers are very large compared to the number of deuterium or helium nuclei. So the exact number is not important.

Following the treatment of Wagoner (1969), we completely turn off the dissociation process if the number of hard photons at a given radius is fewer than about a thousand times that of the deuterium. The same consideration applies for the dissociation of helium nuclei as well.

## 4 Results

Here we discuss a few results assuming matter falls towards a black hole of mass $`10M_{}`$. We choose the polytropic index $`n`$ to be $`3`$. We assume that the accretion flow consists of both Keplerian and sub-Keplerian matter (Chakrabarti & Titarchuk, 1995). We assume that these two rates are $`\dot{m}_d`$ and $`\dot{m}_h`$ respectively and $`\dot{m}_d+\dot{m}_h=2`$ for concreteness. We first compute the electron and proton temperatures in the flow by using two temperature hydrodynamic equations numerically (Chakrabarti & Titarchuk, 1995) for different $`\dot{m}_d`$ and $`\dot{m}_h`$ rates and use the resulting (cooler) temperatures to compute nucleosynthesis in our analytical code. Carrying out of nucleosynthesis work in presence of Comptonization would be otherwise very much prohibitive. We verify a posteriori from the solution with nucleosynthesis that such an approximation is valid. In all the three cases listed below we choose $`\lambda =1.6`$ and outer edge of the flow at $`x_o=1149`$. A flow without nucleosynthesis has a smooth solution passing through the sonic point at $`3.05`$ with these parameters provided it is launched at the outer edge at $`x_o`$ with a radial velocity $`v_o=8.8566\times 10^5`$ and sound speed $`a_o=0.061966`$.

### 4.1 Case I

Here we choose $`\dot{m}_d=0.05`$, $`\dot{m}_h=1.95`$. In this case, the photon number density is $`1.61\times 10^{17}`$/cc. We first check the global behaviour of the infalling matter, such as the position and nature of the sonic points if nucleosynthesis is turned on. Following C89, we plot in Fig. 1a, energy $``$ at the sonic point as a function of the entropy accretion rate $`\dot{}`$ at these points for a fixed specific angular momentum ($`\lambda =1.6`$) of the sub-Keplerian and Keplerian mixture. Dashed curve indicates the behaviour when nucleosynthesis is not included. This behavior is similar to that observed in C89. The solid curve indicates the variation of energy with entropy accretion rate at the sonic point when nucleosynthesis in the flow is included. The dashed curve is shifted vertically by $`0.01`$ in order to see the effects of nucleosynthesis prominently. If the ratio of the number density of photon and the nuclear species to be dissociated in the disk is less than $`10^3`$ the dissociation is turned off because of the deficiency of photons (Wagoner, 1969). For the present set of parameters, it is seen that till $`\dot{}10^5`$ nucleosynthesis does not affect the sonic points. In the range $`\dot{}10^510^6`$, the outer sonic points disappear as the solid curve shows a distinct gap in that range where nucleosynthesis is prominent. Subsequently, a new branch is originated. As the inner sonic point is not affected, we can conclude that in these parameters the stability of the disk is not affected. Matter can pass through either the inner sonic point (I) and or through the outer sonic point (O) if the location of the latter point is far enough. Flows with sonic points located at the intermediate distances are unstable.

In Fig. 1b, we show the variation of the Mach number as a function of radial distance. The inner sonic point is chosen at $`x_{\mathrm{in}}=3.05`$. The velocity of the matter at the outer edge of the disk (assumed at $`1149`$) is chosen to be $`v_o=8.8566\times 10^5`$ as this allows a smooth solution through the inner sonic point when nuclear energy is not included. The long-dashed curve is the resulting solution. The sound speed at the launching point is $`a_o=0.061966`$. All quantities are expressed in geometric units. When nucleosynthesis is included, in order that the flow passes through the same inner sonic point smoothly, we had to raise the sound speed to be $`a_o=0.0791819`$. The solid curve indicates this solution. If the flow with nucleosynthesis is given the same $`a_o=0.061966`$ as that of the smooth solution without nucleosynthesis (long dashed curve), the resulting solution becomes totally sub-sonic and is unphysical (dash-dotted curve) in the context of black hole accretion. For comparison, we present two other solutions (short dashed and dotted) which begin with different sound speed, and are also found to deviate from the true transonic solution (solid curve). Fig. 2 shows the variation of the energy release due to nucleosynthesis as a function of the logarithmic radial distance of the flow (solid curve). The energy release is basically due to $`pp`$ reaction as deuterium burning took place instantaneously after launching the flow. Note that the release is at the most $`10^{17}`$ ergs/gm/sec which the specific energy of the flow with an inner sonic point could be around $`0.001c^2`$ergs/gm. Thus the nuclear cooling time scale is $`10`$s which is large compared to, infall or dynamical time scales $`t_d10^3`$s. Thus, energetically nuclear energy release is not significant.

### 4.2 Case II

Here, we choose a comparatively hotter flow. We use $`\dot{m}_d=0.001`$, $`\dot{m}_h=1.999`$ photon number density is calculated as $`3.02\times 10^{21}`$/cc. The global behaviour and particular solutions are shown in Fig. 3(a-b). They have similar properties as those of Case I. In Fig. 3a, we clearly see the disappearance of the outer sonic points between $`O_1`$ and $`O_2`$. The deviation from a ‘non-nuclear’ solution is higher, and one requires to launch the flow with a higher speed of sound: $`a_o=0.082167`$. In this case, the gas has to be heated at a much higher temperature in order that it may pass through the same sonic point. In that sense the flow is prone to become more unstable than that of Case I. Fig. 2 (dashed curve) shows energy release for this case. A hotter gas releases somewhat higher energy.

### 4.3 Case III

In this case, we choose a cooler flow so that deuterium is dissociated closer to the black hole. We use $`\dot{m}_d=0.4`$, $`\dot{m}_h=1.6`$ and specific angular momentum of the mixture of Keplerian and sub-Keplerian gas is $`\lambda =1.6`$ as before. Since accretion rate is large, there are virtually no Comptonized hard photons. However, there are synchrotron photons which dissociate deuterium. These synchrotron photons are computed from a power-law intensity distribution with a slope of $`1.25`$ joining with blackbody spectrum at roughly the same place where Comptonized power-law joins. Fig. 2 shows the total energy release (dotted upper curve) as well as energy absorption due to deuterium in this case (dotted, marked by $`Q_D`$). As matter comes closer to the black hole, some deuterium is produced by the $`pp`$ reaction which is also instantaneously destroyed. This raises energy absorption closer to a black hole. Note also that only absolute value of $`Q_D`$ is drawn here.

Fig. 4 gives the variation of the abundance of deuterium, proton and neutron with logarithmic radial distance. In this case, deuterium is totally destroyed at around $`100`$ Schwarzschild radii causing a rise in the protons and neutrons.

## 5 Concluding remarks

We studied analytically the effects of nucleosynthesis on black hole accretion flows of constant angular momentum. We took only the basic reactions such as $`pp`$ reaction producing deuterium and destruction of deuterium and helium due to photo-dissociation since (a) it is easier to handle only a few reactions analytically and (b) these reactions primarily control the energetics of the flow. Since we were dealing with an optically thin flow, the radiation could not be treated as a thermal blackbody and the photon number density had to be calculated from the actual thermal Comptonized spectrum and from synchrotron emission from coronal electrons. Photons with energy higher than the threshold value were used to destroy deuterium and helium.

Though, the energy released or absorbed were found to be not very significant compared to the rest mass of the flow it has a strong influence on the stability properties at the transonic region. For instance, we find that when nucleosynthesis in included, the sound speed of the injected matter must be higher by at least $`25`$ to $`30`$ percent in order that the flow passes through inner sonic point. In other words, the temperature of the flow at the injection point must be $`50`$ to $`60`$ percent higher. If a physical process (such as magnetic heating) could not be found to enable this, the flow will not be transonic and would be unstable (See, Figs. 1b and 3b).

So far, in the literature, study of accretion flows around black holes had been done using two-step numerical procedure — pure hydrodynamics is followed by pure nucleosynthesis and results of nucleosynthesis is fed back into hydrodynamics. This procedure was repeated for the whole disk till a reasonable convergence occurs. In the present paper, we combined hydrodynamics and nucleosynthesis for the first time and studied analytically how nucleosynthesis influences topology of the flow. In some region of the parameter space, the outer sonic point is found to be absent. How would a flow behave in this region? It is conjectured that matter would try to pass through ‘some’ outer sonic point, becomes supersonic and subsequently pass through the inner sonic point. Since the outer sonic point is not present, the searching procedure to the true sonic point should give rise to an unstable or oscillating solution. Only a time dependent numerical simulation can tell if this conjecture is correct. This will be attempted in the future.


text:
# Coulomb Explosion and Thermal Spikes

## Abstract

A fast ion penetrating a solid creates a track of excitations. This can produce displacements seen as an etched track, a process initially used to detect energetic particles but now used to alter materials. From the seminal papers by Fleischer et al.  to the present , ‘Coulomb explosion’ and thermal spike models are treated as conflicting models for describing ion track effects. Here molecular dynamics simulations of electronic-sputtering, a surface manifestation of ion track formation, show that ‘Coulomb explosion’ produces a ‘heat’ spike so that these are early and late aspects of the same process. Therefore, differences in scaling are due to the use of incomplete spike models.

Since Thompson and Rutherford there has been interest in the track of excitations produced by an energetic ion penetrating a solid. The fast ion excites the electron cloud producing excitons and electron-hole pairs in an insulator . The electrons expelled from the track create additional excitations, then cool to the lattice, screening and finally recombining with the holes. They eventually return to the ground state by radiative or non-radiative processes. Prior to this decay an excited region with a net repulsive energy persists. Fleischer et al.  (hereafter FPW) proposed that this repulsion could produce displacements altering the material along the particle’s track. When accounting only for the holes the process has been called ‘Coulomb explosion’, suggested to produce amorphized tracks , cratering, and sputtering . Tracks have also been observed in liquids , semiconductors  and even metals  when the excitation density in the track is sufficiently large.

In describing thresholds for track registration in insulators, FPW compared their ‘Coulomb explosion’ model to the often used thermal spike model for defect production . Such a comparison is not only of historical interest but persists in the recent literature . However, one model describes a mechanism for energy-input (electrostatic potential energy into kinetic energy), whereas the other describes the transport of energy out of a ‘heated’ cylindrical region. Here we present a molecular dynamics (MD) description of the response of a solid to a repulsive track. We show that at high excitation densities the primary effect is the production of a cylindrical ‘heat’ spike. Therefore, distinctions made between ‘Coulomb explosion’ and spike models can be artifacts of the approximations used.

MD simulations are carried out for a model solid with a surface because of our interest in electronically-induced sputtering, the surface manifestation of track formation in insulators . As in track registration, both spike models  and ‘Coulomb explosion’ models  have been employed to parametrize laboratory data for electronic sputtering. In spike models the atoms in the ion track are assumed to have a radial temperature profile determined by the energy deposited per unit path length, $`(dE/dx)`$. When the energy of an atom exceeds a barrier, defect production in the bulk or evaporation at the surface can occur. From the time dependent surface temperature, $`T_{surf}`$, the spike sputtering yield, $`Y_S`$, the number of atoms ejected per ion track is calculated integrating the surface flux over time and area. When a radial diffusion equation is used to determine $`T_{surf}`$, then $`Y_S\left(dE/dx\right)_{eff}^2`$ for a fixed track radius at high $`(dE/dx)_{eff}`$ , where $`(dE/dx)_{eff}`$ is the fraction of $`dE/dx`$ going into energetic non-radiative processes . This has been used to describe sputtering produced by a track ‘heated’ by individual repulsive decay events , secondary electrons  and a repulsive track . However, when conditions leading to the $`\left(dE/dx\right)_{eff}^2`$ dependence should apply, the model fails because the energy transport is not diffusive . A melt and a pressure pulse control the energy transport leading to ,

$$Y_S0.18(r_{cyl}/U)(dE/dx)_{eff}$$
(1)
where $`U`$ is the solid’s cohesive energy and the effect depth scales with $`r_{cyl}`$, the initial mean radius of the spike. Here we show that a repulsively-induced ‘heat’ spike is produced in an ion track over a broad range of excitation densities and the transport is not diffusive. Therefore, a correct spike model must be used.

We simulate the evolution of a track of repulsive energy that might be produced by a fast incident ion. Since sputtering  and crater formation in both amorphous and crystalline, atomic or molecular solids can be scaled over a broad range of material properties, we use a model solid made of atoms interacting via Lennard-Jones (LJ) potentials. We use parameters for condensed gas solids, but the results scale with the energy and length parameters. Also, we showed earlier that scaling was maintained using more complex potentials . Interactions occur between all atoms within a cut-off radius $`r_{cut}2.54l`$, where $`l`$ is the mean atom spacing. The (001) layer spacing in the fcc lattice is $`l_s0.8l`$. MD simulations focusing on damage have been made for repulsive energy between a few charges in LiF , a distribution of holes in Si , and using an assumed velocity distribution . Here we follow the conversion of a track of repulsive energy into atomic motion in a solid.

On electronically exciting a condensed gas solid, the track of excitations can produce large sputtering yields at high $`dE/dx`$, a process relevant to icy outer solar system bodies. Ejection is presumed to be due to the net repulsive energy, but a quantitative description is lacking. At high ionization densities the net repulsion can persist if the holes have low mobility and are not sufficiently screened during the time displacements could occur. Even when holes are neutralized, the excited atoms have overlapping charged clouds which act repulsively at high $`dE/dx`$ . Here we describe the net repulsion between ‘excited’ atoms in a track using $`V=\left(e^2/r\right)\mathrm{exp}\left(r/a\right)`$, where $`a`$ is an average screening constant . We do not distinguish between closely spaced, partially screened holes or overlapping excited species, since in both cases the electrons screen the interactions between ‘excited’ neighbors . The interactions of ‘excited’ atoms with unexcited neighbors was left unchanged; including a weak polarization had a small effect. For energy conservation, a larger cut-off, $`r_{cut}^{coul}=7a`$, is used for interactions between ‘excited’ species.

The material is ‘excited’ by instantaneously changing the potentials between atoms in the track. The resulting velocities and positions of all particles are then followed . This is done for a number of excitation densities and screening constants. In a second set of calculations the excitations were quenched (recombination). This was done statistically so the average number of excitations in the track decayed exponentially, $`\mathrm{exp}\left[t/\tau \right]`$ where $`\tau `$ is the quenching (neutralization) time. The number of excitations per unit path length, $`dJ/dx`$, is related to $`dE/dx`$ for fast ions  and the mean radius of the distributions, $`r_{coul}`$, depends on the speed of the incident ion. For given track parameters \[$`dJ/dx`$, $`r_{coul}`$ (here $`l_s`$)\] excited species are chosen randomly. On excitation, atoms are sputtered if they cross a plane $`2r_{cut}^{coul}`$ above the ‘top’, here the (001) surface. The number of atoms ejected in each run is called the yield, which is related to the track of damage in FPW but is more easily measured. Sample depth was chosen to be at least twice $`r_{cut}^{coul}`$. Different boundary conditions did not change the average yields significantly. The sample size $`(3x10^43x10^5)`$ and simulation times ($`1580`$ ps) were adjusted to $`a`$ and $`dJ/dx`$. Extending times by tens of picoseconds or doubling the thickness did not change the average yield. The distribution in the yield is broad, especially for small $`l_s(dJ/dx)`$ where sputtering occurs when two excitations are produced close together . Results, averaged over $`10200`$ excitation distributions, are given in Fig. 1a as a function of charge density, $`l_s(dJ/dx)`$, for 3 values of the screening constant $`a`$. The yield is seen to be quadratic in $`dJ/dx`$ for each $`a`$ and increases non-linearly with $`a`$; lines $`Y=14.\mathrm{ln}\left[1.3\left(a/r_{coul}\right)\right]\left[l_s\left(dJ/dx\right)\right]^2`$ give a good fit . The yields as a function of neutralization time, $`\tau ,`$ for a fixed $`(dJ/dx)`$ and $`a`$ are given in Fig. 1b scaled by the fit in Fig 1a. Because neutralization is treated stochastically, the size of the yield is affected even for relatively large $`\tau `$ and only becomes independent of the neutralization rate for $`\tau 10\tau _D`$, with $`\tau _D`$ the Debye period for the lattice ($`0.5`$ ps) . However, the yield is still seen to be nearly quadratic in $`(dJ/dx)`$ for the full range of $`\tau `$ with an additional weak dependence on excitation density.

In MD simulations details of the energy transfer and sputtering can be extracted. At very low excitation densities $`(Y1)`$ the probability of neighbors being excited near the surface leads to ejection. At the highest excitation densities $`(Y5)`$, the repulsive energy is transferred to neighbors in times $`0.2\tau _D`$ producing a cylindrically heated region. Comparing the evolution of the resulting radial temperature profiles with our earlier simulations of ‘heat’ spikes, the energy transport is found to be similar. That is, we find for the repulsively produced heat spike transport is not diffusive , a melt and a pressure pulse are formed and these control the energy transport. Seiberling et al.  suggested that a repulsive track could produce a heat spike but used the standard diffusive model to obtaining a yield dependence very different from that calculated here.

From the energy distributions for the ejecta at high excitation density two regions can be seen in Fig. 2. At large ejecta energies, $`E`$ $`10U`$, the spectra has peaks due to prompt ejection from the upper surface layers of the initial track. This accounts for $`20`$% of the ejecta at large $`dJ/dx`$ but dominates at very small $`dJ/dx`$. On the other hand, the principal component of the ejecta in Fig. 2 exhibits an energy distribution like that found in our studies of ejection from a narrow cylindrical heat spike shown as the dashed line . That is, there is a broad, quasi-thermal distribution at low ejecta energies, $`E<U`$ which gives way to a non-Maxwellian, $`E^2`$ dependence at $`E>U`$. The latter dependence is characteristic of low energy cascades in a solid  but differs from thermal spike model predictions. Surprisingly, although the ejection process changes in going from the lowest to the highest $`dJ/dx`$, there is no dramatic change in the dependence of the yield on $`dJ/dx`$ in Fig. 1. This is due to the nearly linear dependence of spike yield on $`(dE/dx)_{eff}`$ discussed below.

From the above, the primary effect of the repulsive decay at the high excitation densities is the production of a ‘heat’ spike. This is also shown quantitatively using the MD calculation to determine the repulsive energy driving atomic motion in times short compared to desorption times. At $`0.2\tau _D`$ the energy density deposited via repulsion is $`(dE/dx)_{rep}0.15e^2(dJ/dx)^2`$ for $`a=1.13l_s`$and $`r_{coul}l_s`$. This energy is localized in a cylindrical region of radius $`r_{rep}`$ $`2.6l_s`$ for our values of $`r_{coul}`$. Using $`(dE/dx)_{eff}=(dE/dx)_{rep}`$ and $`r_{cyl}=r_{rep}`$ as initial conditions in the expression for the yield from a cylindrical ‘heat’ spike in Eq. 1 gives a yield that has the same dependence on $`dJ/dx`$ and a size within 10% of those in Fig.1a. Therefore, if a repulsive region is sustained, a ‘heat spike’ is formed at high $`dJ/dx`$ which determines the subsequent energy transport, sputtering and displacements. Therefore, ‘Coulomb explosion’ and spikes are the early and late aspects of the same process, so that contradictions for effects occurring later than $`0.2\tau _D`$ are due to incorrect descriptions of the ‘heat’ spike. Below we examine under what conditions a ‘Coulomb’ track is sustained long enough to produce a ‘heat’ spike.

Since neutralization quenches the energy of a spike, we simulated the neutralization and screening of holes by free electrons. That is, we make a separate simulation in which we track the cooling of the ‘excited’ electrons in the field of the track of positive charges. The lattice is like that in the MD simulations and we choose appropriate mean-free paths for elastic, $`\lambda _s`$, or inelastic, $`\lambda _{in}`$, scattering, both assumed to be isotropic . A time dependent $`a`$ or an average $`a`$ and $`\tau `$ can be determined from such a simulation. Instead, we give in Fig. 3 the instantaneous total potential energy density in the track vs. time. This is the energy available for repulsive heating of the lattice. The interactions between the positive charges in the track and the electrons are Coulombic beyond an atomic radius, $`a_o`$. Inside $`a_o`$ the potentials gradually become flat giving a binding energy of $`8.6`$ eV. We examined a number of potential forms inside $`a_o`$ and initial electron energies and the trends are the same. Electrons are ‘excited’ by receiving, on the average, 15 eV, roughly the energy dissipated to the lattice (W-value minus ionization energy). This does not treat the very fast electrons which would further slow neutralization.

For only elastic scattering, the total energy in the track in Fig. 3 is sustained with a very large effective $`a`$. Examining the radial dependence of the electron density, a fraction of the electrons cool rapidly by electron-electron collisions partially neutralizing the track. This is dielectronic recombination, also called Auger recombination. The remainder of the electron cloud, now ‘hotter’ but more fully screened, expands to a larger average radius. Varying $`\lambda _s`$ affects the electron density distribution, whereas including an inelastic energy loss, $`\mathrm{\Delta }\epsilon _{in}`$, decreases the potential energy in the track. This also heats the lattice, but here we are interested in the repulsive energy. For a $`\mathrm{\Delta }\epsilon _{in}`$ and $`\lambda _{in}`$ roughly corresponding to electrons in a polymer or in ice (dotted curve) , it is seen in Fig. 3 that most of the initial track energy remains during the time it takes to produce a ‘heat spike’ repulsively, $`0.2\tau _D`$. This would correspond to an intermediate $`\tau `$ and relatively large $`a`$ in Figs. 1b and 2a. Therefore, for a reasonable inelastic energy loss, screening is not sufficient to quench the repulsive production of a ‘heat spike’. Increasing $`\mathrm{\Delta }\epsilon _{in}`$ (dot-dashed curve) or decreasing $`\lambda _{in}`$ leads to more rapid but still incomplete quenching over the relevant time. Quantities for specific measurements need to be used and the role of excited states in the track must be included.

In this paper we close the circle on one aspect of an old  but topical  problem, the effect on a solid of a track of excitations produced by a fast penetrating ion. A MD simulation of repulsive explosion in a track of excitations has been carried out and the effect of neutralization has been examined. This was done for narrow screened Coulomb spikes in a condensed-gas solid to constrain the simulation size. However, we showed earlier that the yields scale with $`U`$ and $`(dE/dx)_{eff}`$. Whereas at low excitation densities ejection can occur if neighboring ‘excited’ species are formed near the surface , at the higher excitation densities the repulsive energy in the track produces a ‘heat’ spike. Therefore, ‘Coulomb explosion’ and spike models are the early and late aspects of the repulsive decay process. The spike formed by the repulsive heating of the lattice can lead to displacements and ejection (electronic sputtering), examined here. Because the yield is roughly linear in $`(dE/dx)_{eff}`$ at high excitation densities and $`(dE/dx)_{eff}`$ is proportional to $`(dJ/dx)^2`$, the yield from a repulsive track is seen (Fig.1a) to vary smoothly, having the same quadratic dependence on $`dJ/dx`$ in going from low to high excitation densities although the character of the ejection process changes. At high excitation densities the yield can now be estimated by substituting the repulsive energy deposition into the new, useful expression for spike sputtering in Eq. 1. At higher excitation densities than those studied here, the more distributed ‘heat’ spike produced by the electron cooling to the lattice might also be important . Earlier applications of a ‘heat’ spike for describing track formation or electronic sputtering foundered on the use of a spike model in which the energy transport was described incorrectly. That is, the transport is not diffusive  and the local kinetic energy distribution is not Maxwellian (Fig. 2). Calculations should now focus on accurately describing the interactions in an ion track for times $`0.2\tau _D`$ including electrons, holes and excitons, as the latter also contribute to the net repulsive energy.

###### Acknowledgements.

We thank M. Liu for the neutralization calculations, H. Urbassek, R. Vidal, and R. Baragiola for comments, the NSF Astronomy and Chemistry Divisions for support.