text:
# 𝛾-Ray Lines From Asymmetric Supernovae

## 1 Introduction

Among the many surprises that supernova (SN) 1987A brought astronomers was the early emission of X-rays and $`\gamma `$-rays (X-rays: e.g. Dotani et al. 1987; Sunyaev et al. 1987; $`\gamma `$-rays: e.g. Cook et al. 1988; Mahoney et al. 1988; Matz et al. 1988). This high energy emission, arising from the decay of 56Co, appeared nearly 6 months earlier than was predicted by theoretical models (e.g. Pinto & Woosley 1988a, Arnett et al. 1989) and led theorists to conclude that the 56Ni, produced deep in the core of this exploding star, had mixed into the outer layers of the supernova ejecta (e.g. Pinto & Woosley 1988b, Arnett et al. 1989 and references therein). Additional evidence all seems to support the mixing of SN 1987A’s ejecta: models of the supernova light curves (Arnett 1988; Shigeyama, Nomoto, & Hashimoto 1988; Woosley 1988; Shigeyama & Nomoto 1990) and explanations of spectral line widths (Haas et al. 1990; Spyromilio, Meikle, & Allen 1990; Tueller et al. 1990). SN 1987A is not peculiar in this mixing. Many supernovae show evidence of mixing in their spectra (e.g. Spyromilio 1994, Fassia et al. 1998) and the light curves and spectra of Type Ib SN seem to be best fit by mixed models (Shigeyama et al. 1990; Woosley & Eastman 1997). It appears that mixing is a generic process in core-collapse supernovae.

These results have stimulated a series of multi-dimensional hydrodynamical simulations trying to produce the observed mixing (Arnett, Fryxell, & Müller 1989; Hachisu et al. 1990, Fryxell, Arnett & Müller 1991; Müller, Fryxell & Arnett 1991; Herant & Benz 1991; Herant & Benz 1992; Herant & Woosley 1994; Kifonidis et al. 2000; Kifonidis et al. 2003). Although these simulations seem to be able to explain the mixing in Type Ib supernovae (Kifonidis et al. 2000), none of these results are able to explain the extensive mixing observed in SN 1987A. Possible solutions to this shortcoming have been proposed: (1) perhaps the decay of 56Ni injects enough energy to force additional mixing (Herant & Benz 1992) or (2) convection in the pre-collapse core provides enough seeds to enhance mixing (Herant & Benz 1992). A third possibility is that the supernova explosion itself is asymmetric (Nagataki et al. 1997; Nagataki et al. 1998; Nagataki 2000 and references therein). Nagataki et al. (1998) found that not only could slight asymmetries in the supernova explosion produce the required mixing to explain 1987A, but they could also explain anomalies in the nucleosynthetic yields produced by several supernovae.

The evidence for such global asymmetries in supernova explosions has been growing steadily. One piece of evidence arises from attempts to understand the high space velocities of neutron stars. The high observed velocities of pulsars, along with evidence of neutron star/remnant associations, and the formation scenarios of neutron star binaries all suggest that neutron stars are given strong kicks at birth. These kicks are most easily explained by some asymmetry in the supernova explosion (see Fryer, Burrows, & Benz 1996 for a review). In addition, the most straightforward explanation of the large polarization seen in core-collapse supernovae (see Wang et al. 2001; Leonard & Filippenko 2001 and references therein) is that the explosion driving these supernovae is inherently asymmetric (Höflich 1991).

The asymmetries are believed to have their origin in the explosion mechanism itself. For instance, even if the collapsing star is initially spherically symmetric, some asymmetry can be produced due to convection taking place in, and above, the proto-neutron star (Herant et al. 1994; Burrows, Hayes, & Fryxell 1995; Janka & Müller 1996). To date, these asymmetries are not extreme enough in the theoretical models to explain the mixing. Large asymmetries may occur if the collapsing star is asymmetric due to nuclear burning (Burrows & Hayes 1996; Lai 2000), however, sufficiently large departures from spherical symmetry have only been produced by assuming extremely asymmetric collapsing cores (Burrows & Hayes 1996). Alternatively, rotation can produce strong asymmetries in the supernova explosion (Mönchmeyer & Müller 1989; Janka & Mönchmeyer 1989; Fryer & Heger 2000, Khokhlov et al. 1999). The nature of these asymmetries depends upon the angular momentum profile of the collapsing star and, although most calculations predict jet-like explosions along the rotation axis, some calculations imply that an equatorial explosion could occur (Mönchmeyer & Müller 1989).

In this paper, we follow supernova explosions in 3-dimensions from 100 seconds to 1 year after core-bounce. We model a series of initial explosion conditions with both jet-like “axial” and equatorial asymmetries of varying degree. In §2, we describe these simulations and their results with comparisons to past work. In §3, we discuss the $`\gamma `$-ray emission from these explosions and present calculations of $`\gamma `$-ray spectra as a function of time and viewing angle. We conclude with a disscussion of the observational prospects of these results, highlighting upcoming gamma-ray missions.

## 2 Explosion Simulations

For our hydrodynamic simulations, we have used the 15 Mprogenitor (s15s7b) by Weaver & Woosley (1993). This star has been evolved with a piston-driven explosion to 100 s after bounce, producing 0.24 Mof 56Ni. The total energy of this model is roughly $`1.5\times 10^{51}`$ erg with roughly $`1.0\times 10^{51}`$ erg in kinetic energy. As this explosion moves through the star, the shock hits composition boundaries where strong entropy gradients exist. When the shock hits these boundaries, Rayleigh-Taylor instabilities develop, which can grow and cause the star to mix (Chevalier & Klein 1978; Weaver & Woosley 1980). Our simulations model this mixing and concentrate on the effects that asymmetries have on it.

### 2.1 Numerical Schemes

We model a series of explosions from 100 s to 1 year after the launch of the supernova shock (Table 1) using a 3-dimensional smooth particle hydrodynamics code (see Warren et al. 2002 for details) based on the parallel oct-tree algorithm developed by Warren & Salmon (1993). This lagrangian code tracks the composition exactly. For 2 of our simulations, we included the energy injection from 56Ni and 56Co decay (and hence also trace the abundances of 56Co and 56Fe). For these 2 models, we assume that all of the decay energy is injected into the gas. At late times, this will overestimate the total energy deposited, as a fraction of this energy will escape, but since we would like an upper limit on the effects of 56Ni decay, and because the injection of energy is less important to the explosion dynamics at late times, this assumption is adequate. For decay energy, we assume that the total energy from decay is $`9.3\times 10^{16}`$ erg g-1 with 33% being released from 56Ni decay with a 6.1 d half-life and the other 67% being released from 56Co decay with a 77 d half-life (Colgate, Petschek, & Kriese 1980). This is similar to what Herant & Benz (1992) and Herant & Woosley (1994) used for their models. For an equation of state, we use the same “low-density” equation of state developed by Blinnikov, Dunina-Barkovskaya, & Nadyozhin (1996) that is used in our progenitor model.

We mapped the Weaver & Woosley (1993) model with its 100 s long, spherically symmetric explosion into our 3D SPH code with 2.2 million variably-massed particles. We model the entire sphere, so there are no axis boundaries in this simulation and the neutron star mass at the center is added through an external gravitational force. Asymmetries are added artificially to the velocities only and, for most of the explosions, are added in such a way to insure that the total kinetic energy of the explosion is conserved. The two asymmetries we model are jet explosions:

$$v_{\mathrm{radial}}=(\alpha +\beta \times |z|/r)v_{\mathrm{radial}}^{\mathrm{sym}}$$
(1)
and equatorial explosions:

$$v_{\mathrm{radial}}=(\alpha \beta \times |x|/r)v_{\mathrm{radial}}^{\mathrm{sym}}$$
(2)
where $`v_{\mathrm{radial}}^{\mathrm{sym}}`$ is the velocity from the one dimensional calculation, $`v_{\mathrm{radial}}`$ is the radial velocity for the asymmetric setup, and $`x`$, $`z`$, $`r`$ are the x-position, z-position, and radius of the particle. The values of $`\alpha `$ and $`\beta `$ for each model are given in Table 1, along with the initial kinetic energy of each model. The magnitude of the asymmetries are guided by the results of Fryer & Heger (2000), who found that their rotating core-collapse simulations in 2-dimensions had velocities that were a factor of 2 higher in the pole than in the equator 1.5 s after bounce. We do not use the extreme asymmetric explosions of Khokhlov et al. (1999) which concentrate the explosion energy into a narrow jet. As we shall show (and as Nagataki 2000 argued), such extreme asymmetries are not required to explain outward mixing of nickel in supernovae. Although at the end of their simulation, the asymmetry in polar vs. equatorial velocities is growing, the trend as the shock moves outward will be to spherize the shock. So at this point, it is not clear whether the shock at 100 s will be more or less asymmetric than what Fryer & Heger (2000) found at the end of their simulation.

After mapping these models into our 3-dimensional SPH code, we then run the explosion out to 1 year. When the shock reaches the edge of the star, we assume it is moving through a vacuum. In reality, such a star will have a stellar wind atmosphere surrounding it, but the density of a 15 Mwind is so low, that for the purposes of our simulation, zero density material is appropriate. In addition, photon transport is not modeled in our calculations. However, prior to shock breakout, the photons are essentially trapped in the shock. By the time the shock breaks out, most of the internal energy has already been converted to kinetic energy, so although including photon transport will change the ionization state of the exploding stellar material, it does not affect the kinematics significantly.

The smooth particle hydrodynamic technique automatically adds a level of perturbation in the code. With our initial conditions, these perturbations are random with a maximum $`1\sigma `$ deviation in a shell of 5-7%. Such large deviations only occur at composition boundaries, specifically the silicon and oxygen shells, and match well the deviations arising from explosive oxygen and silicon flashes prior to collapse (e.g. Bazan & Arnett 1998).

### 2.2 Explosion and Nickel Distribution

Although we use the same 15 Mprogenitor that was used by Herant & Woosley (1994), it has a higher explosion energy and we map this model onto our 3-dimensional grid 200 s earlier than Herant & Woosley (1994) mapped their spherically-symmetric explosion simulation onto a 2-dimensional grid. Therefore, although their study has the closest similarities with our work, it is difficult to make direct comparisons to their simulations. Nevertheless, it is interesting to compare the velocity distribution of each chemical element from our 3-dimensional simulations with the 2-dimensional simulations of Herant & Woosley (1994) at 90 days (compare the lower right panel in Fig. 1 of Herant & Woosley 1994 with Fig. 1 in this paper). In Fig.1, the material labelled “hydrogen” includes all material in the hydrogen envelope (as did Herant & Woosley 1994). Similarly, by “nickel” we refer to both the distribution of nickel as well as its decay products (most notably $`{}_{}{}^{56}Co`$). Although the distribution of elements is similar in both the 2 and 3 dimensional simulations, the stronger 3-dimensional explosion causes all of the ejecta to be moving slightly faster than that of the 2-dimensional simulation and it is difficult to compare mixing instabilities.

Comparing the convective instabilities themselves is also difficult. In 3-dimensions, the “mushroom”-like structures formed by Rayleigh-Taylor instabilities are not so well defined, and don’t lie along any one plane. However, 4.3 hours into the explosion, it is clear that instabilities have developed (Fig. 2) and these instabilities ultimately mix nickel knots far out into the star (Fig. 3). At the start of the explosion, 56Ni is found only in the inner 1.6 Mof the star (1.3 Mbecomes the neutron star, so the nickel is limited to the inner 0.3 Mof ejecta). By the end of the simulation, 56Ni has mixed out nearly to 5 M, beyond the $``$4.5 Mboundary that marked the edge of the helium layer (Fig. 4). Unfortunately, this mixing is, if anything, less than the amount of mixing found in the 2-dimensional simulations of Herant & Woosley (1994). The fact that the mixing is less in 3-dimensions vs. 2-dimensions could be due to the lower effective resolution (we only have 2.2 million particles in 3-dimensions vs. 25,000 particles in the 2-dimensional simulations)111Note that the 2-dimensional simulations of Herant & Woosley (1994) have poor mass resolution and low resolution could be a problem in the 2-dimensional simulations as well as the 3-dimensional simulations. However, bear in mind that the turbulent inverse cascade behaves differently in 2 and 3 dimensions (2-dimensional inverse cascades drive energy to large scales whereas 3-dimensional simulations suggest the energy is driven to small scales and dissipated). It is likely that these differences cause the 2-dimensional simulations to produce more extended instabilities. In any event, it appears that neither our spherical 3-dimensional simulations nor the 2-dimensional simulations seem to give enough mixing to explain the observations of supernovae like 1987A.

Nagataki et al. (1998) and Nagataki (2000) found that they required mild asymmetries ($`v_{\mathrm{Pole}}/v_{\mathrm{Equator}}=2`$) to explain SN 1987A. Since we model a 15 Mstar, not a SN 1987A progenitor, it is difficult to both compare with this past work as well as constrain our results with observations of SN 1987A. But we can discuss the basic trends caused by asymmetries. Fig. 5 shows model Jet2 1 year after explosion. Note that although the density distribution has spherized as the shock propogates through the shallow density gradients of the red supergiant envelope (Chevalier & Soker 1989), the distribution of 56Co (the decay product of 56Ni) retains a large asymmetry. We discuss the effects of these asymmetries on the $`\gamma `$ ray emission in §3.

Like Nagataki et al. (1998) and Nagataki (2000), we find that the asymmetries broaden the velocity profile of 56Ni (Fig. 6). However, in our simulations, mild asymmetries ($`v_{\mathrm{Pole}}/v_{\mathrm{Equator}}=2`$) led to only a small increase in th maximum nickel velocity from 2500 km s-1 to 2900 km s-1. For Nagataki et al. (1998), such mild asymmetries increase the maximum nickel velocity from 2200 km s-1 to 3200 km s-1! This difference could be progenitor dependent, an effect of 3-dimensional vs. 2-dimensional convection, or the lack of resolution in our 3-dimensional models. Extracting the true cause of this difference awaits future calculations with similar initial conditions.

However, increasing the amount of asymmetry by another factor of 2 (Models Jet4, Eq4) causes some nickel to be ejected at very high velocities. The amount of mixing in these cases reaches extremes with the nickel well into the hydrogen envelope. Such mixing has decided signatures in both the emergence of the $`\gamma `$-ray line flux and the shape of the $`\gamma `$-ray lines (§3). Note that the energy released from the decay of nickel also helps to mix out the nickel (on par with the effects of mild asymmetries). Clearly, the energy released from the decay of nickel cannot be neglected in any accurate mixing calculation.

This mixing also has important repercussions for nucleosynthetic yields and the mass-cut for the remnant mass. Most black holes are formed in stars which produce supernova explosions that are too weak to throw off all of the stellar envelope and the subsequent fallback produces a black hole (Fryer & Kalogera 2001). In our models (Jet2, Sym+Decay, etc.), more than 10% of the nickel produced is ejected well beyond the helium core (Fig. 7). If this trend holds for more massive stars such as the progenitor of SN 1997D (Turatto et al. 1998), then the entire helium core of such a star ($`>`$8 M) could fall back and still enough nickel would escape to power the observed light curve! Bear in mind, however, that weaker explosions may well produce less mixing, so adapting the results of our simulations to supernovae like 1997D must be taken with some caution.

## 3 High Energy Spectral Calculations

For our spectral calculations, we have used data from the 3-dimensional explosion simulation discussed in the previous section. We input ejecta material properties from five different snapshots in time at 150, 200, 250, 300 and 365 days after explosion. Spectral calculations were carried out for both the Jet2 and Symmetric explosion models. Our analysis of these model spectra concentrates on the differences in total luminosity and line profile shape with the introduction of realistic explosion asymmetries. Since the progenitor star used as input to our simulations was a 15 Mred supergiant, we are unable to directly compare our calculated spectra with the observed high energy spectra of SN 1987A. However, we discuss how our models compare to various spectral trends observed from SN 1987A.

### 3.1 Numerical Schemes

We used a Monte Carlo technique, similar to that described in Ambwani & Sutherland 1988, for modelling $`\gamma `$-ray transport in 3-dimensions. Input models of the supernova ejecta (element abundances, density and velocities) were taken from the “Jet2” and “Symmetric” SPH explosion simulations and mapped onto a 140 $`\times `$140 $`\times `$140 cartesian grid. Escaping photons were tallied into 250 coarse energy bins, with finer binning at the decay line energies to provide line profile information. The emergent photons were also tallied into 11 angular bins ($`\mathrm{\Delta }\theta `$ = 10) along the polar axis (the models investigated in this work are axisymmetric, alleviating the need to tally in azimuthal angle as well.)

The decay of the radioactive species (predominantly 56Ni and its decay product: 56Co) in the supernova ejecta gives rise to the $`\gamma `$-ray line emission. As in Ambwani & Sutherland 1988, we assign the energy of the emitted photon packets according to the decay probabilities (Lederer & Shirley 1978) for the various radioactive species (56Ni, 56Co, 57Co, 44Ti, 44Sc, and 22Na). We include a total of 56 decay lines from these species, but for the explosion times considered, the packets fall predominantly into $``$ 15 decay lines. Roughly $`10^9`$ Monte Carlo photon bundles were generated for each input model in proportion to the mass of radioactive material distributed throughout the ejecta. The material properties of the ejecta were not evolved with photon flight time. However, we found that 99 % of the photons contributing to the observed model spectra have escape times of less than 2 days. This is sufficiently shorter than the timescales for change in the hydrodynamic models that our assumption of a fixed material background should be valid for the time slices considered here. The luminosity weight of each photon packet and the opacities seen by the packet were calculated in the comoving frame of the fluid, but all photon properties were boosted to the observer’s frame before being tallied into spectral observables.

In all models, photoelectric and pair production opacities were calculated for the elements H, He, C, N, O, Ne, Mg, Si, S, Ar, Ca, Ti, Cr, Fe, Co, and Ni which correspond to the elements used in the nucleosynthesis calculations for the progenitor star from Weaver & Woosley (1993). The cross section data for these elements were taken from the LLNL Evaluated Nuclear Data Library (Plechaty, Cullen, and Howerton 1981, revised 1987). The angle- and energy-dependent Compton scattering opacities were calculated assuming that all electrons, bound and free, contribute to the total cross section. Daughter products from the absorption processes were not followed (e.g. positron annihilation photons from pair production and K-shell fluorescence photons from photoelectric absorption.) In order to test the validity of this last approximation, we compared the results of our 3-dimensional Monte Carlo transport code (Maverick) with the 1-dimensional code FASTGAM (Pinto & Woosley 1988a). Good agreement was found between the model spectra whether the daughter products of absorption were included in FASTGAM or not. This justified the decision to disregard the fluorescence and annihilation photons in Maverick.

Figure 8, shows comparison spectra with the 1-dimensional FASTGAM code for model 10HMM (Pinto & Woosley 1988b). This is a logarithmic plot of photon flux (photons/second/MeV/cm2, assuming a distance of 60 kpc) across the energy range 0.3 keV - 4 MeV. The input model (10HMM) is an artificially mixed version of the Woosley 1988 10H model progenitor for SN 1987A. The first comparisons between FASTGAM and Maverick (left panel of Figure 8) showed significant differences in the location of the hard X-ray fall off. This was due to an invalid assumption in FASTGAM of constant electron fraction for all species which contribute to photoelectric and pair production opacities. This assumption is not valid for hydrogen, and results in a lower calculated opacity for the lower energy X-rays. The right panel of Figure 8 shows the hard X- and $`\gamma `$-ray spectra calculated from the two different codes after the absorption opacity correction was made. There is good agreement across the spectrum to within the uncertainties of the Monte Carlo calculation.

### 3.2 Hard X-ray and $`\gamma `$-ray Spectrum

Figure 9 is a logarithmic plot of photon flux in units of photons/second/MeV/cm2 across the energy range investigated with these simulations (0.3 keV - 4 MeV). We have placed this object at the distance of the Large Magellanic Cloud (60 kpc) for easy comparison with flux data from SN 1987A observations. The 5 panels are spectra from the different time slices; in each panel, we plot the spectrum for the Symmetric model, along with polar and equatorial views of the Jet2 model. It can be seen immediately that the hard X-rays emerge earlier from the ejecta with a global explosion asymmetry (Jet2 model). This holds regardless of viewing angle (pole versus equator) towards the explosion.

The fact that the hard X-ray flux in the aspherical explosion model is larger than the symmetric explosion, regardless of line of sight, can be understood in principle from optical depth arguments. In Figure 10, we show a contour plot of density (outer contour) and 56Co number density (inner contour) for the Jet2 and Symmetric models at t = 150 days. Decay of 56Co is the major source of $`\gamma `$-ray photons, so the inner contour essentially traces the surface of the emission region. The horizontal and vertical lines in Figure 10 represent lines of sight from the ejecta surface to the emission source and are labeled with the optical depth along that line of sight. The dominant opacity for the hard X- and $`\gamma `$-rays is Compton scattering off electrons and, since the density contours remain roughly spherical in both models, the optical depth from a given point to the ejecta surface is roughly constant. In the Jet2 model, the 56Ni was mixed out to larger radii in the polar direction, so it makes perfect sense that we see enhanced emission over the Symmetric model spectrum for that viewing angle (the optical depth that the high energy photons must pass through is 10 in the symmetric model versus $`67`$ along the polar line-of-sight in the Jet2 model).

At a first glance, one might expect that, in the Jet2 explosion along the equatorial line-of-sight, the total flux should also be low (the optical depth from the nickel in the equator is also roughly 10). However, this material does not dominate the high energy emission seen along the equatorial line-of-sight. The material ejected along the poles has been mixed far enough out in the ejecta that the optical depth these high energy photons must travel through, even along the equatorial line-of-sight, is quite low ($`7`$). It is this nickel which dominates the hard X-ray emission at all viewing angles. In fact, the optical depth from the “ends” of the 56Co distribution does not differ very much between the polar view and the equator view ($`\tau `$ = 6 & 7 respectively), which explains why the overall hard X-ray flux depends only mildly on viewing angle.

For the later time slices, this mismatch in escaping emission from the “ends” versus the equatorial plane ejecta is less pronounced, and the equator view spectrum has comparable contributions from both regions.

### 3.3 $`\gamma `$-ray Line Profiles

Although the overall hard X-ray emission shows little variation between pole and equator views, a detailed look at the $`\gamma `$-ray line profile shapes and strengths, for the 1.238 and 0.847 MeV 56Co lines, does reveal trends with viewing angle. Figures 11 and 12 show line profiles of these two 56Co lines for both the Symmetric and Jet2 explosion models. The broadening of the line is caused by Doppler velocity shifts resulting from the spatial distribution of radioactive nickel in the homologously expanding ejecta. The 4 panels shown are for days 200, 250, 300, and 365 after explosion. The three lines in the Jet2 spectra represent different viewing angles through the ejecta (along the pole, the equator and an intermediate angle $`45^{}`$.) For the Symmetric spectra, we have plotted these same viewing angles.

#### 3.3.1 Global Asymmetry

Both explosion scenarios (Symmetric and Jet2) show blue-shifted line profiles and this effect is most enhanced in the Jet2 model. These differences can be best understood by examining the physical effects which dictate the formation of the line profile edges. The blue edge to the lines is set by the maximum observed line-of-sight velocity of the 56Co in the ejecta. Because the nickel was mixed further out (and achieves higher velocities) for the asymmetric explosions, the $`\gamma `$-ray line profiles viewed along the line-of-sight of the jet are much more blueshifted for the Jet2 model than the Symmetric model

Note, however, that the line profiles in the asymmetric explosion model (Jet2) depend upon viewing angle. $`90^{}`$ off the jet axis, the nickel producing the line profiles arises from both fast-moving jet ejecta and the ejecta along the equator. Because the fast-moving ejecta is moving perpendicular to the line of sight, the $`\gamma `$-rays from this material are not blueshifted significantly, and the slow moving equatorial material also produces only a modest blueshift. Hence, the line profiles viewed along the equatorial line-of-sight are blueshifted much less than those viewed along the polar line-of-sight.

The red edge of the lines is determined by the escaping emission from 56Co with the smallest line-of-sight velocity in the ejecta. In a Symmetric model, this should be an indication of how deep into the ejecta we can see along a given viewing angle. However, there is a more pronounced effect at play in the asymmetric explosion models. As we mentioned in the previous section, much of the $`\gamma `$-ray emission for the equatorial view arises from the “ends” of the 56Co distribution. This material has a very low line-of-sight velocity for an equatorial observer, since it is being ejected predominantly in the polar direction. This allows for a significantly lower velocity red edge of the equator view lines, even though the optical depth profiles do not vary much between polar and equator viewing angles.

Also note that the line centroids become less blue-shifted with time. As the supernova expands, the opacity in the ejecta drops and emission from material located deeper into the ejecta (and thus at smaller radial velocities) becomes visible. This results in an overall broadening of the line, as well as a redward shift of the line centroid.

#### 3.3.2 Clumping Asymmetry

This basic understanding of the line profile edges addresses only global features of the emitted decay lines. To understand the range of possibilities that the line shape may take between those edges, we must concentrate on the spatial distribution of the radioactive elements in the supernova ejecta.

Because the expansion is basically homologous after 100 days, the line-of-sight velocity of a fluid element in the ejecta is proportional to its distance above the mid-plane of the explosion. Each spectral energy bin in the line profile can be mapped to a unique line-of-sight velocity in the ejecta, which can in turn be mapped to a specific height above the mid-plane. So, given the distribution of cobalt as a function of x, y and z-coordinates ( $`𝒩_{\mathrm{Co}}(x,y,z)`$ , where z is along the line-of-sight to the observer), the line profile shape should be proportional to

$$𝐍_{\mathrm{Co}}(z)=\underset{x,y}{}𝒩_{\mathrm{Co}}(x,y,z).$$
(3)
Optical depth effects, as discussed above, will alter the line shape from this basic picture, but it is important to realize that the clumping, which occurs as the inner material is mixed outward, can cause sizeable irregularities in the profiles. Figure 13 shows a plot of $`𝐍_{\mathrm{Co}}`$ versus line-of-sight velocity for the polar view of model Jet2 at 250 days with higher spatial resolution than the previous run. The emergent $`\gamma `$-ray line flux is plotted with a dotted line and can be seen to peak at velocities of around 1800 km/s. This corresponds to the location in the ejecta at which the exponential attenuation from the optical depth surpasses the flux enhancement from increasing cobalt mass: $`\tau _{MeV}4`$. As the ejecta expand and optical depth drops, the $`\gamma `$-ray line profiles will take on the irregular shape seen in the cobalt distribution.

## 4 Conclusions

In our simulations, we found that asymmetric explosions lead to extensive mixing of the supernova ejecta, placing the products of explosive burning well into the helium layer of the star. Even mild 2:1 asymmetries can mix 10% of the nickel out to the inner edge of the hydrogen layer. If such mixing occurred in weak explosions, as well as the strong explosions presented in this paper, these mildly asymmetric supernova explosions could explain the extensive mixing required in population III stars (Umeda & Nomoto 2002) and black hole binary systems such as Nova Scorpii (Podsiadlowski et al. 2002).

From the Monte Carlo transport simulations, we find that the high energy fluxes at early times ($``$150 days) for the asymmetric explosion model is roughly 4-5 times larger than for the corresponding symmetric explosion model. This suggests that emission from an explosion with mixing and global asymmetry would be observed roughly 50-100 days earlier than a symmetric explosion model with mixing alone. In addition, the line profiles for the asymmetric explosion vary with viewing angle and are brightest for the polar view. Along this viewing angle, the integrated line flux is enhanced by a factor of 2-3 over the symmetric explosion, and the lines peak at roughly 3-5 times brighter than the symmetric explosion line profiles. The integrated flux enhancement for the equator view of the Jet2 model is of order 30 % over the symmetric model. The line centroids observed for the asymmetric explosion are shifted relative to those from the symmetric model (more to the blue for the polar view and slightly more to the red for the equator view.) At later epochs in the SN explosion, the line profiles should reflect the irregular shapes seen in the underlying cobalt distribution, which result from the clumpy nature of the outward mixing of nickel in the ejecta. Table 2. summarizes the calculated luminosities of the different models for various continuum energy bands and the 1.238 and 0.847 MeV 56Co lines.

With the recent launch of the International Gamma-Ray Astrophysics Laboratory (INTEGRAL), it is particularly interesting to analyze our results in terms of the instrumental resolution and sensitivities this satellite can provide across these energy bands. It is important to keep in mind that the nickel mass synthesized in this 15 M model is roughly 3 times larger than the typically observed value ($``$0.08 M) for core-collapse SN explosions. Since the high energy emission should scale linearly with nickel mass, the fluxes observed for this model will be roughly 3 times brighter than we would expect for a more typical core-collapse event. In fact, the high energy continuum for our Symmetric model explosion is of order 3 times brighter than was observed towards SN 1987A. Since 1987A’s nickel yield was fairly typical, this comparison seems to validate the linear scaling of flux with nickel mass.

Our models suggest that $`\gamma `$-ray line observations will be most useful for diagnosing the departure of a core-collapse explosion from spherical symmetry. At energies around 1 MeV, INTEGRAL will have a spectral resolution of 2 keV and a narrow line sensitivity ($`3\sigma `$ in 106 seconds) of $``$5$`\times `$10-6 phot s-1 cm-2 (Hermsen & Winkler 2002). Our model lines are about 5 times broader than this resolution element, so the sensitivity for detecting them is worse by roughly $`\sqrt{5}`$. Using these specifications, INTEGRAL would be able to detect the 56Co lines from the polar view of the Jet2 model at a distance of $``$650 kpc.

We would like to thank P. Pinto, T. Evans, T. Urbatsch, A. Heger, K. Nomoto, S. Woosley, J. Cannizzo and N. Gehrels for their advice and encouragement on this paper. The large hydrodynamical simulations were run on the Avalon cluster at LANL. This work was performed under the auspices of the U.S. Department of Energy by Los Alamos National Laboratory under contract W-7405-ENG-36. Funding was also provided by a Feynman Fellowship at LANL, DOE SciDAC grant number DE-FC02-01ER41176 and Beomax, Inc. These simulations were run on the Avalon and Space-Simulator clusters at LANL.


format: text

text:
# On the existence of a new family of Diophantine equations for 𝛀

## 1. Introduction and Motivation

In his 1975 paper , Gregory Chaitin introduced the number $`\mathrm{\Omega }`$. $`\mathrm{\Omega }`$ is a real number between 0 and 1 which is the *halting probability*: the probability that a randomly chosen program will halt. In order to make this precise, one must choose a particular (recursive) programming language in which the programs are self-delimiting binary strings. This allows us to choose a random program by flipping an unbiased coin to generate each bit, stopping when we reach a valid (self-delimited) program. This method makes longer programs increasingly unlikely and allows $`\mathrm{\Omega }`$ to be a well defined probability. The exact value of $`\mathrm{\Omega }`$ obviously depends upon which programming language we choose, and thus Chaitin  specifies a particular language to make talk of $`\mathrm{\Omega }`$ concrete. However, its important properties are the same regardless and the exact language chosen need not concern us.111The properties possessed by $`\mathrm{\Omega }`$ are numerous and we can only discuss the most pertinent here. For further information on $`\mathrm{\Omega }`$’s exotic attributes, see .

It is well known that if we had access to $`\mathrm{\Omega }`$ we could solve the halting problem . Furthermore, if we had some method of solving the halting problem then we could compute $`\mathrm{\Omega }`$ to any desired accuracy. In this sense there is an equivalence between $`\mathrm{\Omega }`$ and the halting problem: the ability to compute one gives us the ability to compute the other.

It is also well known that both the halting problem and $`\mathrm{\Omega }`$ have analogues in number theory. The halting problem has been shown to be equivalent to the problem of determining whether a given Diophantine equation has solutions, while $`\mathrm{\Omega }`$ has been shown to be equivalent to determining whether a finite or infinite amount of Diophantine equations of a particular type have solutions. However, unlike their counterparts in the theory of computability, these number theoretic versions have no direct connection between them. In this paper, we present a new representation of $`\mathrm{\Omega }`$ within number theory that revolves around a matter of parity rather than finitude and makes clear the link between the number theoretic analogues of $`\mathrm{\Omega }`$ and the halting problem.

In Sections 2 and 3 we provide the formal definitions of $`\mathrm{\Omega }`$ and a related real number $`\tau `$, presenting an efficient method for using the solutions to the halting problem to determine $`\mathrm{\Omega }`$. In Sections 4 and 5, we then introduce Diophantine equations and review the known method of using them to represent $`\mathrm{\Omega }`$, providing proofs when they will facilitate the proofs of our new results. In Sections 6 and 7 we present and prove our results for a new representation of $`\mathrm{\Omega }`$ and discuss their implications.

## 2. Algorithmic Randomness and $`𝛀`$

###### Definition 2.1.

$`\mathrm{\Omega }`$ is formally defined by the following equation, where $`p`$ ranges over all programs in the language specified by Chaitin  and $`\left|p\right|`$ is the number of digits in the binary representation of $`p`$.

(2.1)
$$\mathrm{\Omega }=\underset{p\text{ halts}}{}2^{\left|p\right|}$$

The importance of $`\mathrm{\Omega }`$ lies in its contrasting properties of being both an *algorithmically random* and *recursively enumerable* real number.

As a random real, $`\mathrm{\Omega }`$ is uncomputable in a very strong sense. If we look at the binary expansion of $`\mathrm{\Omega }`$ (the infinite sequence of 1’s and 0’s following the binary point), this sequence is highly incompressible. To see how this is so, let us compare it with the number $`\tau `$ from Copeland and Proudfoot , a non-recursive real that is not random.

###### Definition 2.2.

$`\tau `$ is the real number between 0 and 1 whose $`n`$-th binary digit is 1 if the $`n`$-th program $`p_n`$ halts when given no input and 0 otherwise. In order to pick out a particular ordering of programs, let each program have a binary representation as above and let them be ordered in the usual lexical ordering for binary strings. We could also write:

(2.2)
$$\tau =\underset{p_n\text{ halts}}{}2^n$$

$`\tau `$ in an encoding of the answers to every instance of the halting problem in a single real number. While each bit of $`\tau `$ tells us whether a particular program halts when run with no input, we can easily use the bits of $`\tau `$ to answer the more common formulation of halting problem: ‘does program $`p`$ halt when run on input $`i`$?’. Given $`p`$ and $`i`$, we can (recursively) construct a new program which takes no input and simulates the computation of $`p`$ on $`i`$. The bit of $`\tau `$ corresponding to this new program tells us whether or not $`p`$ halts on $`i`$.

The concept of compressibility for an infinite string can be understood through examining the minimum amount of advice needed to find out $`n`$ bits of that string . For recursive infinite strings (such as the binary digits of $`\pi `$), we can perform massive compression by asking for the bits of a particular program which will generate every bit of $`\pi `$, one by one. In this sense, there are only finitely many bits of information (those of that particular program) in the infinite bitstring that makes up $`\pi `$.

On the other hand, if we wish to determine the bits of $`\tau `$, there is no program that will do this, since the halting problem is undecidable. Therefore, determining $`n`$ bits of $`\tau `$ will require an amount of advice that increases with the value of $`n`$. Clearly we could find out $`n`$ bits of $`\tau `$ with $`n`$ bits of advice as we could be directly given those bits. However, Chaitin  shows how we can do better than this, using only $`\mathrm{log}_2n`$ bits. This is because the $`n`$ bits of $`\tau `$ represent $`n`$ instances of the halting problem and to solve these, we only need to know how many of the $`n`$ programs halt. We can then simulate each program in parallel until this many have halted, being confident that every one that is still running will never halt. Therefore, while the (Turing) non-computability of $`\tau `$ makes it globally incompressible (computing infinitely many bits of $`\tau `$ requires infinitely many bits of advice), it is still locally compressible as $`n`$ bits of $`\tau `$ are computable from $`\mathrm{log}_2n`$ bits of advice.

In contrast, Chaitin  defines a random real as one for which calculating $`n`$ bits of its binary expansion requires more than $`n`$ bits of advice. The reason that *more* than $`n`$ bits are needed is that the advice, just like the program itself, must be self-delimiting. Chaitin  has shown that $`\mathrm{\Omega }`$ satisfies this condition and is thus a random real. In this way, the first $`n`$ bits of $`\mathrm{\Omega }`$ contain $`n`$ bits of *algorithmically incompressible* information.

In addition to recursive incompressibility, random reals are also characterised by *recursive unpredictability* . Consider a ‘predictive’ program that takes a finite initial segment of an infinite bitstring and returns a value indicating either ‘the next bit is 1’, ‘the next bit is 0’ or ‘no prediction’. If any such program is run on all finite prefixes of the binary expansion of a random real and makes an infinite amount of predictions, the limiting relative frequency of correct predictions approaches $`\frac{1}{2}`$. In other words when any program is used to predict infinitely many bits of a random real, such as $`\mathrm{\Omega }`$, it does no better than random — even with information about all the prior bits.

In contrast to these results, $`\mathrm{\Omega }`$ is a recursively enumerable (r.e.) real. For a mathematical object to be r.e. it need not be computable, but there must be a certain method of successively approximating it. For a set of positive integers to be r.e. there must be a program that halts on $`i`$ if and only if $`i`$ is in the set. For an infinite binary sequence to be r.e. there must be a program that halts on $`i`$ if and only if the $`i`$-th bit of the sequence is 1. For a real number $`x`$ to be r.e. there must be a program that takes a positive integer $`i`$ and generates a rational approximation $`x_i`$ where $`\left\{x_i\right\}`$ is an increasing sequence of rationals that converges to $`x`$. In the case of $`\mathrm{\Omega }`$, we can construct such a sequence which we shall denote $`\left\{\mathrm{\Omega }_i\right\}`$, where

(2.3)
$$\mathrm{\Omega }_i=\underset{\stackrel{\left|p\right|i}{p\text{ halts in }i\text{ steps}}}{}2^{\left|p\right|}$$

However, it is also important to ask whether the binary expansion of $`\mathrm{\Omega }`$ is an r.e. bitstring. For $`\tau `$, this is the case because there is a program that halts on $`i`$ if and only if the $`i`$-th bit is 1 — indeed this is just a universal program, which simulates the $`i`$-th program (and thus halts iff the $`i`$-th program does). For $`\mathrm{\Omega }`$ however, there is no such program. This is a rather surprising, but crucial, fact about $`\mathrm{\Omega }`$.

We can see that this must be true given the incompressibility of $`\mathrm{\Omega }`$. If there were a program to enumerate the bits of $`\mathrm{\Omega }`$, we could use the trick described above for compressing $`\tau `$: we could determine $`n`$ bits of $`\mathrm{\Omega }`$ by asking how many of these $`n`$ values make this program halt (and thereby require only $`\mathrm{log}_2n`$ bits). Therefore, as $`\mathrm{\Omega }`$ is incompressible, we see that there is indeed no such program for $`\mathrm{\Omega }`$ and that its sequence of bits is not r.e. Because all random reals share the property of being incompressible, this argument carries over and no random real can have an r.e. sequence of bits.

In the remainder of the paper, we shall use the fact that $`\mathrm{\Omega }`$ is an r.e. real to express it through Diophantine equations and thereby show how algorithmic randomness occurs even in the heart of number theory.

## 3. Computing $`𝛀`$ from $`\tau `$

Having seen that there are interesting connections between $`\mathrm{\Omega }`$ and $`\tau `$, it is natural to ask whether we can compute $`\mathrm{\Omega }`$ given access to $`\tau `$ (or in the language of recursion theory, whether $`\mathrm{\Omega }`$ is *Turing reducible* to $`\tau `$). One way to go about this is to construct a program $`Q`$ that takes two positive integers, $`k`$ and $`N`$, as input and halts iff $`\frac{N}{2^k}<\mathrm{\Omega }`$. Such a program is quite easy to construct, it just needs to enumerate $`\left\{\mathrm{\Omega }_i\right\}`$ (as done in the previous section) and check at each stage whether $`\frac{N}{2^k}<\mathrm{\Omega }_i`$, halting if this is true and continuing otherwise. Since $`\left\{\mathrm{\Omega }_i\right\}`$ approaches $`\mathrm{\Omega }`$ from below, if $`\frac{N}{2^k}<\mathrm{\Omega }`$ then there is some $`i`$ for which $`\frac{N}{2^k}<\mathrm{\Omega }_i`$ and $`Q`$ halts as required. On the other hand, if $`\frac{N}{2^k}\mathrm{\Omega }`$ then there is no such $`i`$ and $`Q`$ will not halt. From $`\tau `$, we can determine whether or not $`Q`$ halts on a given input, and thus whether or not each $`\frac{N}{2^k}`$ is less than $`\mathrm{\Omega }`$.

To determine the first $`k`$ bits of $`\mathrm{\Omega }`$, we just need to determine the greatest value of $`N`$ for which $`\frac{N}{2^k}<\mathrm{\Omega }`$. We will then know that for this $`N`$, $`\frac{N}{2^k}<\mathrm{\Omega }<\frac{N+1}{2^k}`$ and thus that if $`N`$ is expressed as a $`k`$ digit binary number (with leading zeros if required), it will be the first $`k`$ digits of $`\mathrm{\Omega }`$.

We can therefore use $`\tau `$ to determine the first $`k`$ bits of $`\mathrm{\Omega }`$, by checking whether $`Q`$ halts when applied to $`k`$ and $`N`$ for each value of $`N`$ from 1 to $`2^k1`$. This technique allows us to determine the first $`k`$ bits of $`\mathrm{\Omega }`$ from $`2^k1`$ carefully chosen bits of $`\tau `$. Indeed, given these bits, it is so simple to compute the first $`k`$ bits of $`\mathrm{\Omega }`$, that we could do it by simply combining them with a truth table and $`\mathrm{\Omega }`$ is not only Turing reducible to $`\tau `$, but also *truth table reducible*.

However, if we are prepared to sacrifice this extreme algorithmic simplicity, we can get by with even fewer bits of $`\tau `$ by using a *bisection search*. In the bits of $`\tau `$ that we need we know there is a lot of structure: if $`\frac{N}{2^k}<\mathrm{\Omega }`$ then for all $`i<N`$, $`\frac{i}{2^k}<\mathrm{\Omega }`$. Similarly, if $`\frac{N}{2^k}\mathrm{\Omega }`$ then for all $`i>N`$, $`\frac{i}{2^k}\mathrm{\Omega }`$. Thus, all we need to find is the greatest value of $`N`$ for which $`\frac{N}{2^k}<\mathrm{\Omega }`$. This can be done by first trying $`N=\frac{2^k}{2}`$ and checking if it is less than or greater than $`\mathrm{\Omega }`$. If it is greater, we then try the midpoint of 0 and $`\frac{2^k}{2}`$. If it is less, we try the midpoint of $`\frac{2^k}{2}`$ and $`2^k`$. We proceed in this manner until finding the value of $`N`$ for which $`\frac{N}{2^k}<\mathrm{\Omega }`$ and $`\frac{N+1}{2^k}\mathrm{\Omega }`$.

Performing this type of bisection search on $`n`$ items requires only $`\mathrm{log}_2(n+1)`$ queries, rounding up to the nearest integer. As we do not need to test for $`N=0`$ or $`N=2^k`$, there are only $`2^k1`$ remaining values to check and we can therefore get the first $`k`$ bits of $`\mathrm{\Omega }`$ with only $`k`$ bits of $`\tau `$. From the incompressibility of $`\mathrm{\Omega }`$, we can see that this must be optimal: no less than $`k`$ bits of any other real whatsoever can give us $`k`$ bits of $`\mathrm{\Omega }`$.222It may seem like there is a contradiction here because we could use the compression trick for $`\tau `$ again, giving us a method for computing $`k`$ bits of $`\mathrm{\Omega }`$ from only $`\mathrm{log}_2k`$ bits of advice. However, this turns out to be impossible because to use the bisection search, we need to find out the value of one bit of $`\tau `$ before we know which bit to ask for next. There is no time at which we know in advance a complete set of $`k`$ bits of $`\tau `$ that we can ask for in compressed form.

Having shown a method for determining the bits of $`\mathrm{\Omega }`$ from an r.e. sequence of bits ($`\tau `$), we will now examine Diophantine equations and look at how this method can be used to find $`\mathrm{\Omega }`$ in number theory.

## 4. Diophantine Equations and Hilbert’s Tenth Problem

###### Definition 4.1.

A *Diophantine equation* is an equation of the form

(4.1)
$$D(x_1,\mathrm{},x_m)=0$$

where $`D`$ is a polynomial with integer coefficients. The variables are typically allowed to range over the integers or the non-negative integers, but here it will be convenient to restrict them to the positive integers.

The study of Diophantine equations is a central area of number theory. It is well known that some Diophantine equations have solutions (such as $`2x_14=0`$) and some do not (such as $`2x_13=0`$). The task of devising an algorithm for determining whether or not an arbitrary Diophantine equation is solvable is known as Hilbert’s tenth problem and has been a major area of research in 20th Century mathematics (see Matiyasevich ).

###### Definition 4.2.

A *family of Diophantine equations* is a relation of the form

(4.2)
$$D(a_1,\mathrm{},a_n,x_1,\mathrm{},x_m)=0$$

where we distinguish between the variables $`a_1,\mathrm{},a_n`$ which are called parameters and $`x_1,\mathrm{},x_m`$ which are called unknowns. Fixing values of the parameters, specifies one of the individual Diophantine equations that comprise the family.

###### Definition 4.3.

A set of $`n`$-tuples, $`𝒟`$, is *Diophantine* iff it can be represented by a family of Diophantine equations in the following way

(4.3)
$$a_1,\mathrm{},a_n𝒟x_1\mathrm{}x_mD(a_1,\mathrm{},a_n,x_1,\mathrm{},x_m)=0$$

Many sets, such as the square numbers or the pairs $`a,b`$ such that $`a`$ divides $`b`$ are Diophantine. The question of exactly which sets are Diophantine has received considerable attention and was solved through the work of Martin Davis, Hilary Putnam, Julia Robinson and Yuri Matiyasevich .

###### Theorem 4.4 (dprm Theorem).

Every r.e. set $`𝒟`$ of $`n`$-tuples of positive integers is Diophantine, that is:

(4.4)
$$a_1,\mathrm{},a_n𝒟x_1\mathrm{}x_mD(a_1,\mathrm{},a_n,x_1,\mathrm{},x_m)=0$$

###### Proof.

See . ∎

The dprm Theorem provides an excellent tool for showing the existence of Diophantine equations with certain properties. All that is needed is to come up with a set that has the desired property and prove that it is r.e. For example, the following is immediately evident:

###### Corollary 4.5 (Family of Diophantine equations for $`\tau `$).

There is a family of Diophantine equations:

(4.5)
$$T(n,x_1,\mathrm{},x_m)=0$$

which has solutions for a given value of $`n`$ iff the $`n`$-th bit of $`\tau `$ is 1 or, equivalently, iff the $`n`$-th program halts.

###### Proof.

The set of positive integers $`n`$ such that the $`n`$-th bit of $`\tau `$ is 1 is clearly an r.e. set. Thus, by the dprm Theorem, there is a family of Diophantine equations with a parameter, $`n`$, which has a solution for a given value of $`n`$ iff the $`n`$-th bit of $`\tau `$ is 1. ∎

From this it is clear that Hilbert’s tenth problem must be recursively undecidable. No program could decide whether a given Diophantine equation has a solution because this would allow a program to compute the bits of $`\tau `$ — a task that is known to be non-recursive.333Indeed, it was long known that the recursive undecidability of Hilbert’s Tenth Problem would follow immediately from the dprm Theorem and this was the main motivation for its proof — the Diophantine representations for all other r.e. sets being a nice corollary. We therefore have undecidability in number theory, with direct analogues of the halting problem and $`\tau `$. In what follows we show how we can also find algorithmic randomness and, in doing so, add $`\mathrm{\Omega }`$ to this list.

## 5. Chaitin’s Expression of $`𝛀`$ Through Diophantine Equations

From the dprm Theorem, we can see that $`\mathrm{\Omega }`$ cannot be directly represented by a Diophantine set (one that has $`k`$ as a member iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1) because the bits of $`\mathrm{\Omega }`$ are not r.e. However, taking a slightly different approach, Chaitin  shows that the bits of $`\mathrm{\Omega }`$ can be found in another property of Diophantine equations.

###### Theorem 5.1 (Family of Diophantine Equations for $`\mathrm{\Omega }`$ via Finitude).

There is a family of Diophantine equations

(5.1)
$$\chi _1(k,N,x_1,\mathrm{},x_m)=0$$

such that for a given $`k`$ there are infinitely many values of $`N`$ for which it has a solution if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1 and finitely many values of $`N`$ for which it has a solution if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 0.

###### Proof.

There is a program $`P`$ that takes two inputs, $`k`$ and $`N`$, generating $`\mathrm{\Omega }_N`$ and returning the value of its $`k`$-th bit. As $`N`$ increases, the $`k`$-th bit of $`\mathrm{\Omega }_N`$ may change between 0 and 1 many times, but since $`\left\{\mathrm{\Omega }_N\right\}`$ is monotonically increasing, each bit can only change finitely many times (at most $`2^k1`$). Therefore, the $`k`$-th bit of $`\mathrm{\Omega }_N`$ will eventually settle on either a 0 or a 1 and because $`\left\{\mathrm{\Omega }_N\right\}`$ converges to $`\mathrm{\Omega }`$, this final value must be the same as the $`k`$-th bit of $`\mathrm{\Omega }`$. Thus, the $`k`$-th bit of $`\mathrm{\Omega }_N`$ will only differ from the $`k`$-th bit of $`\mathrm{\Omega }`$ for finitely many values of $`N`$.

Let $`𝒞_1`$ be the set of pairs $`k,N`$ such that $`P`$ returns 1 when applied to $`k`$ and $`N`$. If the $`k`$-th bit of $`\mathrm{\Omega }`$ is 0, $`P`$ will return 1 (the correct value for $`\mathrm{\Omega }_N`$ but incorrect for $`\mathrm{\Omega }`$) for only finitely many values of $`N`$. On the other hand, if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1, $`P`$ will return 1 for infinitely many values of $`N`$. In other words, for a given $`k`$, there are infinitely many pairs $`k,N𝒞_1`$ if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1 and finitely many if it is 0. Since $`𝒞_1`$ is recursive, it is r.e. and thus via the dprm Theorem there is a Diophantine equation, $`\chi _1(k,N,x_1,\mathrm{},x_m)=0`$, that has a solution iff $`k,N𝒞_1`$. Thus, for a given value of $`k`$, this equation has solutions for infinitely many values of $`N`$ iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1. ∎

Chaitin takes this result further and provides a simpler example of $`\mathrm{\Omega }`$ occurring in number theory. In doing so he moves to *exponential Diophantine equations*. Where the polynomial in a Diophantine equation consists of variables and integer constants composed together with addition and multiplication, exponential Diophantine equations merely add exponentiation to this list, with the proviso that no negative constants can appear in the exponents. In proving the dprm Theorem, Matiyasevich  showed that all sets which are exponential Diophantine are Diophantine as well. However, Chaitin’s move to exponential Diophantine equations allows the use of another key result due to Matiyasevich  which concerns not just the existence of solutions to an equation, but the quantity of solutions as well.

###### Theorem 5.2 (Existence of Singlefold Exponential Diophantine Equations).

For any Diophantine set, $`𝒟`$, there is a family of exponential Diophantine equations

(5.2)
$$E(a_1,\mathrm{},a_n,x_1,\mathrm{},x_m)=0$$

that is a *singlefold* representation of $`𝒟`$: it has exactly one solution for $`a_1,\mathrm{},a_n𝒟`$ and no solutions otherwise.

###### Proof.

See . ∎

###### Theorem 5.3 (Family of Exponential Diophantine Equations for $`𝛀`$ via Finitude).

There is a family of exponential Diophantine equations

(5.3)
$$\chi _1^e(k,x_0,x_1,\mathrm{},x_m)=0$$

with infinitely many solutions for $`k`$ if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1 and finitely many if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 0.

###### Proof.

Let $`\chi _1^e(k,N,x_1,\mathrm{},x_m)=0`$ be the family of singlefold exponential Diophantine equations for $`\chi _1`$ of Theorem 5.1. The existence of $`\chi _1^e`$ is ensured by Theorem 5.2. Given a value of $`k`$, there is a single solution for each of infinitely many values of $`N`$ iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1. If we treat $`N`$ as another unknown instead of a parameter, we get $`\chi _1^e(k,x_0,x_1,\mathrm{},x_m)=0`$ which is a family of exponential Diophantine equations with infinitely many solutions for $`k`$ iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1. ∎

Thus, while the question of whether a Diophantine equation has solutions is undecidable in general, the question of whether an exponential Diophantine equation has infinitely many solutions is much worse. The task is no longer r.e. (since solving it gives the bits of $`\mathrm{\Omega }`$) and as a single parameter is varied, the results fluctuate in an algorithmically random manner: there is absolutely no recursive pattern to be found.

Chaitin  went even further than we have here by actually constructing an exponential Diophantine equation $`\chi _1^e(k,x_0,x_1,\mathrm{},x_m)=0`$. This equation was automatically generated from a complex register machine program and is very large, with approximately 17,000 unknowns. While it is has been shown that this can be reduced to just three , doing so would be a very challenging task.

## 6. A New Expression of $`𝛀`$ Through Diophantine Equations

Using the method of computing $`\mathrm{\Omega }`$ from $`\tau `$ discussed in Section 3 we can now present our main results. While the bits of $`\mathrm{\Omega }`$ occur in Theorems 5.1 and 5.3 through the distinction between the finite and the infinite, the presentation below remains within the realms of the finite, with the bits of $`\mathrm{\Omega }`$ occurring in the changes of parity.

###### Theorem 6.1 (Family of Diophantine Equations for $`𝛀`$ via Parity).

There is a family of Diophantine equations

(6.1)
$$\chi _2(k,N,x_1,\mathrm{},x_m)=0$$

such that for a given $`k`$:

* there are less than $`2^k`$ values of $`N`$, taken from $`\{1,\mathrm{},2^k1\}`$, for which there is a solution,
* and there are an odd number of values of $`N`$ for which there is a solution iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1.

###### Proof.

Let $`𝒞_2`$ be the set of pairs of positive integers $`k,N`$ such that $`\frac{N}{2^k}<\mathrm{\Omega }`$. From Section 3, we see that $`𝒞_2`$ is r.e. and so, by the dprm Theorem, there exists a Diophantine equation, $`\chi _2(k,N,x_1,\mathrm{},x_m)=0`$, which has a solution iff $`k,N𝒞_2`$. Since $`0<\mathrm{\Omega }<1`$, $`k,N𝒞_2`$ implies that $`N\{1,\mathrm{},2^k1\}`$, and so these are the only values of $`N`$ for which $`\chi _2(k,N,x_1,\mathrm{},x_m)=0`$ can have solutions.

If the greatest value of $`N`$ for which $`\frac{N}{2^k}<\mathrm{\Omega }`$ is odd, then the $`k`$-th bit of $`\mathrm{\Omega }`$ must be 1 and similarly, if the greatest value of $`N`$ for which $`\frac{N}{2^k}<\mathrm{\Omega }`$ is even, then the $`k`$-th bit of $`\mathrm{\Omega }`$ must be 0. As the greatest value of $`N`$ for which $`\frac{N}{2^k}<\mathrm{\Omega }`$ is equal to the number of values of $`N`$ for which $`\frac{N}{2^k}<\mathrm{\Omega }`$, it follows that for a given $`k`$, there are an odd number of values of $`N`$ for which $`\chi _2(k,N,x_1,\mathrm{},x_m)=0`$ has solutions iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1. ∎

Unlike the representation of $`\mathrm{\Omega }`$ through finitude, this representation would allow $`\mathrm{\Omega }`$ to be directly computed by determining whether certain Diophantine equations have solutions (i.e. by solving several instances of Hilbert’s tenth problem). This mirrors the relationship between $`\mathrm{\Omega }`$ and the halting problem in computability theory. Even the efficiency of the computation is preserved: the bisection approach could also be used here to determine the $`k`$-th bit of $`\mathrm{\Omega }`$ from the solutions to just $`k`$ instances of Hilbert’s tenth problem.

###### Theorem 6.2 (Family of Exponential Diophantine Equations for $`𝛀`$ via Parity).

There is a family of exponential Diophantine equations

(6.2)
$$\chi _2^e(k,x_0,x_1,\mathrm{},x_m)=0$$

which for a given $`k`$,

* has less than $`2^k`$ solutions, where $`x_0`$ takes distinct values from the set $`\{1,\mathrm{},2^k1\}`$,
* and has an odd number of solutions iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1.

###### Proof.

Let $`\chi _2^e(k,N,x_1,\mathrm{},x_m)=0`$ be the singlefold family of exponential Diophantine equations for $`𝒞_2`$. Given a value of $`k`$, there is a single solution for each of an odd number of values of $`N`$ iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1. If we treat $`N`$ as another unknown instead of as a parameter, we get $`\chi _2^e(k,x_0,x_1,\mathrm{},x_m)=0`$. For a given value of $`k`$, this has less than $`2^k`$ solutions and in these solutions, $`x_0`$ takes distinct values from $`\{1,\mathrm{},2^k1\}`$. Furthermore, the number of solutions is odd iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1. ∎

In addition to the above results, we can use both the finitude and parity methods to produce a new source of algorithmic randomness in number theory. Using the following lemma , we can design polynomial expressions with variables and a parameter $`k`$ ranging over the positive integers that exhibit randomness in the number of positive values they assume as $`k`$ is varied.

###### Lemma 6.3.

Given a family of Diophantine equations with two parameters

(6.3)
$$D(k,N,x_1,\mathrm{},x_m)=0$$

we can construct a polynomial with integer coefficients, one parameter and one additional variable $`x_0`$ as follows

(6.4)
$$P(k,x_0,x_1,\mathrm{},x_m)x_0\left(1(D(k,x_0,x_1,\mathrm{},x_m))^2\right).$$

If we restrict the variables to values in the positive integers, then the set of positive values that this polynomial assumes for a given $`k`$ is exactly the set of all $`N`$ such that $`D(k,N,x_1,\mathrm{},x_m)=0`$ has solutions.

###### Proof.

For all $`x_0,x_1,\mathrm{},x_m`$ such that $`D(k,x_0,x_1,\mathrm{},x_m)=0`$, the polynomial $`P(k,x_0,x_1,\mathrm{},x_m)`$ will assume the value $`x_0`$. For all $`x_0,x_1,\mathrm{},x_m`$ such that $`D(k,x_0,x_1,\mathrm{},x_m)0`$, $`P(k,x_0,x_1,\mathrm{},x_m)`$ will assume a non-positive value. Thus, $`P(k,x_0,x_1,\mathrm{},x_m)`$ will assume all values of $`N`$ such that $`D(k,N,x_1,\mathrm{},x_m)=0`$ has solutions and no other positive values. ∎

###### Theorem 6.4 (Polynomial Expression for $`𝛀`$ via Finitude).

(6.5)
$$P_1(k,x_0,x_1,\mathrm{},x_m)x_0\left(1(\chi _1(k,x_0,x_1,\mathrm{},x_m))^2\right)$$

is a polynomial with integer coefficients and a parameter $`k`$. If we restrict the variables to values in the positive integers, then this polynomial assumes infinitely many positive values iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1.

###### Proof.

By Lemma 6.3, the set of positive values that $`P_1`$ assumes for a given $`k`$ is precisely the set of $`N`$ such that $`\chi _1(k,N,x_1,\mathrm{},x_m)=0`$ has solutions. By Theorem 5.1, this set has infinitely many members iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1, so $`P_1`$ must assume a countable infinity of positive values iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1. ∎

###### Theorem 6.5 (Polynomial Expression for $`𝛀`$ via Parity).

(6.6)
$$P_2(k,x_0,x_1,\mathrm{},x_m)x_0\left(1(\chi _2(k,x_0,x_1,\mathrm{},x_m))^2\right)$$

is a polynomial with integer coefficients and a parameter $`k`$. If we restrict the variables to values in the positive integers, then

* this polynomial assumes less than $`2^k`$ values in the positive integers, taken from $`\{1,\mathrm{},2^k1\}`$,
* and the number of positive values that it takes is odd iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1.

###### Proof.

From Lemma 6.3, the set of positive integer values that $`P_2(k,x_0,x_1,\mathrm{},x_m)`$ assumes for a given $`k`$ is precisely the set of $`N`$ such that $`k,N𝒞_2`$. By Theorem 6.1, this set only includes values from $`\{1,\mathrm{},2^k1\}`$. Theorem 6.1 also tells us that there are an odd number of $`N`$ such that $`k,N𝒞_2`$ iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1, so $`P_2(k,x_0,x_1,\mathrm{},x_m)`$ must assume an odd number of positive values iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1. ∎

While Theorems 6.1, 6.2 and 6.5 are presented in terms of producing the $`k`$-th bit of $`\mathrm{\Omega }`$, the ideas behind them have a somewhat more general application. In each theorem, a set of values is discussed and its size is considered. For Theorem 6.1, this is the number of values of $`N`$ for which $`\chi _2=0`$ has solutions, for Theorem 6.2 it is the number of solutions for $`\chi _2^e=0`$ and for Theorem 6.5 it is the number of positive values taken by $`P_2`$. For a given value of $`k`$, each of these quantities will equal the same number which we shall call $`q_k`$.

Previously, we just looked at the parity of $`q_k`$ and used this single bit of information to determine a single bit of $`\mathrm{\Omega }`$. Since $`q_k`$ may take on $`2^k`$ different values, it contains $`k1`$ additional bits of information and these can directly provide all prior bits of $`\mathrm{\Omega }`$. When $`q_k`$ is expressed in binary (with enough leading zeros to give $`k`$ digits) it comprises the first $`k`$ bits of $`\mathrm{\Omega }`$. For example, if for $`k=5`$, there are six solutions to $`\chi _2^e=0`$, then the first five bits of $`\mathrm{\Omega }`$ are 00110.

In this way, we can see that the first $`k1`$ bits of $`q_k`$ are just the $`k1`$ bits of $`q_{k1}`$. So, given $`q_{k1}`$, there are only two possibilities for $`q_k`$: either $`q_k=2q_{k1}`$ or $`q_k=2q_{k1}+1`$. Thus, while the sequence $`\left\{q_kmod2\right\}`$ is algorithmically random, the sequence $`\left\{q_k\right\}`$ is highly structured.

We can also extend our results to the expression of $`\mathrm{\Omega }`$ in different bases. If we consider the representation of $`\mathrm{\Omega }`$ in base $`b`$, we can generate an equation for finding its $`k`$-th digit by replacing all references to $`2^k`$ with $`b^k`$. Instead of checking the parity of the appropriate quantity, the $`k`$-th digit is simply the remainder when this quantity is divided by $`b`$. Similarly, if we represent the quantity in base $`b`$ with enough leading zeros, it forms the first $`k`$ digits in a base $`b`$ representation of $`\mathrm{\Omega }`$ and $`q_k`$ must be one of $`\{bq_{k1},bq_{k1}+1,\mathrm{},bq_{k1}+(b1)\}`$.

## 7. Concluding Remarks

In this paper we report on our findings on the existence of a new family of Diophantine equations for determining the digits of $`\mathrm{\Omega }`$. Our equations differ from the previous families of Diophantine equations in several important features. Previously, a particular bit of $`\mathrm{\Omega }`$ was represented by whether some Diophantine equations have solutions for a finite or infinite number of values of a given parameter or, in the exponential Diophantine case, whether the number of solutions is finite. In contrast, our formulation always has a finite number of parameter values for which solutions could occur and a finite number of solutions in the case of the exponential Diophantine equations. The value of the corresponding bit of $`\mathrm{\Omega }`$ depends only on whether this number is even or odd. Thus, in the relatively mundane switching between an odd or an even number of solutions as a parameter is varied, the full subtlety of algorithmic randomness is felt: the incompressibility and the unpredictability.

These new families of Diophantine equations for $`\mathrm{\Omega }`$ also provide a symmetry to the relationship between randomness and undecidability in the fields of computability and number theory. Just as we can directly compute $`\mathrm{\Omega }`$ (in its guise as the halting probability) from solutions to instances of the halting problem, so can we compute $`\mathrm{\Omega }`$ (in its guise as a property of a family of Diophantine equations) from solutions to instances of Hilbert’s tenth problem. Furthermore, the translation into the domain of number theory preserves the efficiency of the computation, producing $`n`$ bits of $`\mathrm{\Omega }`$ from either the solutions to $`n`$ instances of the halting problem or the solutions to $`n`$ instances of Hilbert’s tenth problem.

Finally, we provide a further example of how algorithmic randomness occurs in number theory: a polynomial with positive integer variables and a parameter $`k`$ that takes on an odd number of positive values iff the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1. Along with the other results of this paper, this hints at the variety of ways in which this randomness can occur. While the polynomials and exponentials that give rise to algorithmic randomness are not of the type that are likely to occur in classical research in number theory, their importance lies in showing that in some places there are facts with no recursive pattern at all. This is not to say that these facts are *completely* patternless — on the contrary, we have shown intricate dependencies between the expressions of $`\mathrm{\Omega }`$ and those of the halting problem — but in the face of so subtle a pattern, computer programs can perform no better than coin tossing. Without some means to transcend the fundamental limits of our current computers, whether through some kind of mathematical insight or radical new technology444Examples and discussions of such non-classical computation have been presented by both authors  and a very thorough survey can be found in Copeland ., they will remain completely beyond our grasp.

## Acknowledgements

We would like to thank Gregory Chaitin for helpful discussion. TDK wishes to acknowledge the continuing support of Peter Hannaford.


text:
# Evolution of bimodal accretion flows

## 1 Introduction

The standard model for an optically thick, Keplerian accretion disc (Shakura & Sunyaev 1973) explains the continuous emission observed from a variety of accretion sources both on the galactic and extragalactic scale. This so-called SSD model cannot explain, however, the hard spectral component observed in most compact accreting systems, e.g., in X-ray binaries or in AGNs. This component is generally fitted by a power-law and ascribed to a thin, radiatively inefficient plasma. Such conditions are met in an advection dominated accretion flow (ADAF, for a review see Narayan et al. 1998) where the heat, that is liberated locally by a disc annulus due to viscous dissipation, is advected with the radial bulk flow and eventually disappears into the central black hole, without being emitted towards the observer.

The presence of reflection components and the Fe K$`\alpha `$ line observed in X-ray spectra of BHXB’s (Ebisawa et al. 1996) strongly suggest that cold gas from a disc and hot gas must coexist in a small volume. A possible configuration is the radial transition from an outer SSD to an inner ADAF. A bimodal disc model along this line, with the transition radius as a free parameter and with the ADAF and SSD parts treated as separate entities, allows to fit the spectra of persistent X-ray sources (Narayan et al. 1996; Manmoto et al. 1997; Esin et al. 2001). Furthermore, this model does also explain the different spectral states observed in some soft X-ray transients by varying the location of the transition radius (Esin et al. 1998). A number of these sources do also show high-frequency quasi-periodic oscillations (QPO) in their fourier-power spectra (Nowak & Lehr 1998). These examples demonstrate the necessity for a time-dependent treatment.

In a consistent theory, the disc transition radius is obtained naturally as a function of flow parameters and boundary conditions. For the stationary limit, this issue was addressed by Manmoto & Kato (2000). If the transition region between the inner ADAF and the outer SSD is sufficiently narrow (i.e., $`\delta R/R1`$), specific connection conditions can be derived which in turn significantly constrain the parameter space available for global solutions. The most restrictive condition is that the ADAF must support a minimum outward energy flux into the SSD, which is then liberated as radiation. The nature of this energy flux is still a matter of debate, in particular within the one-dimensional approximation (Abramowicz et al. 2000). Honma (1996) and Manmoto & Kato (2000) specifically considered energy transport by turbulent convection, and conclude that the transition radius for stationary flows is then constrained to a radial domain with boundaries depending on the model parameters.

## 2 The disc model

We investigate time-dependent transition models by solving the full Navier-Stokes equations for a vertically integrated axisymmetric disc. In the energy equation for the gas we include dissipation of a turbulent convective energy flux, and in addition we solve a radiative transfer equation in the diffusion approximation. Our disc model thus generalises the stationary transition model of Manmoto & Kato (2000) to the non-stationary case.

Assuming vertical hydrostatic equilibrium, the vertical pressure scale height is $`H=c_\mathrm{s}/\mathrm{\Omega }_\mathrm{K}`$, with $`c_\mathrm{s}`$ the isothermal sound speed and $`\mathrm{\Omega }_\mathrm{K}`$ the Keplerian angular velocity associated with the pseudo-Newtonian gravitational potential for the space-time of a non-rotating black hole with mass $`M`$ (Paczynsky & Wiita 1980). As a result of the vertical averaging procedure (e.g., Narayan & Yi 1995a), the thermodynamic variables which appear in the model equations are all vertically integrated quantities, such as surface density $`\mathrm{\Sigma }=2H\rho `$, total pressure $`P`$, internal energy density $`E`$ and radiation energy density $`U`$.

The continuity equation and the momentum equations in radial and azimuthal direction are

$$_t\mathrm{\Sigma }+\mathbf{}\mathrm{\Sigma }𝒗=0$$
(1)
$$_t\mathrm{\Sigma }v_R+\mathbf{}\mathrm{\Sigma }v_R𝒗=_RP+\mathrm{\Sigma }(\mathrm{\Omega }^2\mathrm{\Omega }_\mathrm{K}^2)R+(\mathbf{}\tau )_R$$
(2)
$$_t\mathrm{\Sigma }l+\mathbf{}\mathrm{\Sigma }l𝒗=(\mathbf{}\tau )_\varphi R,$$
(3)
where $`R`$ is the cylindrical radius and $`l`$ is the specific angular momentum. Viscous transport is modeled as a Navier-Stokes stress, i.e., the vertically integrated stress tensor $`\tau =2\mu \sigma `$, with $`\sigma `$ the shear tensor and $`\mu `$ the vertically integrated dynamical viscosity. We adopt the $`\alpha `$-viscosity law of Shakura & Sunyaev (1973) as a model for anomalous turbulent viscosity, $`\nu =\alpha c_\mathrm{s}H`$, which results in $`\mu =\mathrm{\Sigma }\nu =\mathrm{\Sigma }\alpha c_\mathrm{s}H`$.

The thermal disc structure is governed by the magnitude of the various heating and cooling terms which contribute to the energy equation for both the gas and radiation field. The gas energy equation is written as

$$_tE+\mathbf{}E𝒗+P_\mathrm{g}\mathbf{}𝒗=Q_{\mathrm{vsc}}^++Q_{\mathrm{trb}}^+Q_{\mathrm{rad}}^{}.$$
(4)
The vertically integrated viscous heating rate is given by $`Q_{\mathrm{vsc}}^+=(\tau \mathbf{})𝒗`$, with $`\tau `$ as given above.

As the entropy gradient in ADAFs decreases outward, ADAFs are convectively unstable (Honma 1996). Thus we additionally consider a bulk energy transport due to a convective energy flux $`𝑭_{\mathrm{trb}}`$ with the associated net heating rate given by $`Q_{\mathrm{trb}}^+=\mathbf{}𝑭_{\mathrm{trb}}`$. A natural assumption is to set $`F_{\mathrm{trb}}`$ proportional to the local entropy gradient, i.e.,

$$𝑭_{\mathrm{trb}}=\kappa _\mathrm{T}\mathrm{\Sigma }T\mathbf{}s,$$
(5)
with $`\kappa _\mathrm{T}`$ an effective diffusion coefficient. Following (Honma 1996; Manmoto et al. 2000; Manmoto & Kato 2000) we parametrise $`\kappa _\mathrm{T}`$ in analogy with the $`\alpha `$-viscosity prescription as $`\kappa _\mathrm{T}=\alpha _\mathrm{T}c_\mathrm{s}H`$. In general, $`\alpha _\mathrm{T}\alpha `$ is expected, since convective turbulence contributes to the viscosity, but not all sources of viscosity produce a bulk energy transport (Ryu & Goodman 1992). The rate at which heat is advected with the flow is $`Q_{\mathrm{adv}}^{}=\mathrm{\Sigma }T(𝒗\mathbf{})s=\mathbf{}E𝒗+P_\mathrm{g}\mathbf{}𝒗`$. The radiative cooling rate is $`Q_{\mathrm{rad}}^{}=2c\rho H\kappa _P(a_RT^4U/2H)`$, with $`\kappa _P`$ the Planck-mean opacity for free-free absorption.

Finally, the radiation energy equation is

$$_tU+Q_{\mathrm{diff}}^{}=Q_{\mathrm{src}}^+Q_{\mathrm{esc}}^{},$$
(6)
where $`Q_{\mathrm{src}}^+=Q_{\mathrm{rad}}^{}`$. Sinks for the radiative energy density originate from radiative losses through the upper and lower disk surfaces $`Q_{\mathrm{esc}}^{}`$ and from radial diffusion $`Q_{\mathrm{diff}}^{}=\mathbf{}𝑭_{\mathrm{rad}}`$. Following Manmoto & Kato (2000), we write

$$𝑭_{\mathrm{rad}}=\frac{cH}{3\rho H\kappa _R}_RU\widehat{𝒆}_𝑹,$$
(7)
where $`\kappa _R`$ is the Rosseland-mean opacity including electron scattering and free-free absorption. The equations are closed by an equation of state for the pressure, $`P=P_\mathrm{g}+P_{\mathrm{rad}}`$. We apply $`P_\mathrm{g}=(\gamma 1)E`$, with $`\gamma `$ the ratio of specific heats, together with $`P_{\mathrm{rad}}=U/3`$.

We discretise the equations on a non-uniform radial mesh, and use a time semi-implicit, operator-split scheme to advance the solution forward in time. The most restrictive constraint on the time-step in explicit schemes is due to radiative and viscous terms as a result of the associated time-scales. We avoid this problem by treating these terms implicitly. Implicit methods for finite volume discretisation of our equations result in a large sparse matrix, where only the neighbourhood of the diagonal is populated. The coupling between individual equations gives raise to a few off-diagonal entries. If this coupling is weak, e.g. if the time increment is small, these off-diagonal terms may be neglected at low order and the resulting matrix can be inverted very efficiently. In fact, we first solve each equation independently, then update the dynamical variables and finally recompute the global defects with the full system. The procedure is now repeated until global convergence is reached. This scheme was originally proposed by Hujeirat & Rannacher (1998) and is naturally advantageous for moderately coupled systems. If the coupling is too strong, the efficiency in terms of achievable physical simulation time over computational CPU-time is very low, and may even drop below that of an explicit scheme.

In this letter we restrict attention to the model parameters $`(\dot{m}_i,\alpha ,\alpha _\mathrm{T})`$, where $`\dot{m}_i`$ is the mass accretion rate of the initial model in Eddington units. The black hole mass and the adiabatic index are held fixed at $`M=10M_{}`$ and $`\gamma =1.5`$ in all models. As initial data we take a self-similar ADAF profile of Narayan & Yi (1994) with the advection efficiency parameter $`f`$ a step-function, i.e., with $`f`$ dropping from unity (initial ADAF) to zero (initial SSD) at an arbitrarily chosen initial transition radius. At the outer boundary, located at $`2000\mathrm{R}_\mathrm{g}`$ – where $`\mathrm{R}_\mathrm{g}=GM/c^2`$ is the gravitational radius – we fix the accretion rate to $`\dot{m}_i`$ and the specific angular momentum to the local Keplerian value. Vanishing gradient boundary conditions are applied to the remaining variables. At the inner boundary at $`2.2\mathrm{R}_\mathrm{g}`$, we apply free outflow boundary conditions.

## 3 A fiducial solution

The transition radius $`R_{\mathrm{tr}}`$ is defined as the location where the vertical optical depth is equal to unity. Variables at the transition radius are denoted by the subscript ’tr’, and the subscript ’K’ refers to Keplerian motion. Time is measured in units of the Keplerian orbital period taken at the initial transition radius, $`P_\mathrm{K}^{\mathrm{tr}}`$, unless explicitly noted otherwise. Our unit of length is the gravitational radius $`\mathrm{R}_\mathrm{g}`$. In Fig. 1 we plot the contributions to the energy balance for a fiducial model with parameters $`(\dot{m}_i,\alpha ,\alpha _\mathrm{T})=(0.021,0.4,0.4)`$. The transition in the initial conditions was located at $`R_{\mathrm{tr}}=55\mathrm{R}_\mathrm{g}`$. These parameters were chosen to facilitate comparison with other published results.

We generally identify three domains within the flow according to the magnitude of the different terms in the stationary energy equation (4), even though our solutions are intrinsically time-dependent and cannot be directly compared with steady state profiles. With increasing radial distance from the central black hole we identify the inner ADAF where the optical depth is very low. Next we define the transition region as the zone where the rotation is super-Keplerian. Super-Keplerian rotation is present whenever the transition region is narrow (cf. Abramowicz et al. 1998). Between the outer edge of the transition region and the outer edge of the disc we observe the SSD where the optical depth is very large. The radial profiles are never stationary, as discussed below. Nevertheless, the time-averaged profiles are rather similar to the stationary solutions of Manmoto & Kato (2000).

## 4 Evolution of the transition radius

The transition radius $`R_{\mathrm{tr}}`$ derived from our solutions is a function of time. We find that $`R_{\mathrm{tr}}`$ is always located within the transition region. After a short initial relaxation phase, we generally observe a quasi-periodic time variation of $`R_{\mathrm{tr}}`$, superimposed on a much slower outward directed drift. This behaviour is shown in Fig. 2(a).

The mass flow rate through the transition region is not constant in time, but modulated with the same quasi-periodic pattern as $`R_{\mathrm{tr}}`$. The transition region tends to abruptly deposit cold gas into the ADAF. Whenever this happens, the transition region grows at the expense of the ADAF.

On top of these rapid oscillations we observe a slow outward drift of the transition radius and consequently, the SSD shrinks at this very rate while the ADAF grows. For the fiducial model introduced above, $`R_{\mathrm{tr}}`$ grows from initially $`55\mathrm{R}_\mathrm{g}`$ to finally $`123\mathrm{R}_\mathrm{g}`$ at $`t=92P_\mathrm{K}^{\mathrm{tr}}`$. The associated drift velocity decreases with increasing $`R_{\mathrm{tr}}`$ or, accordingly, with increasing time. A drifting behaviour of this type is perceived in all our simulations. However, the overall bimodal disc structure is never disrupted and the transition region is always narrow in the sense defined above.

At $`t=92P_\mathrm{K}^{\mathrm{tr}}`$ the time-step of our simulation dropped drastically to very small values, and the simulation was halted. This behaviour was seen in all our simulations (Gracia 2002). Interestingly, the values we measure for the transition radius at that time is nearly equal to the maximum transition radius predicted for a stationary bimodal disc model with the same parameters and boundary conditions (Manmoto & Kato 2000). Near this maximum transition radius the equations become very stiff, and thus the convergence of our numerical method decreases while the time-step freezes. The issue of the drift as such is further discussed in section 4.2.

### 4.1 The oscillation cycle

The period of an oscillation cycle $`\mathrm{\Delta }t_{\mathrm{osc}}`$ increases with time. In Fig. 2(a) we plot the innermost location and the outermost location which the transition radius takes during the complete evolution of our fiducial model. The outer location is treated as the generic transition radius and fitted by an analytical function $`_{\mathrm{tr}}(t)`$. The innermost location was in turn considered as the perturbed transition radius in the course of the oscillation cycle and used to derive the oscillation period $`\mathrm{\Delta }t_{\mathrm{osc}}`$.

As shown in Fig. 2(b), a strong correlation is found between $`\mathrm{\Delta }t_{\mathrm{osc}}`$ and the generic transition radius $`_{\mathrm{tr}}(t)`$. Furthermore, $`\mathrm{\Delta }t_{\mathrm{osc}}`$ is very close to (within a factor of order unity) the Keplerian orbital period at the generic transition radius, $`P_\mathrm{K}(_{\mathrm{tr}}(t))`$, independent of model parameters $`(\dot{m}_i,\alpha ,\alpha _\mathrm{T})`$. However, the width of the transition region is a function of the adiabatic index $`\gamma `$ (Manmoto & Kato 2000) and therefore we expect it to modify the proportionality factor between $`\mathrm{\Delta }t_{\mathrm{osc}}`$ and $`P_\mathrm{K}(_{\mathrm{tr}}(t))`$.

A necessary condition for a narrow transition region is super-Keplerian rotation within the transition region, i.e. $`l>l_\mathrm{K}`$ (Abramowicz et al. 1998). In our calculations the specific angular momentum profile decreases outward and smoothly matches the Keplerian rotation profile of the SSD. However, if $`l`$ decreases with increasing radius, $`_Rl<0`$, the flow is Rayleigh unstable. Consequently, a narrow ADAF-SSD transition region is the location of putative instability.

Outward decreasing specific angular momentum profiles have previously been obtained by e.g. Honma (1996) and Manmoto & Kato (2000) as seen from their figures. The specific angular momentum profile of Abramowicz et al. (1998) does not decrease outward. They used an adiabatic equation of state, which yields broader transition regions, i.e. $`l`$ changes slower near the outer edge of the transition region and matches the SSD profile without actually decreasing.

Local stability analysis of a rotating flow predicts linear perturbations to oscillate with frequency $`\omega `$ that is bounded by the epicyclic frequency $`\kappa `$. Perturbations with $`\omega ^2\kappa _{\mathrm{max}}^2`$ will preferably grow or decay (Kato & Manmoto 2000; Kato 2001, for a review). If $`\kappa ^2`$ is negative the flow is dynamically unstable. If $`\kappa ^2`$ is positive the perturbation will oscillate with frequency $`\omega ^2`$ near $`\kappa _{\mathrm{max}}^2`$. Kato & Manmoto (2000) further argue that the oscillations will be trapped within the transition region, since they are reflected at the transition region boundaries.

In our solutions we find an upper limit $`\kappa _{\mathrm{max}}^2/\mathrm{\Omega }_\mathrm{K}^22`$, i.e., the oscillation period $`\mathrm{\Delta }t_{\mathrm{osc}}/P_\mathrm{K}^{\mathrm{tr}}`$ has a lower limit $`1/\sqrt{2}`$ of order unity, consistent with our numerical results and with our proposition that the oscillations have a physical nature. On the other hand, we note that $`\kappa _{\mathrm{max}}^2`$ sensitively depends on the angular momentum gradient in the transition region, which is rather poorly resolved in our current models. Therefore, we cannot conclude from our simulations that $`\mathrm{\Delta }t_{\mathrm{osc}}`$ is always bounded by the above limit.

We mentioned earlier that perturbations generated within the transition region are expected to be effectively trapped within the same. In reality, this simple picture does not survive the confrontation with numerical experiments entirely unharmed. Perturbations from within the transition region do back-react with the flow outside. Fig. 3 illustrates the dynamical evolution of the flow during an oscillation cycle at high temporal resolution.

The oscillation cycle appears to be initiated by a perturbation as seen, e.g., in the radial momentum flux density, approaching the inner boundary of the transition region. The perturbation leaks into the ADAF, cold matter with high angular momentum detaches from the transition region and penetrates the ADAF as seen in Fig. 3. The ejected blob carries nearly constant angular momentum, such that the region of super-Keplerian rotation, i.e., the transition region, grows far into the ADAF. When the perturbation reaches the sonic radius, it is partly reflected as an outgoing wave and sweeps up hot plasma towards the transition. In continuation the ADAF relaxes slowly back into the unperturbed state. This causes high variability not only of the transition region, but also of the inner ADAF, where most of the X-ray flux is expected to originate. Hence the evolution of our models is strongly affected by the existence of oscillations in the transition region.

### 4.2 Radial drift or the strong ADAF principle

The outward drift of the transition region which we see in our simulations does also require interpretation. It is well known that, for a given accretion rate, the radial extent of the ADAF (if existent) is restricted by an upper limit $`R_\tau `$ (e.g. Chen et al. 1995; Abramowicz et al. 1995). At every radius $`R>R_\tau `$, the optical depth exceeds unity and radiative cooling is efficient, such that the ADAF solution branch is no longer available. For certain $`R<R_\tau `$, the ADAF and SSD-like solutions are both available in principle. From entropy considerations we expect, that nature realizes the ADAF whenever possible, a proposition referred to as the strong ADAF principle (Narayan & Yi 1995b).

Our models do indeed seem to support this hypothesis. The transition radius $`_{\mathrm{tr}}(t)`$ does always drift outward, thereby increasing the total entropy of the entire solution. The drift is stopped once $`R_{\mathrm{tr}}`$ has reached the outermost location, $`R_{\mathrm{trb}}`$, that is allowed from stationary analysis. Note that $`R_{\mathrm{trb}}`$ is smaller than $`R_\tau `$ because turbulent convective energy transport is taken into account in our models.

For a given set of parameters, the specific entropy of the ADAF is, at any radius of interest, larger than the specific entropy of the SSD. The particular form applied for turbulent convective energy transport (namely along the negative entropy gradient), implies that the transition region is always supplied with energy from the inner ADAF. As long as this energy surplus cannot be completely shed by radiative losses, the transition region is heated to nearly the local virial temperature and the flow switches to the ADAF mode. As a consequence, the SSD is slowly evaporated radially at its inner edge. This process does not stop until radiative energy shedding in the transition region completely dominates energy supply by convection – exactly the definition of the outermost allowed transition radius $`R_{\mathrm{trb}}`$. Our models do indeed evolve in general within the parameter space of Manmoto & Kato (2000, see their Figs. 8, 9 and 10).

## 5 Concluding remarks

We have modelled the time evolution of bimodal accretion discs where an outer cooling dominated disc matches an inner ADAF-type flow through a narrow transition region. Our models show an intrinsically time-dependent behaviour that can be understood from the dynamical instabilities predicted for stationary ADAF-SSD transitions. Specifically, oscillations with frequency close to $`\mathrm{\Omega }_\mathrm{K}`$ are excited in the transition region. We associate these oscillations with a mutual local instability that originates from the radially decreasing specific angular momentum profile, a characteristic feature of any narrow-transition bimodal disc.

We observe that the waves excited in the transition region leak into the inner ADAF, but do not essentially propagate outward into the SSD. These oscillations give rise to rapid periodic modulations of the ADAF and the transition region. We estimated the total disk luminosity by integrating the radiative cooling rate over the disk surface and again find a periodic modulation with a variability of a few percent. This could possibly explain the QPOs observed in X-Ray transient systems (e.g. Nowak & Lehr 1998, and references therein). Kato & Manmoto (2000) estimated the oscillation frequency in bimodal disc transitions to a small fraction $`10^3`$ of $`\mathrm{\Omega }_\mathrm{K}`$, albeit under the assumption of small frequency and using a rather crude estimate for the disc structure. Our numerical results show frequencies very close to $`\mathrm{\Omega }_\mathrm{K}`$ at the transition radius.

We furthermore observe that the initial transition radius does slowly drift to larger radii, i.e., even though a stationary transition is allowed for a wide range of radii, the disc appears to favour the outermost possible transition radius. This fits well into the strong ADAF principle hypothesis (Narayan & Yi 1995b). Within our modelling, convective energy transport is the key to make this drift possible, i.e., to evaporate the SSD radially until the transition reaches its outermost allowed location. But even if we neglect convective energy transfer as modelled in this paper, the large entropy gradient across the transition still persists. Significant transfer of heat from the ADAF to the SSD is a necessary condition for the existence of a stationary transition front. Whether turbulence is the relevant physical mechanism for this heat flux and whether the shape of the transition front is actually close to cylindrical cannot be addressed by 1-dimensional modeling. The resolution of this issue does require at least 2-dimensional simulations in the parameter range that allows for a transition in the 1-dimensional stationary case. This region of the parameter space, however, has not been addressed in a 2-dimensional simulation by Hujeirat & Camenzind (2000). The gap between vertical disc evaporation models (Meyer & Meyer-Hofmeister 1994; Meyer et al. 2000) and the present 1-dimensional radial transition physics therefore remains to be closed in a future work.

## Acknowledgements

Part of this work was supported by the *Deutsche Forschungsgemeinschaft, DFG* (Sonderforschungsbereich 439 and Sonderforschungsbereich 359).