format: text

text:
# COSMOS: A Radiation-Chemo-Hydrodynamics Code for Astrophysical Problems

## 1 Introduction

This is the second in a series of papers using a new code (Cosmos) that we have developed for a broad range of astrophysical problems, including scalar- and radiation-field dominated processes in the early universe, cosmological structure formation, black-hole accretion, neutron-star binaries, astrophysical jets, and multiphase galactic dynamics. Such an ambitious range of topics requires that the code be able to accurately model a wide variety of physical processes, in both the relativistic and Newtonian regimes. Our first contributions therefore present the results of tests designed to examine the abilities of Cosmos to compute physical processes relevant for astrophysical problems.

Cosmos is a collection of massively parallel, multi-dimensional, multi-physics solvers utilizing the MPI paradigm for parallel computing of both Newtonian and general relativistic astrophysical problems. It currently includes several different options for computational fluid dynamics (CFD) methods, equilibrium and non-equilibrium primordial chemistry with 27 atomic and molecular reactions, various radiative cooling processes, nonequilibrium radiation flux-limited diffusion, radiation pressure, relativistic scalar fields, Newtonian external and self-gravity, arbitrary spacetime curvature in the form of a generic background metric, and viscous stress in a fully covariant formulation. In general, Cosmos assumes no particular symmetry in the equations, and is therefore designed to run on structured Cartesian meshes. However, due to the covariant formulation adopted for the hydrodynamics equations, Cosmos can also be run on various grid geometries (e.g., Cartesian, cylindrical, spherical) for problems using the relativistic fluid dynamics solvers.

Numerical methods used to solve the hydrodynamics equations include a total variation diminishing (TVD) Godunov solver for Newtonian flows using Roe’s (Roe, 1981) approximate Riemann solver and a third order Runge-Kutta time-marching scheme (Shu & Osher, 1988). For either Newtonian or general relativistic systems, two artificial viscosity (vonNeumann & Richtmyer, 1950) algorithms are available: a non-staggered grid method in which all variables are located at the zone-center, and a staggered grid and field centering method similar to the Zeus code (Stone & Norman, 1992) in which scalar (vector) quantities are located at the zone (face) centers. Two additional options for CFD methods are included for both Newtonian and relativistic problems that are based on non-oscillatory central difference schemes (Jiang et al., 1998), also differentiated by grid centering: staggered versus centered in time as well as space. The relativistic hydrodynamics algorithms have been presented in Anninos & Fragile (2003) along with several tests of the code, including relativistic shock tube, wall shock, and dust accretion problems. Here we emphasize the Newtonian and multi-physics descriptions and code tests of primordial chemistry, radiative cooling, and radiation diffusion coupled together with hydrodynamics.

We proceed in §2 by describing the basic equations used by Cosmos, and in §3 by presenting tests of the code. The tests are designed to examine the ability of the code to follow shocks, blast waves and dynamical instabilities, to advect materials, to resolve heating and cooling flows, to transfer radiation, to simulate chemical networks, and to calculate self-gravitating gas distributions.

## 2 Basic Equations

The Newtonian multi-species mass, momentum and energy (hydrodynamic and radiation diffusion) continuity equations are written in an Eulerian frame as

$`{\displaystyle \frac{\rho }{t}}+_i(\rho v^i)`$ $`=`$ $`0,`$ (2-1)
$`{\displaystyle \frac{\rho ^{[m]}}{t}}+_i(\rho ^{[m]}v^i)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N_s}{}}}{\displaystyle \underset{j=1}{\overset{N_s}{}}}k_{ij}(T)\rho ^{[i]}\rho ^{[j]}+{\displaystyle \underset{i=1}{\overset{N_s}{}}}I_i(\nu )\rho ^{[i]},`$ (2-2)
$`{\displaystyle \frac{(\rho v^k)}{t}}+_i(\rho v^kv^i)`$ $`=`$ $`_k(P+P_R)\rho _k\varphi ,`$ (2-3)
$`{\displaystyle \frac{e}{t}}+_i(ev^i)`$ $`=`$ $`P_iv^i+c\rho (\sigma _aE\sigma _pa_rT^4)+\mathrm{\Lambda }(T,\rho ^{[m]}),`$ (2-4)
$`{\displaystyle \frac{E}{t}}+_i(Ev^i)`$ $`=`$ $`_i\left({\displaystyle \frac{c}{3\rho \sigma _r}}^iE\right){\displaystyle \frac{E}{3}}_iv^ic\rho \sigma _aE+c\rho \sigma _pa_rT^4,`$ (2-5)
where $`v^k`$ is the fluid velocity assumed to be the same for each of the chemical species, $`e`$ is the fluid energy density, $`\rho ^{[m]}`$ are the species densities satisfying $`\rho =_m\rho ^{[m]}`$ for the total density, $`\varphi `$ is the gravitational potential obtained from Poisson’s equation $`^2\varphi =4\pi G\rho `$, $`P`$ is the fluid pressure, $`P_R=E/3`$ is the radiation pressure, $`a_r(=4a/c)`$ is the radiation constant, $`a`$ is Stephan-Boltzmann’s constant, $`c`$ is the speed of light, $`E`$ is the radiation energy, $`\sigma _r`$, $`\sigma _a`$, and $`\sigma _p`$ are the Rosseland, absorption and Planck mean opacities, $`k_{ij}(T)`$ are the rate coefficients for the 2-body reactions which are functions of fluid temperature $`T`$, and $`I_j`$ are frequency-integrated photoionization and dissociation rates. The summations in (2-2) are over the $`N_s`$ atomic and molecular species included in the chemistry model, up to nine: H I, H II, He I, He II, He III, $`e^{}`$, H-, H2, H$`{}_{}{}^{+}{}_{2}{}^{}`$. A total of 27 chemical reactions are included in the full network, which we summarize in Table 1 for convenience, but refer the reader to Abel et al. (1997) and Anninos et al. (1997) for more complete descriptions of the reactions, rate coefficients, and numerical methods.

Assuming that reasonable models are provided for the mean opacities of the fluid, radiation energy is coupled to the net fluid momentum and internal energy, accounting for non-equilibrium heating and cooling effects in the single temperature and gray (spectral average) approximations. However, we have not currently coupled radiation in a self-consistent fashion to the chemistry solver and ionization states of the fluid composition. This is a much more difficult task due to grid resolution requirements, computational time constraints, multi-species interactions, multi-group transport, and non-LTE effects expected to be important in some of our applications. We will address these problems in future developments of Cosmos.

Three radiative cooling and heating models are implemented in the optically thin limit for $`\mathrm{\Lambda }(T,\rho ^{[m]})`$ in (2-4), depending upon whether the chemistry is known. First, if chemistry is not solved, the cooling function is set proportional to the square of the total number density with an empirical model for the ionization fraction that approximates linearly the equilibrium result (2-8) described below. In particular,

$$\mathrm{\Lambda }(T,\rho )=\left[\underset{i}{}\dot{e}_i(T)(f_I\rho )^2+\dot{e}_M(T)(f_M\rho )^2+\mathrm{\Lambda }_C\right]\times \{\begin{array}{cc}e^{(TT_{min})/\delta T}\hfill & \text{if}TT_{min},\hfill \\ 1\hfill & \text{otherwise},\hfill \end{array}$$
(2-6)
where $`\dot{e}_M(T)`$ is the temperature-dependent cooling rate from metals (including carbon, oxygen, neon, and iron lines taken from Bohringer and Hensler (1989), with a metallicity dial to scale relative to solar abundance), $`\dot{e}_i(T)`$ is the cooling rate from hydrogen and helium lines, $`f_M`$ is the mass fraction of metals, $`f_I`$ is an estimate of the ionization fraction defined as $`\text{min}(1,\text{max}(0,(T_{eV}T_c)/3))`$ with $`T_c=1eV`$ to match roughly the upper and lower bounds in the equilibrium model described below. Also,

$$T=\frac{m_p(\mathrm{\Gamma }1)}{k_B(1+f_I)}\frac{\mu e}{\rho }$$
(2-7)
is the gas temperature in Kelvin, $`\mu `$ is the mean molecular weight, $`m_p`$ is the proton mass, $`k_B`$ is Boltzmann’s constant, and the exponential is introduced to suppress cooling at low temperature ($`T_{min}`$ is typically set to $`10^4K`$, with width $`\delta T500K`$). $`\mathrm{\Lambda }_C`$ represents Compton cooling (or heating) due to the interaction of free electrons with the cosmic microwave background.

The second model applies when the chemistry is solved in equilibrium and the ionization fraction is determined from equating the dominant hydrogen recombination and collisional rates (Bond et al., 1984; Anninos et al., 1994)

$$\frac{f}{1f}=3.2\times 10^4T_{eV}^{1.22}e^{I_\mathrm{H}/T_{eV}},$$
(2-8)
where $`T_{eV}`$ is the temperature in electron-volts, $`I_\mathrm{H}=13.6eV`$ is the ionization energy of hydrogen. Assuming the electrons and ions have the same temperature, the gas pressure and energy can be written

$`p`$ $`=`$ $`p_I+p_e=(1+f){\displaystyle \frac{\rho k_BT}{\mu m_p}},`$ (2-9)
$`e`$ $`=`$ $`{\displaystyle \frac{1+f}{\mathrm{\Gamma }1}}{\displaystyle \frac{\rho k_BT}{\mu m_p}}+{\displaystyle \frac{\rho fk_B\stackrel{~}{k}I_\mathrm{H}}{\mu m_p}},`$ (2-10)
with equation of state

$$p=\left(\mathrm{\Gamma }1\right)\left(e\frac{\rho fk_B\stackrel{~}{k}I_\mathrm{H}}{\mu m_p}\right),$$
(2-11)
where $`\stackrel{~}{k}(=11605K/eV)`$ is a numerical constant converting between electron-volts and Kelvin. The cooling function takes the same form as (2-6), but in this case the ionization fraction $`f_I`$ is computed iteratively from the equilibrium model $`f_I=f`$ of equation (2-8), and the exponential cutoff is not applied since this behavior is contained implicitly in $`f`$.

Equilibrium chemistry is solved by setting the time derivatives in equations (2-2) to zero (but only for the local chemistry update, transport is still accounted for with operator splitting) and solving the following algebraic equations, after first making a guess for the electron number density

$`n_{\mathrm{HI}}`$ $`=`$ $`{\displaystyle \frac{k_2n_\mathrm{H}}{k_1+k_2+k_{20}/n_e}},`$ (2-12)
$`n_{\mathrm{HII}}`$ $`=`$ $`n_\mathrm{H}n_{\mathrm{HI}},`$ (2-13)
$`n_{\mathrm{HeII}}`$ $`=`$ $`{\displaystyle \frac{yn_\mathrm{H}}{1+k_4/(k_3+k_{21}/n_e)+k_5/k_6+k_{22}/(k_6n_e)}},`$ (2-14)
$`n_{\mathrm{HeI}}`$ $`=`$ $`{\displaystyle \frac{k_4n_{\mathrm{HeII}}}{k_3+k_{21}/n_e}},`$ (2-15)
$`n_{\mathrm{HeIII}}`$ $`=`$ $`n_{\mathrm{HeII}}\left({\displaystyle \frac{k_5}{k_6}}+{\displaystyle \frac{k_{22}}{k_6n_e}}\right),`$ (2-16)
$`n_e`$ $`=`$ $`n_{\mathrm{HII}}+n_{\mathrm{HeII}}+2n_{\mathrm{HeIII}},`$ (2-17)
where $`f_\mathrm{H}`$ is the mass fraction of hydrogen, $`n_\mathrm{H}=f_\mathrm{H}\rho /m_p`$, and $`y=(1f_\mathrm{H})/(4f_\mathrm{H})`$ parameterizes the relative helium concentration. Equations (2-12) - (2-17) are solved iteratively in the order they are written until the electron density converges to a specified tolerance, typically less than $`10^{10}`$. The subscripts used for the rate coefficients $`k_i`$ refer to a particular chemical reaction ordered as in Table 1.

The third cooling model is used when chemistry is solved dynamically with the full nonequilibrium equations (2-2). In this case, the various ionization states and concentration densities $`\rho ^{[i]}`$ of each element are calculated from the chemistry equations with a stable ordered backwards differencing scheme (Anninos et al., 1997) and used explicitly in the cooling function as

$$\mathrm{\Lambda }(T,\rho ^{[m]})=\underset{i=1}{\overset{N_s}{}}\underset{j=1}{\overset{N_s}{}}\dot{e}_{ij}(T)\rho ^{[i]}\rho ^{[j]}+\underset{i=1}{\overset{N_s}{}}J_i(\nu )\rho ^{[i]}+\dot{e}_M(T)(f_M\rho )^2+\mathrm{\Lambda }_C,$$
(2-18)
where $`\dot{e}_{ij}(T)`$ are the cooling rates from 2-body interactions between species $`i`$ and $`j`$, and $`J_i`$ represents frequency-integrated photoionization and dissociation heating. The equation of state in this case is given by

$$T=\frac{(\mathrm{\Gamma }1)}{k_B}\frac{e}{_i^{N_s}n_i}.$$
(2-19)
We account for a total of eight different cooling and heating mechanisms: collisional-excitation, collisional-ionization, recombination, bremsstrahlung, metal-line cooling (dominantly carbon, oxygen, neon, and iron), molecular-hydrogen cooling, Compton cooling or heating, and photoionization heating.

Finally we note that the following conservation equations for the chemical concentrations of hydrogen, helium, and charge

$`\rho _\mathrm{H}+\rho _{\mathrm{H}^+}+\rho _\mathrm{H}^{}+\rho _{\mathrm{H}_2^+}+\rho _{\mathrm{H}_2}=\rho f_\mathrm{H},`$
$`\rho _{\mathrm{He}}+\rho _{\mathrm{He}^+}+\rho _{\mathrm{He}^{++}}=\rho (1f_\mathrm{H}),`$
$`\rho _{\mathrm{H}^+}\rho _\mathrm{H}^{}+{\displaystyle \frac{1}{2}}\rho _{\mathrm{H}_2^+}+{\displaystyle \frac{1}{4}}\rho _{\mathrm{He}^+}+{\displaystyle \frac{1}{2}}\rho _{\mathrm{He}^{++}}=m_pn_e,`$
are enforced after each computational cycle update.

## 3 Numerical Methods and Tests

Cosmos uses structured, block-ordered meshes for both spatial finite differencing and finite volume discretization methods. Depending on the hydrodynamic algorithm, state variables are either defined all at the zone centers (for one of the artificial viscosity schemes and for the total energy conserving methods - TVD Godunov and non-oscillatory central difference schemes), or on a staggered mesh for a second artificial viscosity method in which scalar and tensor quantities are zone-centered and vector variables are face-centered. Periodic, reflection, constant-in-time, user-specified, outgoing, and flat (vanishing first derivative) boundary conditions are implemented. Both the relativistic and Newtonian hydrodynamic equations are solved using single or multiple step time-explicit, operator-split methods with second-order spatial finite differencing. The gravitational potential and nonequilibrium radiation diffusion equations (discretized implicitly without the transport and compressive terms) are solved using a collection of linear matrix solvers from the Hypre software package developed at Lawrence Livermore National Laboratory (Falgout et al., 2002). Hypre includes options for different conjugate gradient and multigrid methods with various preconditioners to optimize performance. We use a multigrid solver for the gravity and diffusion tests in sections §3.5, §3.6, and §3.7.

Since the main emphasis in this paper is on Newtonian hydrodynamics, the code tests presented in the following sections are designed to verify our code only in that regime, along with multi-physics (e.g. radiation and chemistry) coupling. We refer the interested reader to (Anninos & Fragile, 2003) for discussions of the relativistic tests and for more explicit details of the numerical algorithms. It is anticipated that problems run in the Newtonian regime using Cosmos shall include microphysics (chemistry, heating, and cooling). Such problems are most appropriately solved with one of the artificial viscosity methods, which are written in an internal energy formulation. All of the hydrodynamics tests shown here are therefore computed using the staggered-mesh artificial viscosity method. We have, however, confirmed that results from the other algorithms yield comparable accuracy for most of the tests. The exceptions being the non-oscillatory central difference schemes which are more diffusive, in general, but particularly so for the Rayleigh-Taylor and spherical polytrope tests. Accuracy comparable to the artificial viscosity approach can be achieved for those tests if the grid resolution is roughly doubled.

### 3.1 Shock Tube

We begin testing with one of the standard problems in fluid dynamics, the shock tube or Sod problem. This test is characterized initially by two different fluid states separated by a membrane. At $`t=0`$ the membrane is removed and the fluid evolves in such a way that five distinct regions appear in the flow: an undisturbed region at each end, separated by a rarefaction wave, a contact discontinuity, and a shock wave. This problem provides a good test of the shock-capturing properties of the code since it has an exact solution (Sod, 1978) against which numerical results can be compared.

The initial state is specified by $`\rho _L=1`$, $`P_L=1`$, and $`V_L=0`$ to the left of the membrane, and $`\rho _R=0.125`$, $`P_R=0.1`$, and $`V_R=0`$ to the right. The fluid is assumed to be an ideal gas with $`\mathrm{\Gamma }=1.4`$, and the integration domain extends over a unit grid from $`x=0`$ to $`x=1`$, with the membrane located at $`x=0.5`$. The results presented here were run using scalar artificial viscosity with a quadratic viscosity coefficient $`k_{q2}=2.0`$, linear viscosity coefficient $`k_{q1}=0.3`$, and Courant factor $`k_c=0.6`$, in the notation of Anninos & Fragile (2003). Figure 1 shows spatial profiles of the results at time $`t=0.2`$ for the 64-zone 1D case and the $`64^3`$-zone 3D case (along the main diagonal). Table 2 summarizes the errors in $`\rho `$, $`P`$, and $`V`$ for different grid resolutions using the $`L`$-1 norm (i.e., $`E(a)_1=_{i,j,k}\mathrm{\Delta }x_i\mathrm{\Delta }y_j\mathrm{\Delta }z_k|a_{i,j,k}^nA_{i,j,k}^n|`$, where $`a_{i,j,k}^n`$ and $`A_{i,j,k}^n`$ are the numerical and exact solutions, respectively, and $`j=k=\mathrm{\Delta }y_j=\mathrm{\Delta }z_k=1`$ for 1D problems). The convergence rates for these tests are just under first order, as expected for shock-capturing methods. The slightly higher $`L`$-1 norm errors in the 3D case are due to the fact that the error calculation is computed globally across the whole mesh, and so suffers from reflection effects at the grid boundaries.

### 3.2 Sedov Blast Wave

The next problem we consider is the Sedov blast wave in which energy is released at $`t=0`$ in the form of an explosion into an initially undisturbed, uniform gas. In 3D, this results in a spherical shock wave (or blast wave) expanding from the explosion, such that $`r_St^{2/5}`$, where $`r_S`$ is the radius of the shock and $`t`$ is the elapsed time (Sedov, 1959). This problem encompasses a number of useful tests for our code, as it determines how well the code can follow a spherical shock wave as well as testing energy conservation.

The initial state is specified by $`\rho _0=1`$, $`P_0=3.33\times 10^{11}`$, and $`V_0=0`$. The fluid is assumed to be an ideal gas with $`\mathrm{\Gamma }=4/3`$, and the integration domain extends over a unit cube. The blast wave is initiated by significantly increasing the energy density (relative to the background) in an approximately spherical region with a half-width at half-maximum of two zones and maximum cutoff radius of five zones using a Gaussian profile. The initial energy density contrast between the peak of the Gaussian and the background is $`6.6\times 10^{11}`$. The results presented here were run using a tensor artificial viscosity with a quadratic viscosity coefficient $`k_{q2}=1.0`$, linear viscosity coefficient $`k_{q1}=0.3`$, and Courant factor $`k_c=0.4`$. Figure 2 shows the shock radius as a function of time fit with a curve of the form $`t^{2/5}`$. Figure 3 shows spatial profiles of the density along the $`x`$, $`y`$, and $`z`$ axes, as well as the diagonal at time $`t=2.1`$ for a $`64^3`$-zone octant, demonstrating that the blast wave maintains its self-similar solution along all three major axes and the diagonal. The energy loss in this problem was about 19%, although much of this was at the beginning of the simulation, and the total energy approaches a constant value after a time $`t=0.28`$.

### 3.3 Rayleigh-Taylor Instability

The growth of a classical Rayleigh-Taylor instability has been modeled in two-dimensions. These models serve as tests of the ability of the code to follow the growth of a classical instability in both the linear and non-linear phases, and of its ability to cleanly advect material across the grid during the non-linear growth.

The system which we model has physical dimensions that run from 0 to 0.1 in $`x`$, and from -0.35 to 0.15 in $`z`$. The resolution is 128 zones in $`x`$ and 1280 in $`z`$. The constant gravitational field, $`g=1`$, points toward negative $`z`$. At $`z=0`$, the density and isothermal sound speed of the dense (upper) fluid are $`\rho _u=1`$ and $`(c_s)_h^2=2.4`$, while those of the light (lower) fluid are $`\rho _l=0.1`$ and $`(c_s)_l^2=24`$, such that pressure is continuous across the fluid interface. Away from the interface, the individual fluids are isothermal, and their densities vary so as to maintain hydrostatic equilibrium in the gravitational field, i.e.

$$\rho (z)=\rho (0)\mathrm{exp}\left[\varphi (z)/c_s^2\right],$$
(3-1)
where $`\varphi (z)`$ is the gravitational potential, and we take $`\varphi (0)=0`$. To approximate incompressibility, the fluids are taken to be ideal gases, with large adiabatic index, $`\mathrm{\Gamma }=10`$.

The initial distribution of the fluids, in which quantities depend only upon $`z`$, are perturbed by introducing a vertical shift of the form

$$\delta z(x)=A\mathrm{cos}\left(2\pi x/\lambda \right),$$
(3-2)
where $`\lambda `$ is the wavelength of the perturbation, and $`A`$ is its amplitude. We have examined models having two different values of $`A/\lambda `$. In the first, we choose $`\lambda =0.1`$, and $`A=0.01`$, such that the perturbation is linear, though not strongly so. Resolution requirements and computational time constraints prevent us from modeling a system having a substantially smaller value of $`A`$. In the second model, $`A=\lambda =0.1`$, such that this model examines the growth of initially non-linear perturbations. In introducing the perturbations, density and pressure are not altered, and so the fluids remain in pressure equilibrium with each other (thus preventing the initial growth of sound waves), but are slightly out of equilibrium with the gravitational field. Reflective boundary conditions are applied at the top and bottom of the grid, while periodic boundaries are used at the left and right. Artificial viscosity is not used in these calculations.

Figure 4 shows snapshots of the evolution of the model having $`A/\lambda =0.1`$ at the times $`t=0,0.2,0.5,0.7,`$ and $`1.1`$. Shown, in grayscale, is the tracer material initially placed in the dense upper fluid. The initial perturbation can be seen in the first panel. As can be seen in the subsequent panels, the discretization of the grid creates short-wavelength structure superimposed upon the long-wavelength perturbation. The shortest wavelengths saturate quickly, but longer wavelengths persist, and can be seen superimposed upon the classical single-mode rollup at $`t=0.5`$. Although the small scale vorticity spreads the tracer material widely across the grid, little diffusiveness is apparent.

The early growth of the longest wavelength modes of both non-linear and linear initial perturbations are shown in Figure 5. Also shown is the prediction of linear theory for incompressible, constant density fluids, for which the instability is predicted to grow as $`e^{\omega t}`$, with

$$\omega ^2=\frac{2\pi g}{\lambda }\left(\frac{\rho _u\rho _l}{\rho _u+\rho _l}\right),$$
(3-3)
(Chandrasekhar, 1961). The overall agreement between theory and the model with the linear initial perturbation is good. The computed growth rate is found to be 14% slower than the analytic value. The slow growth rate may be due to the relatively large initial amplitude of the imposed perturbation, a surmise supported by the curvature seen in the growth rate of the numerical results. Both curves stand in contrast to that of the initially non-linear perturbation, which has a much smaller growth rate, as predicted by theory (see below).

In the non-linear regime, the amplitude of the Rayleigh-Taylor instability is known to behave asymptotically as

$$A(t)=\alpha _sg\frac{\rho _u\rho _l}{\rho _u+\rho _l}t^2,$$
(3-4)
(e.g. Youngs, 1994; Glimm et al., 2001). The late-time behavior of the models is shown in Figure 6. In the figure are shown the penetration amplitude of the dense fluid as a function of $`t^2`$ for both the initially linear and initially non-linear simulations. As can be seen in the figure, both simulations display the expected late-time behavior, with $`\alpha _s0.17`$ for the initially linear model, and $`\alpha _s0.04`$ for the initially non-linear model, where $`\alpha _s`$ is computed for the initial Atwood number of the fluids. The differences between the two models are most likely due to compressibility effects, which become more apparent at late times.

### 3.4 Radiative Shock Waves

The initial data for this test is characteristic of pre-shock flows expected from large-scale (galactic-type) structures

$`\rho `$ $`=`$ $`4.72\times 10^{25}\mathrm{g}\mathrm{cm}^3,`$
$`e`$ $`=`$ $`1.0\times 10^{30}\mathrm{g}\mathrm{cm}^1\mathrm{s}^2,`$
$`v^x`$ $`=`$ $`u_{\mathrm{in}}=1.7\times 10^7\mathrm{cm}\mathrm{s}^1,`$
corresponding to a uniform flow of gas along the $`x`$ direction. Reflection boundary conditions are imposed at $`x=0`$ and we use 100 zones to resolve a spatial extent of $`L=2.43\times 10^4`$ Mpc. A shock wave forms at the wall at $`x=0`$ and propagates to the right at velocity $`v_su_{\mathrm{in}}/3`$. As the heated gas cools radiatively, the shock begins to lose pressure support and slows down. Eventually the shock collapses and re-establishes a new pressure equilibrium closer to the wall. As gas continues to accrete, the shock front moves outwards again to repeat the cycle of oscillations as shown in Figure 7, where the shock position, $`x_s`$, is plotted as a function of time in units where the grid length is set to unity and the unit of time is $`10^{15}`$ seconds.

Figure 7 shows results from two calculations: one with no chemistry in which $`\mathrm{\Lambda }(T,\rho )\rho ^2T^{1/2}`$, and a second 6-species equilibrium chemistry model that approximates the number densities of the dominant coolants $`n_{\mathrm{HI}},n_{\mathrm{HII}},n_{\mathrm{HeI}},n_{\mathrm{HeII}},n_{\mathrm{HeIII}},n_e`$ from equilibrium assumptions. We have also run a third calculation to test the 6-species nonequilibrium chemistry model which explicitly solves the non-linear differential kinetics equations. The results in this case are nearly identical to the equilibrium calculations so we do not include them in Figure 7. We do, however, show the mass fraction distribution of each of the chemical species at the final time of the simulation for the nonequilibrium case in Figure 8. Mass fractions in the hot phase, where collisional ionization and recombination effects are expected to be in equilibrium, agree with those computed from the equilibrium model. For all these calculations we assume a perfect fluid model for the gas with adiabatic index $`\mathrm{\Gamma }=5/3`$ and a cooling function dominated by bremsstrahlung effects.

The accuracy of these calculations is evaluated by comparing the fundamental frequency of oscillations to the perturbation results of Chevalier & Imamura (1982) who defined the normalized frequency as $`\omega _I=(2\pi /P)(\overline{x}_s/u_{in})`$, where $`P`$ is the period of the oscillations, $`\overline{x}_s`$ is the average shock position, and $`u_{in}`$ is the inflow velocity. We find $`\omega _I=`$ 0.319 and 0.315 for the no-chemistry and 6-species chemistry cases respectively. These results compare nicely with the perturbation estimate of 0.31.

### 3.5 Marshak Wave

In the following test, we consider the penetration of radiation from a hot source into cold material (Marshak, 1957). Because radiative transport is very efficient in the problems considered here, significant penetration can occur on a timescale much shorter than the timescale for motion of the gas. We therefore ignore hydrodynamic transport and consider only supersonic radiation diffusion through an ambient gas in thermal equilibrium with the radiation and held initially at temperature $`T_0`$, but subject to the boundary condition $`T=T_1>T_0`$ at one end.

The energy diffusion equation

$$\frac{E}{t}=(DE),$$
(3-5)
where $`D=c/(3\rho \sigma _r)`$ is the diffusion coefficient, and $`\sigma _r=\kappa _0(\rho /\rho _0)^\gamma (T/T_0)^m`$ is the Rosseland mean opacity, can be solved approximately for this simple case to give (Long & Tahir, 1986)

$$E(x,t)E_1\left\{\xi _0^2\left[\frac{\overline{n}}{\overline{n}+1}\left(1\xi /\xi _0\right)^{(\overline{n}+1)/\overline{n}}\left(1\frac{1\xi /\xi _0}{2\overline{n}+1}\right)\right]\right\}^{1/\overline{n}},$$
(3-6)
where

$`\xi `$ $`=`$ $`\left({\displaystyle \frac{\overline{n}x^2}{2\eta ^2Dt}}\right)^{1/2},`$ (3-7)
$`\xi _0`$ $`=`$ $`\left({\displaystyle \frac{\overline{n}x_f^2}{2\eta ^2Dt}}\right)^{1/2}=\left({\displaystyle \frac{(\overline{n}+1)(\overline{n}+0.5)}{\overline{n}^2}}\right)^{1/2},`$ (3-8)
$`\overline{n}=(m+4)/4`$, and $`\eta `$ is a numerical fitting factor of order unity. From (3-8) we clearly see that the radiation front location, $`x_f`$, should propagate as $`t^{1/2}`$.

In this work we consider two cases: a constant opacity with $`\sigma _r=\kappa _0`$, and a temperature-dependent opacity with $`\sigma _r=\kappa _0(T/T_0)^m`$ and $`m=3`$. The computational grid is set to approximately 50 times the mean-free path ($`l=1/3\rho \sigma _r`$), and an appropriate stopping time $`t_{final}`$ is defined using equation (3-8). The problem is initialized with $`\kappa _0=1`$ cm2 g-1, $`\rho =1`$ g cm-3, $`T_0=10^4`$ K, $`T_1=10^6`$ K, and the grid is discretized with 400 zones. Figure 9 shows the profile of the radiation front for the constant opacity case normalized to the analytic radiation front position with $`\eta =1.1`$. The numerical results, plotted at (0.25, 0.5, 0.75, 1)$`\times t_{final}`$, nicely display the self-similar nature of this solution. Figure 10 plots the numerical radiation front location as a function of time. Here we define the numerical front to be located at the half-maximum, although our results do not depend strongly upon this choice. The data is well fit with a $`t^{1/2}`$ curve. Figures 11 and 12 present similar results for the temperature-dependent opacity. The analytic curve in Figure 11 is computed with $`\eta =1.04`$.

### 3.6 $`\mathrm{\Gamma }=2`$ Polytrope

Here we test the linear matrix methods used in solving Poisson’s equation for the gravitational potential of a compact self-gravitating source, and the ability of the code to maintain a balance between gravitational and pressure support forces in three dimensions. For this test we adopt an adiabatic polytrope star with $`\mathrm{\Gamma }=2`$ in spherical symmetry and hydrostatic equilibrium with solutions for the gas pressure, density and gravitational potential

$`P`$ $`=`$ $`2\pi G\alpha ^2\rho ^2,`$ (3-9)
$`\rho `$ $`=`$ $`{\displaystyle \frac{\alpha \rho _c}{r}}\mathrm{sin}\left({\displaystyle \frac{r}{\alpha }}\right),`$ (3-10)
$`\varphi `$ $`=`$ $`{\displaystyle \frac{GM}{R_s}}{\displaystyle \frac{4\pi G\rho _c\alpha ^3}{r}}\mathrm{sin}\left({\displaystyle \frac{r}{\alpha }}\right),`$ (3-11)
where $`\alpha =(M/(4\pi ^2\rho _c))^{1/3}`$ for radii $`rR_s`$, and $`R_s=\alpha \pi `$ is the outer surface radius. Also, $`\varphi =GM/r`$ with negligible density and pressure ($`<10^4`$ of the maximum central values) outside the star at radii $`r>R_s`$. These tests are carried out for characteristic neutron star parameters with total mass $`M=1.4M_{}`$, central density $`\rho _c=2.5\times 10^{15}`$ g cm-3, and radius 5.6 km.

We ran a sequence of three simulations at different grid resolutions over two sound crossing times, $`t=2T_s=2R_s/\sqrt{2\pi G\alpha ^2\rho _c}`$, with monopole boundary conditions to find mean relative errors of 0.0687, 0.0209, and 0.0101 for resolutions $`12^3`$, $`24^3`$, and $`48^3`$ cells respectively. Monopole boundary conditions are implemented by computing the total mass $`M_T`$ and mass centroid coordinates $`x_c^k=_{cells}\rho x^k\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z/M_T`$ in the computational domain. The potential field is then set to $`\varphi =GM_T/\sqrt{_k(x_o^kx_c^k)^2}`$ at the outer boundary $`x^k=x_o^k`$. Figure 13 shows spatial profiles in density for the $`48^3`$ case along the $`x`$-axis at the initial time (the analytic solution), and along four separate directions at $`t=2T_s`$: the $`x`$, $`y`$, $`z`$ and diagonal lines running through the origin. The data is displayed in dimensionless code units with length scale $`\stackrel{~}{L}=R_s/N`$ where $`N=48`$ is the number of interior cells along an axis, and density scale $`\stackrel{~}{D}=M_{}/\stackrel{~}{L}^3`$. We observe no significant break in symmetry in the solutions, and the central peak density is maintained to good accuracy.

### 3.7 $`\mathrm{\Gamma }=2`$ Polytrope with Radiation Pressure

The hydrostatic polytrope solution in §3.6 can be generalized to include a radiation field and radiation pressure, thus providing a useful test of the coupling between gravity, fluid pressure, and radiation pressure. It also exercises the multiphysics operator splitting scheme in testing the ability of the code to maintain a balance between the three different self-consistently generated forces in three dimensions.

The solution (3-9) - (3-11) is easily extended to include the effect of radiation pressure by assuming $`P_R=\beta P`$, where $`\beta `$ is a constant. This simplifies the solution considerably and allows for an effective (hydrodynamic plus radiation) pressure to counteract gravity with the same radial dependence as in the pure hydrodynamics case. The radiation energy equation (2-5) reduces, in the spherical, static, and thermal equilibrium limits, to $`r^2_r(r^2D_rE)=0`$, which can be solved trivially if the diffusion coefficient is set to $`D=k_0/(r^2_rE)`$, where $`k_0`$ is a constant. The complete solution including radiation pressure and radially dependent opacity is

$`P`$ $`=`$ $`{\displaystyle \frac{2\pi G\alpha ^2}{1+\beta }}\rho ^2,`$ (3-12)
$`\rho `$ $`=`$ $`{\displaystyle \frac{\alpha \rho _c}{r}}\mathrm{sin}\left({\displaystyle \frac{r}{\alpha }}\right),`$ (3-13)
$`\varphi `$ $`=`$ $`{\displaystyle \frac{GM}{R_s}}{\displaystyle \frac{4\pi G\rho _c\alpha ^3}{r}}\mathrm{sin}\left({\displaystyle \frac{r}{\alpha }}\right),`$ (3-14)
$`P_R`$ $`=`$ $`{\displaystyle \frac{E}{3}}={\displaystyle \frac{2\pi G\beta \alpha ^2}{1+\beta }}\rho ^2,`$ (3-15)
$`\sigma _r`$ $`=`$ $`{\displaystyle \frac{4\pi Gc\beta \rho _c\alpha ^3}{k_0(1+\beta )}}\left({\displaystyle \frac{r}{\alpha }}\mathrm{cos}\left({\displaystyle \frac{r}{\alpha }}\right)\mathrm{sin}\left({\displaystyle \frac{r}{\alpha }}\right)\right),`$ (3-16)
where $`\sigma _r=c/(3\rho D)`$ is the Rosseland mean opacity, and the density and gravitational potential are unchanged.

We ran a sequence of three simulations at different grid resolutions over half a sound crossing time and in an octant to properly specify the constant radiation boundary conditions across the star profile. The parameters are the same as the tests in §3.6, but with the addition of $`k_0=1`$ and $`\beta =1`$ to give equal significance to the radiation and hydrodynamic pressures. We find mean relative errors of 0.148, 0.0405, and 0.0138 for resolutions $`6^3`$, $`12^3`$, and $`24^3`$ cells respectively, demonstrating roughly second order convergence. Figure 14 shows spatial profiles in density for the $`24^3`$ case at the initial time (the analytic solution), and along the $`x`$, $`y`$, $`z`$ and diagonal directions running through the origin at $`t=T_s/2`$.

### 3.8 Astrophysical Jets

In this last section we perform simulations of astrophysical jets as a final test of our code. A beam of low density material ($`n_J=10^3\text{cm}^3`$) with radius $`R_J`$ is injected into a homogeneous higher density ($`n_A=10n_J`$) ambient medium along the $`+z`$ axis. Within a cylindrical radius $`\sqrt{x^2+y^2}R_J=500`$ parsecs in the $`z=0`$ plane, a constant inflow velocity of $`V_J=10^9`$ cm s-1 is maintained as a boundary condition. Simple outflow conditions are imposed at all other external boundaries. The jet material and ambient medium are initialized in pressure equilibrium, with an ambient temperature of $`10^7K`$. An ideal gas equation of state with adiabatic index $`\mathrm{\Gamma }=5/3`$ is used for both materials. The Mach number for these parameters is $`M=V_J/c_s36`$, where $`c_s`$ is the sound speed of the background material into which the beam flows.

Figure 15 shows results at time $`t=54R_J/V_J`$ from two calculations: one using a low resolution grid ($`64\times 64\times 128`$ cells) with 6 zones to cover the jet radius, and a high resolution grid ($`128\times 128\times 256`$ cells) with 12 zones/$`R_J`$. The physical box dimensions in both cases are set to $`10R_J\times 10R_J\times 20R_J`$. The results, particularly the high resolution case, clearly show all the morphological elements of astrophysical jets (Norman et al., 1982): a supersonic beam that ends in a bow shock, a cocoon composed of shock heated jet material, a working surface separating jet and shocked ambient gas, internal shock interactions, and the growth of Kelvin-Helmholtz instabilities.

Equating the ram pressure from the jet front and the equivalent pressure from the ambient medium yields

$$\overline{V}_s=\frac{V_J}{1+\sqrt{\rho _A/\rho _J}},$$
(3-17)
for an estimate of the velocity of the bow shock through the ambient medium, neglecting multi-dimensional effects. The sensitivity of the leading shock velocity $`V_s`$ to grid resolution is clear from Figure 15. In particular, we find, by tracking the bow shock position, velocities of $`V_s=0.32V_J`$ and $`0.26V_J`$ for the low and high resolution runs respectively. In comparison, equation (3-17) predicts $`\overline{V}_s=0.24V_J`$. Higher resolution allows more accurate modeling of the bow shock and working surface which effectively broadens the jet, and resolves to a greater extent 3D instabilities and internal shock interactions, all of which contribute to slowing the shock. These results are generally consistent with those of (Massaglia et al., 1996), who find jets with similar hypersonic Mach numbers and density ratios propagate at near (and greater than) unit efficiencies, defined by $`V_s/\overline{V}_s1`$.

## 4 Summary

We have developed a new multidimensional, multiphysics code (Cosmos) that can be applied to a broad range of astrophysical problems, from highly relativistic scalar-field dominated applications in early universe cosmology, to black-hole accretion and multiphase Newtonian galactic dynamics. In this contribution, we presented the radiation-chemo-hydrodynamics equations solved in the Newtonian limit, along with numerous tests to gauge the accuracy and stability of the code in various multiphysics modes. In a companion paper (Fragile et al., 2002), Cosmos is applied to the problem of supernova-enrichment in dwarf-spheroidal galaxies, which utilizes many of the capabilities (robust shock capturing, radiative cooling flows, multiphase fluids, hydrodynamic instabilities, and gravitational potentials) presented and tested here. This work also complements an earlier paper (Anninos & Fragile, 2003) in which we presented the general relativistic equations solved in Cosmos. There we provided more detailed descriptions of our numerical methods for the different energy and algorithmic formulations, together with numerical tests for highly relativistic hydrodynamical systems and black hole accretion.

This work was performed under the auspices of the U.S. Department of Energy by University of California, Lawrence Livermore National Laboratory under Contract W-7405-Eng-48.


text:
# Explosion Energies, Nickel Masses, and Distances of Supernovae of Type IIP

## 1 Introduction

Plateau Type II supernovae (SNe IIP) are believed to come from the explosion of massive supergiant stars whose envelopes are rich in hydrogen. Their light curves are easy to identify by a long plateau (sometimes up to 120–150 d) which is the result of the propagation of a cooling-and-recombination wave (CRW) through the supernova envelope that is in a state of free inertial expansion $`(u=r/t)`$. The CRW physics is discussed in detail by Imshennik & Nadyozhin (1964), Grassberg, Imshennik & Nadyozhin (1971), and Grassberg & Nadyozhin (1976). The CRW propagates supersonically downward through the expanding supernova envelope and separates almost recombined outer layers from still strongly ionized inner ones. During the plateau phase, the photosphere sits on the upper edge of the CRW front. Since the CRW downward speed turns out to be close to the velocity of the outward expansion, the photospheric radius changes only slowly during the plateau phase. If one takes into account that also the effective temperature does not change appreciably (it approximately equals the recombination temperature 5000–7000 K), the approximate constancy of the luminosity becomes obvious.

The supernova outburst properties are determined mainly by three physical parameters: the explosion energy $`E`$, the mass $``$ of the envelope expelled, and the initial radius $`R`$ of the star just before the explosion (presupernova). Litvinova & Nadyozhin (1983, 1985) have undertaken an attempt to derive these parameters from a comparison of the hydrodynamical supernova models with observations. They constructed simple approximation formulae which allow to estimate $`E`$, $``$, and $`R`$ from the observations of individual SNe IIP. Their results were confirmed by an independent semi-analytical study (Popov 1993). At that time, only one or two supernovae were sufficiently observed to apply these formulae. At present, there exist detailed observational data for 14 such supernovae, including in 12 cases expanding photosphere (EPM) distances, which we use in section 2 to estimate $`E`$, $``$, and $`R`$ by means of these formulae.

In section 3, we propose a new method of distance determination and employ it to 9 individual SNe IIP which are well-observed both at the plateau and radioactive-tail phases. The method is based on the assumption of a correlation between the explosion energy $`E`$ and the mass of 56Ni in the supernova envelope. In section 4 we compare physical parameters and distances of SNe IIP as derived from the new method with those obtained previously from the EPM method and discuss also other aspects of our results. Concluding remarks are given in section 5.

The preliminary results of this study were reported to the Workshop on physics of supernovae held at Garching, Germany, July 2002 (Nadyozhin 2003).

## 2 A comparison of hydrodynamic models with observations

Figure 1 shows a schematic SNe IIP light curve. The plateau is defined as part of the light curve on which the supernova brightness remains within $`1`$ mag of the mean value. For some supernovae, the plateau begins almost immediately after the onset of the explosion $`(t=0)`$ whereas for others a short luminosity peak can precede the plateau. The peak either appears owing to a shock wave breakout in the case of presupernovae of not very large initial radii $`(R<\mathrm{\hspace{0.25em}1000}R_{})`$ or, according to Grassberg et al. (1971), originates from the emergence of a thermal wave precursor for presupernovae of very large radii $`(R(20005000)R_{})`$ and of moderate explosion energies $`(E<1\times 10^{51}\mathrm{erg})`$, or at last it may occur as a result of interaction between the supernova envelope and a dense stellar wind (Grassberg & Nadyozhin 1987). For some SNe IIP the peak duration $`\delta t`$ lasts only a few days and is difficult to observe (shock wave breakout), for others it could be as large as 10–20 days (thermal wave or dense wind) – examples for the latter may be such supernovae as SNe 1988A, 1991al, and 1992af (see below).

It is quite clear that the middle of the plateau is to be used as the main reference point to compare the theoretical models with observations. Litvinova & Nadyozhin (1983, 1985) (LN83 and LN85, hereafter) calculated a grid of the supernova models for $`E`$, $``$, and $`R`$ within limits of $`(0.182.91)\times 10^{51}`$erg, $`(116)_{}`$, and $`(3005000)R_{}`$. They found $`E`$, $``$, $`R`$ to be strongly correlated with the plateau duration $`\mathrm{\Delta }t`$, and the middle-plateau values of the absolute $`V`$-magnitude $`M_V`$, and the expansion velocity $`u_{\mathrm{ph}}`$ at the level of the photosphere (Fig. 1). According to LN85, the following approximate relations can be used to derive $`E`$, $``$, and $`R`$ from observations:

$`\mathrm{lg}E`$ $`=`$ $`0.135M_V+2.34\mathrm{lg}\mathrm{\Delta }t+3.13\mathrm{lg}u_{\mathrm{ph}}4.205,`$ (1)
$`\mathrm{lg}`$ $`=`$ $`0.234M_V+2.91\mathrm{lg}\mathrm{\Delta }t+1.96\mathrm{lg}u_{\mathrm{ph}}1.829,`$ (2)
$`\mathrm{lg}R`$ $`=`$ $`0.572M_V1.07\mathrm{lg}\mathrm{\Delta }t2.74\mathrm{lg}u_{\mathrm{ph}}3.350,`$ (3)
where $`E`$ is expressed in units of $`10^{51}`$erg, $``$ and $`R`$ are in solar units, $`\mathrm{\Delta }t`$ in days, and $`u_{\mathrm{ph}}`$ in $`1000\mathrm{km}\mathrm{s}^1`$. Here $`M_V`$ can be expressed through the apparent $`V`$ magnitude by the relation:

$$M_V=VA_V5\mathrm{lg}(D/1\mathrm{M}\mathrm{p}\mathrm{c})25,$$
(4)
where $`D`$ is the distance to a supernova and $`A_V`$ is the total absorption on the way to the supernova. One can find from Eqs. (1)–(3) that $`E`$, $``$, and $`R`$ scale with the distance as:

$$ED^{0.675},D^{1.17},RD^{2.86}.$$
(5)
Thus, it is very important to know $`D`$ with as high accuracy as possible. We have selected 14 SNe whose observational data are collected in Table$`\mathrm{\hspace{0.17em}1}`$. The entries are: the heliocentric recession velocities $`v_0`$ (from the NED: NASA/IPAC Extragalactic Database) in column 3, the total absorption $`A_V`$ in column 4, the apparent $`V`$ magnitude of the mid-point of the plateau in column 5, the duration $`\mathrm{\Delta }t`$ of the plateau in column 6, and the photosphere expansion velocity $`u_{\mathrm{ph}}`$ in column 7. The references are in column 8.

In order to check the extrapolative capability of Eqs. (1)–(3), we have included the SN 1987A in our analysis. It is well known, that the SN 1987A presupernova radius was as small as $`50R_{}`$ – i.e. outside the interval of $`(3005000)R_{}`$ encompassed by the above equations. Moreover, the major part of the SN 1987A plateau (about 70 of 110 d) was powered by the 56Co-decay (see the review of Imshennik & Nadyozhin 1989 and references therein).

Derived properties of the 14 SNe IIP are in Table 2. Column 2 is the recession velocity $`v_{220}`$ of the supernova (column 1) corrected for a self-consistent Virgocentric infall model with a local infall vector of $`220\mathrm{km}\mathrm{s}^1`$ as described by Kraan-Korteweg (1986). Column 3 gives the distance $`D_H=v_{220}/H_0`$ assuming arbitrarily a value of $`H_0=60\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. For comparison, column 4 gives the distance $`D_{\mathrm{EPM}}`$ obtained with the use of the expanding photosphere method in the references listed at the bottom of the table. The SNe 1991al and 1992af are the exception. Owing to the incompleteness of the observational data, it is hard to determine the EPM distance to the SN 1991al (Hamuy 2001). For the same reason, the EPM distance of 55 Mpc for the SN 1992af obtained by Schmidt et al. (1994a) seems to be quite uncertain as pointed out by Hamuy (2001). For these two SNe, we present in column 4 the distances calculated by Hamuy (2001) from the CMB redshifts and the Hubble constant $`H_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. Columns 5–8 are the absolute magnitude $`M_V`$ of the mid-point of the plateau, the explosion energy $`E`$, the mass expelled $``$, and the presupernova radius $`R`$ – all derived from Eqs. (4), (1)–(3) for the $`D_H`$ distances listed in column 3. Column 9 gives the mass of 56Ni, ejected by some supernovae, which was estimated by reducing the radioactive-tail luminosities, measured by Hamuy (2001), to the distances $`D_H`$ given in column 3.

For SN 1987A, the resulting values of $`E`$ and $``$ (Table 2) differ no more than by a factor of 1.5 from current estimates based on a detailed study. However, the presupernova radius turned out to be too large. This happened because the LN85 approximations do not take into account the radioactive heating. An advanced study (Grassberg & Nadyozhin 1986) demonstrates that the radioactive heating influences only weakly on $`E`$ and $``$ furnished by Eqs. (1)–(3), whereas the $`R`$-values can be overestimated by a factor of 3. In this connection, one should have in mind that for some supernovae the $`R`$-values from Table 2 can be larger than actual presupernova radii.

According to Table 2, the resulting values of $`E`$, $``$, and $`R`$ seem to be reasonable enough: the expelled mass, explosion energy, and presupernova radius remain approximately in limits $`(1030)_{}`$, $`(0.62.6)\times 10^{51}`$ erg, and $`(2001300)R_{}`$, respectively. Hamuy (2001) assumed that SNe 1991al, 1992af were discovered several weeks after the explosion. Their plateaus, therefore, could have lasted for $`\mathrm{\Delta }t110`$d. It is quite probable, however, that their peak duration was $`\delta t20`$d for the reasons mentioned above. Having this in mind, we have chosen in Table 1 $`\mathrm{\Delta }t=90`$d which results in $`u_{\mathrm{ph}}=6000`$km s-1. In the case of $`\mathrm{\Delta }t=110`$d we would have to assume $`u_{\mathrm{ph}}=7000`$km s-1 and would obtain very large values of $`E`$ and $``$ for both supernovae: $`E7\times 10^{51}`$erg and $`40_{}`$. No other special adjustments of the observational data given in Table 1 were made.

## 3 Plateau-Tail Distance Determination

The SN IIP light curve tails are believed to be powered by the 56Co decay. The temporal behavior of the bolometric luminosity is given by (see, e.g. Nadyozhin 1994):

$$L=1.45\times 10^{43}\mathrm{exp}\left(\frac{t}{\tau _{\mathrm{Co}}}\right)\frac{_{\mathrm{Ni0}}}{_{}}\text{erg s}^1,$$
(6)
where $`t`$ is measured from the moment of explosion $`(t=0)`$, $`_{\mathrm{Ni0}}`$ is the total mass of 56Ni at $`t=0`$ which decays with a half-life of 6.10 d into 56Co, and $`\tau _{\mathrm{Co}}=111.3\mathrm{d}`$.

Equation (6) can be written in the form

$$_{\mathrm{Ni0}}=\frac{D^2}{145}Q,QF_{41}(t)\mathrm{exp}\left(\frac{t}{111.3}\right),$$
(7)
where $`_{\mathrm{Ni0}}`$ is in $`_{}`$, $`t`$ in days and $`D`$ in Mpc. The quantity $`F_{41}(t)`$ is the bolometric tail luminosity measured at time $`t`$ in units $`10^{41}\mathrm{erg}\mathrm{s}^1`$ under the assumption that the supernova is at distance $`D=1`$Mpc. Equation (7) contains a single observational parameter $`Q`$ which is independent of $`t`$ and also makes no assumption on $`D`$ as long as $`F_{41}`$ is fixed by observations.

Thus, it is irrelevant at which $`t`$ the luminosity is actually measured – one has only to be sure that the supernova really entered its tail phase. Columns $`810`$ of Table 3 give $`t`$ and corresponding values of $`F_{41}(t)`$ and $`Q`$ derived from Hamuy’s (2001) Figures 5.7 and 5.8 except SN 1999gi for which the values were calculated from the data of Leonard et al. (2002b).

If the value of $`_{\mathrm{Ni0}}`$ was known, one could easily find the distance $`D`$ from Eq. (7). So, we have to look for a way to estimate $`_{\mathrm{Ni0}}`$ independently. It seems reasonable to assume that the supernova explosion energy $`E`$ should correlate with $`_{\mathrm{Ni0}}`$ produced during the explosion. This means that

$$E=f\left(_{\mathrm{Ni0}}\right)=f\left(\frac{D^2}{145}Q\right),$$
(8)
where $`f`$ represents a statistically admissible correlation function rather than a strict mathematical relation. Inserting this expression for $`E`$ into Eq. (1) and using Eq. (4) for $`M_V`$, we obtain an equation which can be solved for $`D`$ when $`VA_V`$, $`u_{\mathrm{ph}}`$, $`\mathrm{\Delta }t`$, and $`Q`$ are known from observations. Then for given $`D`$, we can find $`E`$, $``$, $`R`$, and $`_{\mathrm{Ni0}}`$ from Eqs. (1)–(3), Eq. (4), and Eq. (7), respectively.

What can be said about the function $`f\left(_{\mathrm{Ni0}}\right)`$ at present, when the details of the SN II mechanism remain still ambiguous? First of all, it is reasonable to assume that a good fraction of $`E`$ comes from the recombination of free neutrons and protons into 56Ni just at the bottom of the envelope to be finally expelled (Nadyozhin 1978, Bethe 1996). The hydrodynamical modelling of the collapse (Nadyozhin 1978) have indicated that under favourable conditions a neutron-proton shell could be accumulated just under the steady accreting shock wave. When the mass of such a shell reaches some critical value (presumably of the order of $``$$`0.1_{}`$) the shell can become unstable in respect to recombining into the ”iron group” elements (specifically into 56Ni) to supply the stalled shock wave with the energy of $`10^{51}`$erg necessary to trigger the supernova. Here, there is a physical analogy with the origin of planetary nebulae from red giants where the energy from the recombination of hydrogen and helium causes the expulsion of a red giant rarefied envelope. The recent study (Imshennik 2002, and references therein) of the ”neutrino crown” – the region enclosed within neutrinosphere and accreting shock, turns out to be in line with such a picture of the supernova mechanism . However, some Ni can be produced through the explosive carbon-oxygen burning induced by the outgoing shock wave. In this case the energy release per unit Ni mass is lower by an order of magnitude than for the neutron-proton recombination.

The energy released by the neutron-proton recombination, producing a 56Ni mass of $`_{\mathrm{Ni0}}`$, is given by

$$E(\mathrm{np}\mathrm{Ni})=1.66\times 10^{52}\frac{_{\mathrm{Ni0}}}{_{}}\mathrm{erg}.$$
(9)
Thus, the production of only $`0.06_{}`$ of 56Ni is sufficient to provide the standard explosion energy of $`10^{51}`$erg. The current hydrodynamic models of the SN II explosions (Woosley & Weaver 1995; Rauscher et al. 2002) do not show a correlation between $`E`$ and $`_{\mathrm{Ni0}}`$ because in these models 56Ni comes from explosive silicon and carbon-oxygen burning near to the envelope bottom and its yield is sensitive to the mass cut point. The photometrical and spectroscopical properties of the SN models are virtually independent of the mass cut. On the contrary, the nucleosynthesis yields are very sensitive to the mass cut. In the current SN models the explosion is usually simulated by locating a piston at the internal boundary $`m=_{\mathrm{cut}}`$. The piston moves with time according to a prescribed law, $`R_{\mathrm{pis}}(t)`$, with the velocity $`(\dot{R})`$ amplitude being chosen to ensure the final kinetic energy of the expelled envelope of the order of $`10^{51}`$erg. There are two major uncertainties at this point. First, for a given velocity amplitude the resulting nuclear yields are still sensitive to the form of the function $`R_{\mathrm{pis}}(t)`$. Second, the presupernova structure (especially chemical composition) in the vicinity of $`m=_{\mathrm{cut}}`$ will always remain ambiguous until the detailed mechanism of the SN disintegration onto the collapsed core and thrown envelope is established. The point is that such 2D effects as rotation and large-scale mixing can result in the presupernova structure different from that predicted by the spherically symmetrical models. Under such circumstances, it is difficult to find a serious argument against the possibility to expel a noticeable amount of 56Ni from the recombination of the neutron-proton shell. Thus, we propose a neutron-proton layer which is located somewhat deeper than the value of $`_{\mathrm{cut}}`$ assumed in the current SN models. This layer recombines into 56Ni providing the energy sufficient to convert a steady-state accretion shock into the outgoing blast wave. In this case a good correlation between $`E`$ and $`_{\mathrm{Ni0}}`$ is to be expected.

The proposed correlation can have a complex nature. It is quite probable that the function $`f`$ in Eq. (8) depends also on $``$ since the supernova mechanism is expected to be sensitive to the presupernova mass. For us only the existence of some correlation is important which in combination with Eqs. (1)–(3) allows to determine the distance independently.

To demonstrate how such a method can work we make the simplest assumption that $`E`$ is proportional to $`E(\mathrm{np}\mathrm{Ni})`$. Then one can write:

$$E=\xi E(\mathrm{np}\mathrm{Ni})=\mathrm{\hspace{0.17em}16.6}\xi _{\mathrm{Ni0}}=\mathrm{\hspace{0.17em}0.1145}\xi D^2Q,$$
(10)
where, as usual, $`E`$ is in $`10^{51}\mathrm{erg}`$, $`_{\mathrm{Ni0}}`$ in $`_{}`$ and $`D`$ in Mpc. This equation implies that the function $`f`$, introduced in Eq. (8), reads as $`f(x)=16.6\xi x`$ where $`\xi `$ is an adjustable parameter which can be either less or larger than 1. If there is a noticeable contribution to $`_{Ni0}`$ from the explosive carbon-oxygen burning then $`\xi <1`$; if a noticeable contribution to the explosion energy comes from other source rather than the neutron-proton recombination then $`\xi >1`$.

Inserting $`E`$ from Eq. (10) and $`M_V`$ from Eq. (4) into Eq. (1) and solving for $`D`$, we obtain:

$`\mathrm{lg}D`$ $`=`$ $`0.374\mathrm{lg}(\xi Q)+0.0504\left(VA_V\right)+0.875\mathrm{lg}\mathrm{\Delta }t`$ (11)
$`+1.17\mathrm{lg}u_{\mathrm{ph}}2.482,`$
where $`D`$ is in Mpc, $`\mathrm{\Delta }t`$ in days, and $`u_{\mathrm{ph}}`$ in $`1000\mathrm{km}\mathrm{s}^1`$. We will refer to distances derived from Eq. (11) as ‘plateau-tail distances’, $`D_{\mathrm{P}\mathrm{T}}`$, hereafter. The results are given in Table 3 for nine supernovae selected from Table 2. We did not include SNe 1992am and 1999cr in our analysis because their last available observations may not yet reflect the radioactive tail phase. Specifically, there are only two observations of SN 1992am at the post-plateau phase of the light curve. Since the observations are separated by a short time interval of 3 days, it is difficult to derive the inclination of the bolometric light curve with a required accuracy to be shure that SN 1992am is already in the radioactive-tail phase. Moreover, one has to remember that in addition to the Co-decay the tail luminosity can also be contributed by the ejecta-wind interaction (see Chugai 1991 and references therein). SN 1992am is suspicious in this respect because its presupernova radius seems to be larger than 1000 $`R_{}`$ (Table 2). Hence, the $`_{\mathrm{Ni0}}`$-values for these SNe in Table 2 could be actually upper limits.

The different columns of Table 3 give the following quantities: (2) the distance $`D_{\mathrm{P}\mathrm{T}}`$ from Eq. (11) setting $`\xi =1`$; (3) the corresponding absolute $`V`$-magnitude of the mid-point of the plateau $`M_V`$; (4)–(7) the quantities $`E`$, $``$, $`R`$, and $`_{\mathrm{Ni0}}`$ as in Table 2, but now using the distance $`D_{\mathrm{P}\mathrm{T}}`$ as in column (1); the columns (8)–(10) are explained above.

The values of $`E`$, $``$, $`R`$, and $`_{\mathrm{Ni0}}`$ for the $`\xi `$-values different from 1 can be found using the following scaling relations which result from Eqs. (5), (7), and (11):

$$E\xi ^{0.252},\xi ^{0.438},R\xi ^{1.07},_{\mathrm{Ni0}}\xi ^{0.748}.$$
(12)
For a fixed $`Q`$, the dependence of the distance $`D_{\mathrm{P}\mathrm{T}}`$, defined by Eq. (11), on extinction $`A_V`$ proves to be very weak: an error in $`A_V`$ of $`\pm 1`$mag changes $`D_{\mathrm{P}\mathrm{T}}`$ by only $`\pm 12\%`$. However, if the tail luminosity $`F_{41}`$ is derived from the $`V`$ measurements (just the case of Hamuy’s $`F_{41}`$-values we use here) then the $`\mathrm{lg}F_{41}`$, and consequently $`\mathrm{lg}Q`$, scale as $`0.4A_V`$ and $`\mathrm{lg}D_{\mathrm{P}\mathrm{T}}`$, derived from Eq. (11), actually varies with $`A_V`$ in a standard way, as $`0.2A_V`$. If the tail luminosity were derived from infra-red measurements then the resulting $`D_{\mathrm{P}\mathrm{T}}`$ distances would be largely independent of extinction. Note also rather weak dependence on $`\xi Q`$: $`D_{\mathrm{P}\mathrm{T}}(\xi Q)^{0.374}`$. For instance, the decrease in $`\xi Q`$ by a factor of 2 results in an increase of $`D_{\mathrm{P}\mathrm{T}}`$ by 30% only.

The random errors typically of $`\pm 10\%`$ for the $`\delta t`$ and $`u_{\mathrm{ph}}`$ values assumed in Table 1 result in the uncertainty factor of $`1.2`$ for $`D_{\mathrm{P}\mathrm{T}}`$ and $`1.5`$ for $`_{\mathrm{Ni0}}(D^2)`$ given in Table 3. However, one has to keep in mind two main sources of systematic errors: $`(i)`$ probable deviation of the theoretical models (which Eqs. 13 are based on) from real SNe, and $`(ii)`$ the presentation of the $`E_{\mathrm{Ni0}}`$ correlation in the form of the straight proportionality (Eq. 10). Both the types of systematic errors are difficult to estimate at present. Although the SN models calculated in LN83 and LN85 rest upon a very simplified presupernova structure, they consistently take into account the ionization and recombination of hydrogen and helium thereby remaining still useful. When a new grid of the SN models, based on modern evolutionary presupernova structure, is created the systematic error $`(i)`$ certainly will be reduced. The reduction of the systematic error $`(ii)`$ requires a more profound knowledge of the SN mechanism. Empirically, this problem can be solved by adjusting the factor $`\xi `$ for each individual SN. It is necessary, however, to collect a much more rich statistics (at least by a factor of 3) than that available nowadays (only 9 SNe in Table 3).

## 4 Discussion

The plateau-tail distances derived in section 3 and listed in column 2 of Table 3 are plotted in a Hubble diagram in Fig. 2 (except SN 1987A which is not in the Hubble flow). The eight SNe IIP define a Hubble line with $`H_0=55\pm 5\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. Also shown in Fig. 2 are the eleven SNe IIP for which EPM distances have been published (column 3 of Table 2). They define a Hubble line of $`H_0=70\pm 4\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, i.e. the EPM distances are smaller than the plateau-tail distances by 25% on average.

At this point it is not possible to decide which of the two results is more nearly correct. Both methods, the plateau-tail distances and the EPM distances, depend on assumptions which are difficult to verify. The EPM method faces the problem of the dilution factor in an expanding atmosphere and the definition of the photospheric radius which depends on the uncertainties connected with the opacity of an expanding medium. However, it may be noted that the EPM distance of SN 1987A agrees well with the generally adopted distance of LMC of 50 kpc (Eastman, Schmidt, & Kirshner 1996) and the EPM distance of SN 1968L is indistinguishable from the Cepheid distance of NGC 5236 (M83) (Thim et al. 2003).

The main assumption which affects the plateau-tail distances concerns the nature of the proposed $`E_{\mathrm{Ni0}}`$ correlation. For our simplified example of such a correlation, all the uncertainties turn out to be cumulated in the proportionality factor $`\xi `$ between the explosion energy $`E`$ and the nickel mass $`_{\mathrm{Ni0}}`$. In Table 3 we have adopted a plausible value of $`\xi =1`$, but it cannot be excluded that $`\xi `$ is as low as 0.5 or as high as 2. Since the Hubble constant scales as $`H_0\xi ^{0.374}`$, an average value as high as $`\xi =1.9`$ would be needed to bring the plateau-tail distances in general accord with the EPM distances. Such a high average value of $`\xi `$ is, however, not supported by SNe 1987A and 1999gi. If the $`D_{\mathrm{P}\mathrm{T}}`$ distance of SN 1987A from Table 3 is scaled to the canonical LMC distance of 50 kpc, $`\xi `$ becomes 0.75. And if the host galaxy NGC 3184 of SN 1999gi with a $`D_{\mathrm{P}\mathrm{T}}`$ distance of 14.53 Mpc is a member of the same group as NGC 3198 and NGC 3319, for which Freedman et al. (2001) give a mean Cepheid distance of 13.5 Mpc, $`\xi `$ becomes 1.2. Eventually additional SNe IIP with large distances, where the influence of peculiar motions are negligible, will better determine the scatter of the Hubble diagram and allow a meaningful determination of the actual range of $`\xi `$.

We have considered three sets of the physical supernova parameters $`E`$, $``$, and $`R`$: $`(i)`$ for the Hubble distances $`D_H`$ with $`H_0=60\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ (Table 2, column 3); $`(ii)`$ for the EPM distances $`D_{\mathrm{EPM}}`$ (Table 2, column 4); $`(iii)`$ for the plateau-tail calibrated distances $`D_{\mathrm{P}\mathrm{T}}`$ (Table 3, column 2).

Although the above parameters derived from the EPM-distances are not presented in Table 2, the corresponding $`E`$ and $``$-values can be read out of Figs. 4 and 4 which compare $`E`$ and $``$ for sets $`(i)`$ and $`(ii)`$ with those for set $`(iii)`$. For seven SNe $`E`$ and $``$ are rather insensitive to the adopted distances. However for SNe 1986L and 1990E, labelled in Figs. 4 and 4, the deviations from the (P-T)-values are rather large, especially in case of the envelope mass $``$. These SNe differ from others by having a long plateau of $`(110120)\mathrm{d}`$ in combination with still a substantial expansion velocity of 4000 km s-1. As a result, their envelope masses $``$, derived from the distances defined by the $`D_H`$ and $`D_{\mathrm{EPM}}`$ values, exceed those for other SNe. Such a discrepancy for these two SNe is considerably weakened if $`\xi 2`$. Such a high value of $`\xi `$ implies that half of the explosion energy is supplied by a source different from the neutron-proton recombination. This may indicate that for massive SNe the envelope mass $``$ (in addition to $`_{\mathrm{Ni0}}`$) is involved in the correlation given by Eq. (8).

The random errors of $`E`$ and $``$ from our approximate Eqs. (1)–(3) are estimated to be about $`\pm 30`$% . Observational errors especially in the expansion velocity $`u_{\mathrm{ph}}`$ and the plateau duration $`\mathrm{\Delta }t`$ can modify $`E`$ and $``$ by another factor of 1.3. Thus it seems reasonable to assume a random uncertainty of a factor of $`1.5`$ for the individual values of $`E`$ and $``$ in Tables 2 and 3. The presupernova radii $`R`$ are very sensitive to distance errors (cf. Eq. 5) and may carry random errors of a factor of 2. The radii of SNe with large nickel masses like SN 1991al, 1992af and perhaps 1992am may carry additional systematic errors because Eqs. (1)–(3) do not take into account the radioactive heating in a consistent way.

The expelled masses $``$ are plotted against the explosion energies $`E`$ in Fig. 5 for two cases, i.e. based on EPM and plateau-tail distances. In case of the $`D_{\mathrm{P}\mathrm{T}}`$ distances the mean mass of the eight SNe IIP is $`16_{}`$ with an rms deviation of only $`3_{}`$. This narrow mass range is contrasted by a wide range of explosion energies of $`(0.52.7)\times 10^{51}`$erg. The conclusion that there is no correlation between the expelled mass – which is only $`(1.42)_{}`$ smaller than the presupernova mass – and the explosion energy is somewhat weakened by the values of $``$ and $`E`$ based on the EPM distances suggesting a marginal correlation between $``$ and $`E`$ which is mainly due to only two SNe: 1986L and 1990E.

One can think of a number of parameters which may explain the wide range of explosion energies. It could be rotation and magnetic fields inherited by the collapsing stellar core. It could be also nonspherical jet-like perturbations of a random nature arising from the macroscopic neutrino-driven advection below the accretion shock. Such perturbations could launch the outgoing blast wave earlier when the recombination nuclear energy stored in a hot neutron-proton gas was not yet as large as it should be in the case of spherical symmetry. If this is correct, one may expect that the asphericity of the explosion anticorrelates with the explosion energy.

Recently, a promising project has been started (Van Dyk et al. 1999; Smartt et al. 2001, 2002; and references therein) with the ultimate aim to identify the supernova progenitors (presupernovae) or at least to impose conclusive constraints on their masses by inspecting the prediscovery field of nearby supernovae. In particular, Smartt et al. derived upper mass limits of $`12_{}`$ and $`9_{}`$ for the progenitors of the SNe 1999em and 1999gi, assuming distances $`D`$ for the host galaxies NGC 1637 and NGC 3184 of 7.5 Mpc and 7.9 Mpc, respectively. Note that these upper limits depend on $`D`$ and have to be adjusted for other values of $`D`$ to $`12_{}`$$`\left(D/7.5\mathrm{Mpc}\right)^{0.6}`$ for SN 1999em and $`9_{}`$$`\left(D/7.9\mathrm{Mpc}\right)^{0.6}`$ for SN 1999gi. This follows from the fact that the mass-luminosity relation can be approximated as $`L^{3.3}`$ in the mass interval $`(1015)_{}`$. For SN 1999em at $`D_{\mathrm{P}\mathrm{T}}=11.08`$Mpc (Table 3) follows $`15.2_{}`$ as the upper mass limit for the SN 1999em progenitor. Hence, our result $`=15.0_{}`$ (Table 3) does not contradict to the observations as long as $`D(1999\mathrm{e}\mathrm{m})>\mathrm{\hspace{0.25em}10}`$Mpc. The situation of SN 1999gi is similar. The upper mass limit for $`D(1999\mathrm{g}\mathrm{i})=14.53`$Mpc (Table 3) is $`<9\times (14.53/7.9)^{0.6}=13.0_{}`$, i.e. not in significant contradiction with the $``$-value of $`14.5_{}`$ from Table 3. There is no contradiction either with the upper mass limit of $`15_3^{+5}_{}`$ for the SN 1999gi progenitor imposed recently by Leonard et al. (2002b).

Equations (1)–(3) by LN85, derived from a grid of 23 SNe IIP models covering a wide parameter space, imply a correlation between the absolute magnitude $`M_V`$ (and hence luminosity $`L`$ – both measured at the mid-point of the plateau) and the expansion velocity $`u_{\mathrm{ph}}`$. The correlation is shown in Figs. 7 and 7 where 23 grid models are shown by black dots; the straight lines are the least-squares fits. In Fig. 7 are also shown the eight observed SNe IIP from Table 3 marked by open circles, their absolute magnitudes $`M_V`$ (Table 3, column 2) being calculated from Eq. 4, where the plateau-tail distances $`D_{\mathrm{P}\mathrm{T}}`$ were used from Table 3, column 2. These real SNe follow about the slope of the models, but at a fixed value of $`u_{\mathrm{ph}}`$ they are fainter by $`0.6`$mag on average.

Empirically, Hamuy & Pinto (2002) have also found, using the CMB redshift-based distances, such a correlation. The slopes of their least-squares fits are virtually the same as shown in Figs. 7 and 7. Thus our models confirm their finding.

The main conclusion one can draw from Figs. 7 and 7 is that our three-parametric grid of only 23 SNe IIP properly chosen models is ample enough to reproduce the main features of the real SNe.

## 5 Conclusions

Model calculation by LN83 and LN85 of SNe IIP, leading to Eqs. (1)–(3), are combined with available EPM distances and velocity distances $`(H_0=60)`$ to derive the explosion energy $`E`$, the ejected mass $``$, and the presupernova radius $`R`$ of 14 SNe IIP. Only the apparent, absorption-corrected magnitude $`V`$ and the expansion velocity $`u_{\mathrm{ph}}`$ at the mid-point of the plateau together with its total duration $`\mathrm{\Delta }t`$ are needed as additional input parameters. The results are presented in Table 2.

Instead of using EPM or velocity distances it is also possible to use the bolometric fluxes observed during the SN IIP tail phase to determine the Ni mass and hence new, independent distances called here plateau-tail distances $`D_{\mathrm{P}\mathrm{T}}`$ (cf. Eq. 11). The $`D_{\mathrm{P}\mathrm{T}}`$ distances yield new values of $`E`$, $``$, and $`R`$ given in Table 3 for nine SNe which were observed both during their plateau and tail phases. The values of $`E`$ and $``$, based on EPM and P-T distances agree well, with the exception of SNe 1986L and 1990E whose masses $``$ coming from P-T distances are by a factor of 2 lower than from EPM distances (see Fig. 4).

The P-T distances are larger than the EPM distances by $`25`$% on average. The former suggests a value of $`H_0=55\pm 5`$. The main uncertainty of this result comes from the assumption that $`\xi =1`$, where $`\xi `$ is the ratio between the total explosion energy and the energy liberated by the neutron-proton recombination into 56Ni (cf. Eq. 10). To reduce the P-T distances to the level of the EPM distances, which correspond to $`H_0=70`$, an average value of $`\xi =1.9`$ is required. The consequence that about half of the total energy $`E`$ comes from other sources than the neutron-proton recombination into 56Ni seems rather extreme. In fact it is not supported by two SNe IIP (1987A and 1999gi) with independent distance information, which suggest that $`\xi `$ is of order of unity. Moreover, very recently Leonard et al. (2003) have obtained a Cepheid distance of $`11.7\pm 1`$Mpc to NGC 1637 – the host galaxy of SN 1999em, which is by a factor of 1.4 larger than the EPM distance (Table 2). Our $`D_{\mathrm{P}\mathrm{T}}`$ distance of 11.1 Mpc to SN 1999em (Table 3) is in a good agreement with this result. If it happens that the same factor is applicable also to the EPM distances to SNe 1986L and 1990E, there will be no need to resort to large $`\xi `$-values, such as $`\xi 2`$ (section 4), to remove the discrepancy between $`D_{\mathrm{P}\mathrm{T}}`$ and $`D_{\mathrm{EPM}}`$ for these SNe.

In conclusion we emphasize the necessity of constructing a new grid of hydrodynamic SN IIP models based on current evolutionary presupernova models and taking into account 56Ni as an additional parameter in a consistent way. Such a grid would allow to create more precise analytic approximations for a number of correlations between the physical parameters of SN IIP and their observable properties.

The ‘plateau-tail’ method of distance determination needs, of course, further critical analysis requiring a close collaboration between astronomers observing supernovae and theorists modelling their explosions. If the proposed $`E`$$`_{\mathrm{Ni}}`$ correlation is confirmed it promises to become a tool to explore the mechanism of SN II with the aid of optical and spectroscopical observations.

## Acknowledgments

It is a pleasure to express my deep gratitude to the Max-Planck-Institut für Astrophysik for hospitality and financial support. The work was supported also by the Swiss Science National Foundation (Grant 2000/061822) and the Russian Foundation for Basic Research (Project 00-02-17230).

I am grateful to G.A. Tammann for fruitful discussions and a big help and to W. Hillebrandt for useful suggestions. I thank B. Parodi for acquainting me with the FORTRAN code converting $`v_0`$ into $`v_{220}`$ and M. Hamuy for sending a copy of his Thesis where a good portion of the observational data used in this work was taken from. The anonymous referee is gratefully acknowledged for the constructive critical comments which helped to improve the paper.


text:
# Galois currents and the projective kernel in Rational Conformal Field Theory

## 1 Introduction

A most important characteristic of a (Rational) Conformal Field Theory is the associated modular representation . Not only does it prescribe the transformation of the genus one characters of the primary fields under the modular group $`\mathrm{\Gamma }\left(1\right)=SL(2,)`$, but it does also determine the central charge (modulo 8), the fractional part of the conformal weights, the fusion rules via Verlinde’s formula, etc.. All these properties are linked to the fact that the modular representation of an RCFT provides part of the defining data of a modular tensor category .

A fundamental property of the modular representation, conjectured by many authors over the years  and finally proved in , is the congruence subgroup property: there exists a natural number $`N`$ (which turns out to equal the order of the matrix $`T`$ representing the Dehn-twist $`\tau \tau +1`$), such that the kernel of the modular representation contains the principal congruence subgroup

$$\mathrm{\Gamma }\left(N\right)=\left\{\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\mathrm{\Gamma }\left(1\right)\right|a,d1\left(\mathrm{mod}N\right),b,c0\left(\mathrm{mod}N\right)\}.$$

In other words all modular transformations which belong to $`\mathrm{\Gamma }\left(N\right)`$ are represented by the identity. This result has far reaching implications both in theory and practice, e.g. it is the basis of powerful algorithms to compute arbitrary modular matrices.

Once we know that the kernel is a congruence subgroup, we can look for a complete description of it in terms of a small set of easy-to-compute invariants. As it turns out, this won’t work for the kernel, but the reason is that the kernel is not the right notion to look at. It is the projective kernel, i.e. the set of modular transformations represented by scalar matrices, that is the natural object, and for it one does indeed get a simple description. Along the way emerges the notion of Galois currents, which are simple currents  related to the Galois action . The basic result about Galois currents is uniqueness: a model can have at most one non-trivial Galois current. This leads at once to the partition of all RCFTs into two classes, according to whether they have a non-trivial Galois current or not. As we shall see, while the latter are the generic ones, many important examples of RCFTs - like the Ising model - fall into the first class.

This paper aims to give a brief survey of the above results, concentrating on the conceptual issues, and leaving the technical details to a separate publication. Besides introducing the relevant notions and stating the main results, a couple of simple examples are included in order to illustrate the general theory.

## 2 Modular matrices and the Galois action

As indicated in the introduction, to each rational CFT is associated a finite dimensional representation $`D`$ of the modular group $`\mathrm{\Gamma }\left(1\right)`$. This representation is most conveniently described by a pair of matrices $`S`$ and $`T`$, which represent the generators $`\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)`$ and $`\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)`$ of the modular group. These matrices have to satisfy the defining relations

$`STS`$ $`=`$ $`T^1ST^1,`$
$`S^4`$ $`=`$ $`1`$
of the modular group. It should be noted that the matrices $`S`$ and $`T`$ carry much more information than a mere linear representation of $`\mathrm{\Gamma }\left(1\right)`$, for the individual matrix elements have an invariant meaning in themselves, in other words the representation comes equipped with a distinguished basis formed by the genus 1 characters of the primary fields. The properties of the matrix elements with respect to this distinguished basis are summarized in Verlinde’s theorem :

1. The matrix $`T`$ is diagonal and of finite order;
2. The matrix $`S`$ is symmetric;
3. The fusion rules of the primaries are given by Verlinde’s formula

$$N_{pqr}=\underset{s}{}\frac{S_p^sS_q^sS_r^s}{S_0^s},$$

where $`0`$ labels the vacuum (the conformal block of the identity).

One may show that the modular representation is quite special in many respects, the basic result being the congruence subgroup property mentioned in the introduction, i.e. that the kernel $`\mathrm{ker}D`$ of the modular representation is a congruence subgroup of level $`N`$ . One may also show that the *conductor* $`N`$, which is equal to the order of the matrix $`T`$, is bounded by a universal function of the dimension (= number of primary fields). Moreover, if we denote by $`K`$ the projective order of $`T`$, i.e. the smallest positive integer such that $`T^K`$ is a scalar matrix, then the integer $`e=\frac{N}{K}`$ is a divisor of $`12`$ . The integers $`e`$ and $`K`$ are important characteristics of the model.

Because the individual matrix elements have an invariant meaning, it is meaningful to look at their arithmetic properties, which is the content of the theory of the Galois action . The basic idea is to look at the field $`F`$ obtained by adjoining to the rationals the matrix elements of all modular matrices (it is enough to adjoin those of the generators $`S`$ and $`T`$). It turns out that the field $`F`$ equals the cyclotomic field $`\left[\zeta _N\right]`$, where $`N`$ is the conductor of the model. By a well known result of algebraic number theory, the Galois group of the cyclotomic field $`\left[\zeta _N\right]`$ over $``$ is isomorphic to the group $`𝒢_N=\left(/NZ\right)^{}`$ of prime residues modulo $`N`$, where to the prime residue $`l𝒢_N`$ corresponds the map $`\sigma _l`$ sending $`\zeta _N`$ to its $`l`$-th power $`\zeta _N^l`$. As $`\sigma _l`$ maps $`F`$ into itself, it is meaningful to consider its action on the matrix elements of $`S`$ and $`T`$. The basic result , which follows from Verlinde’s theorem, is

$$\sigma _l\left(S_q^p\right)=\epsilon _l\left(q\right)S_{\pi _lq}^p,$$

where $`\pi _l`$ is a permutation of the primaries - the *Galois permutation* associated to $`l𝒢_N`$ -, while $`\epsilon _l\left(q\right)=\pm 1`$ is a sign. The above Galois permutations define a permutation action of $`𝒢_N`$ on the set of primaries, because they satisfy $`\pi _{lm}=\pi _l\pi _m`$.

## 3 Galois currents

According to the traditional definition , a simple current is a primary field $`\alpha `$ whose quantum dimension is $`1`$, i.e. such that $`S_0^\alpha =S_0^0`$. From this definition follows that the fusion product of a simple current $`\alpha `$ with a primary $`p`$ is again a primary field, denoted $`\alpha p`$, i.e. simple currents induce permutations of the primaries. Of course, the permutations arising this way are quite special, and it is a challenging problem to characterize those permutation groups that may correspond to the action of simple currents on the primary fields of some RCFT. In any case, this induced permutation action allows us to consider for each simple current $`\alpha `$ the set

$$L_\alpha =\left\{l𝒢_N\right|\alpha p=\pi _lp\},$$

i.e. the set of prime residues mod $`N`$ such that the Galois permutation $`\pi _l`$ is the same as the permutation induced by $`\alpha `$ (this set may be empty). We call a simple current $`\alpha `$ a *Galois current* when $`L_\alpha \mathrm{}`$. The set of Galois currents may be shown to be closed under the fusion product, in other words they form a group.

Note that the vacuum is always a Galois current, as it is a simple current for which the corresponding permutation is the identity, consequently $`1L_0`$. The interesting question is whether there exists models with non-trivial simple currents, so let’s give a couple of simple examples.

Our first example is the Ising model, the Virasoro minimal model $`(4,3)`$ of central charge $`c=\frac{1}{2}`$. As it is well known, this model has 3 primaries: the vacuum $`0`$, the energy operator $`\epsilon `$, and the spin operator $`\sigma `$. One has $`K=16`$ and $`e=3`$ in this model. A simple computation reveals that $`L_\epsilon =\{5,11,13,19,29,35,37,43\}`$, i.e. $`\epsilon `$ is a non-trivial Galois current. More generally, one may show that the only Virasoro minimal models with a non-trivial Galois current are those of the form $`(4,p)`$ with $`p`$ odd - the Ising model corresponding to $`p=3`$ -, in which case the primary field with Kac-label $`(3,1)`$ is the non-trivial Galois current (for $`p>3`$).

There are many more examples of theories with non-trivial Galois currents, e.g. some of the Ashkin-Teller models (= $`_2`$ orbifolds of the compactified boson, the Galois current being the marginal perturbation), many WZNW models, products and orbifolds of the above, etc. As these examples show, Galois currents appear in many types of RCFTs. Closer examination of the above examples leads to the following observations:

* the group of Galois currents is small (in all of the above examples its order is either 1 or 2);
* theories with non-trivial Galois currents are sparse: in two-parameter families of models they form one parameter subfamilies.

As it turns out, the above observations reflect general properties of RCFT, for one may show - the proof is quite lengthy and technical \- that

1. If there is a non-trivial Galois current, then it is unique, i.e. the group of Galois currents has either 1 or 2 elements;
2. If a theory has a non-trivial Galois current, then $`e`$ is odd and $`K`$ is a multiple of $`16`$.

With respect to this last result, we note that the reverse implication is by no means true, for there exist examples of RCFTs with odd $`e`$ and $`K`$ a multiple of $`16`$, but no non-trivial Galois current. On the other hand, in those classes of RCFTs for which one has an explicit expression for $`e`$ and $`K`$, those with odd $`e`$ and $`K`$ a multiple of $`16`$ are indeed sparse, in accord with the second observation above.

As to the first result, it exhibits an important dichotomy for RCFTs: a model has either one non-trivial Galois current, or it has none. As we shall see in the next section, the existence of a non-trivial Galois current has important consequences regarding the structure of the modular representation, so the above mentioned dichotomy is not only meaningful, but also relevant. In this respect one might say that the existence of a non-trivial Galois current for the Ising model is a good signal of its special nature.

## 4 Ordinary vs projective kernel

Given a linear representation of a group, it is natural to look at the kernel of the representation, i.e. the set of group elements that are represented by the identity operator. As we have mentioned in section 2, in case of the modular representation associated to a RCFT, the kernel is a congruence subgroup whose level equals the order $`N`$ of $`T`$. Knowing this fundamental result, the next natural step is to look for a simple description of the kernel. As it turns out the natural object to look at is not the kernel itself, but rather the *projective kernel* $`𝒦=\left\{m\mathrm{\Gamma }\left(1\right)\right|D\left(m\right)=\xi \left(m\right)1\}`$, which consists of those modular transformations that are represented by scalar matrices, where $`\xi :𝒦`$ is a linear character of $`𝒦`$, which we term the *central character*.. Let’s explain why this is so.

Consider a composite system denoted $`𝒞_1𝒞_2`$, which is made up of two independent and non-interacting subsystems $`𝒞_1`$ and $`𝒞_2`$. Note that if both subsystems $`𝒞_1`$and $`𝒞_2`$ are RCFTs, then so is their composite $`𝒞_1𝒞_2`$. As the subsystems are completely independent, one expects that it should be possible to determine unambiguously the value of any ”natural” quantity for the composite system from the knowledge of the corresponding values for the subsystems.

In this respect, the partition function and all correlators are ”natural”, for they factorize in the composite system into a product of corresponding quantities for the subsystems. Similarly, the modular representation is ”natural”, for the modular representation associated to $`𝒞_1𝒞_2`$ is the tensor product of the modular representations associated to $`𝒞_1`$ and $`𝒞_2`$. On the other hand, the kernel is not ”natural”, for one cannot determine the kernel of a tensor product from the sole knowledge of the kernels of the factors. The ”natural” notion is that of the projective kernel, because the projective kernel of a tensor product is simply the intersection of the projective kernels of the factors. This is the conceptual explanation for looking at the projective kernel instead of the ordinary one. We note that the same ”naturality” argument explains why it is the parameter $`K`$ and not the conductor $`N`$ that enters most of our results, for the former is ”natural” in the above sense, while the later is not.

After these preliminaries, let’s describe the structure of the projective kernel. To this end, let’s first introduce the notation

$$\mathrm{\Gamma }(K,𝔤)=\left\{\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\mathrm{\Gamma }\left(1\right)\right|a,d𝔤,b,c0\left(\mathrm{mod}K\right)\}$$

for a subgroup $`𝔤<𝒢_K`$. Clearly, $`\mathrm{\Gamma }(K,\left\{1\right\})`$ equals the principal congruence subgroup $`\mathrm{\Gamma }\left(K\right)`$, from which follows that $`\mathrm{\Gamma }(K,𝔤)`$ is a congruence subgroup for any $`𝔤<𝒢_K`$.

Remember the Galois permutations $`\pi _l`$ from section 2, and consider the kernel of the permutation action $`\pi `$, i.e. the set $`\mathrm{ker}\pi =\left\{l𝒢_N\right|\pi _l=1\}`$. Clearly, $`\mathrm{ker}\pi `$ is a subgroup of $`𝒢_N`$, and reducing all elements of $`\mathrm{ker}\pi `$ modulo $`K`$ we get a subgroup $`𝔥`$ of $`𝒢_K`$. It does follow from general properties of the Galois action that the exponent of $`𝔥`$ is $`2`$, in other words $`l^21\left(\mathrm{mod}K\right)`$ for all $`l𝔥`$ . The basic result about the projective kernel is that $`\mathrm{\Gamma }(K,𝔥)`$ is a normal subgroup of $`𝒦`$ \- in particular $`𝒦`$ is a congruence subgroup of level $`K`$ -, and there is an explicit isomorphism between the factor group

$$𝒦/\mathrm{\Gamma }(K,𝔥)$$

and the group of Galois currents! This means that the knowledge of $`K`$, the subgroup $`𝔥`$ and the group of Galois currents does completely determine the projective kernel. Note that the index of the principal congruence subgroup $`\mathrm{\Gamma }\left(K\right)`$ in $`𝒦`$, which equals the number of Galois currents times the order of $`𝔥`$, is always a power of $`2`$ by the above, a result that seems fairly non-trivial.

As an example, let’s consider once again the Ising model $`(4,3)`$. As discussed in section 3, one has $`e=3`$ and $`K=16`$, and there is a non-trivial Galois current (the energy operator $`\epsilon `$). In this case, the group $`𝔥`$ consists of the prime residues $`\{1,7,9,15\}`$ modulo $`16`$. As to the projective kernel, it has the coset decomposition

$$𝒦=\mathrm{\Gamma }(K,𝔥)\left(\begin{array}{cc}3& 8\\ 8& 11\end{array}\right)\mathrm{\Gamma }(K,𝔥)$$

## 5 Discussion

As we have seen, the notion of Galois currents is not only meaningful in the sense that there exists non-trivial examples, but also relevant to the analysis of the properties of RCFTs. While we have concentrated on their impact on the structure of the projective kernel, it is quite plausible that they play an important role for other aspects of the theory as well, e.g. fusion rules and modular invariants. Elucidating these connections seems to be a rewarding task for the future.

The simple characterization of the projective kernel described in section 4 should be regarded only as a first step. It should be followed by an understanding of the central character, the map $`\xi :𝒦`$ appearing in the definition of the projective kernel. Another interesting point would be to find out which subgroups $`𝔥<𝒢_K`$ are allowed for a given value of $`K`$.

Finally, one might speculate on the relevance of the above for the classification of RCFTs. The knowledge of the projective kernel might be an important ingredient of classification attempts. From another point of view, one may hope that the special properties of RCFTs with non-trivial Galois currents could lead to a classification of this special class of model.


text:
# XAssist: A System for the Automation of X-ray Astrophysics Analysis

## 1. Introduction

We are currently in a renaissance for X-ray astronomy, with two major missions, Chandra and XMM-Newton, currently operating. These missions are producing large amounts of archival data, which is supplementing existing databases from missions such as ROSAT and ASCA. Historically, only X-ray “experts” usually attempted the analysis of X-ray data. This is because there were fundamental differences in the analysis of X-ray data compared with other bandpasses, most notably the fact that individual photons are detected as opposed to the accumulation of (only) spectra or images. Most modern X-ray detectors are imaging spectrometers so each observation results in a photon list from which images, spectra and light curves can be extracted. In general the numbers of photons are small so Poissonian rather than Gaussian statistical methods must be used. The spectral and spatial resolution of most detectors is moderate at best and forward-fitting convolution methods are needed to properly fit the data. All of these factors limit the accessibility of the X-ray data to non-experts. In addition each mission tends to have its own unique software package for the reduction and analysis of the data, and X-ray data often require reprocessing as the calibration improves. These latter two factors also limit the ability of experts to take advantage of all available data for a given project, particularly large-scale surveys.

We have developed a software package to address these concerns. XAssist is capable of performing data reduction and preliminary analysis for ROSAT, Chandra and XMM-Newton data. It is fully automatic making it well-suited for surveys, as well as for the reprocessing of existing data. Below we will discuss its capabilities and prospects for the future.

## 2. Capabilities

XAssist currently has the following capabilities:

* Downloads data
* Reprocesses data
* Creates exposure maps and detector masks (if possible)
* Detects sources (using built-in routine for ROSAT and ASCA, CIAO wavdetect for Chandra, etc.)
* Fits each source with ”simple” (i.e., not including point-spread function) model to establish source extent and (Poisson-correct) significance (using the stand-alone python program ximgfit111http://www.xassist.org/xassist/Download.jsp;)
* Flags extended, confused and problematic sources
* Computes median (or mean) background level
* Excludes times of high background
* Extracts spectra, “postage stamp” images, and light curves of each source for more detailed analysis (a simple power-law model is fit to sources with more than 100 sources)
* Optionally performs chip-by-chip analysis (relevant just for Chandra right now)
* Analysis can be restricted to an energy band
* Large emphasis on detailed reporting
* Looks for correlations of X-ray sources with astronomy databases and provides links on the detailed source web reports to query Simbad (see below)

Note that in the case of Chandra analysis, most of the data reduction steps are based on the “threads” reported in the Chandra web site222http://cxc.harvard.edu;. Most parameters controlling XAssist are read from IRAF/FTOOLS-style parameter file, and can be set on the command-line as well (allowing for the automated setup of XAssist for surveys). XAssist can be run (and configured) interactively (with a text-based interface).

## 3. Sample Output

Figure 1 shows images created as part of the report for the processing of the Chandra observation NGC 1569. Figure 2 shows the detailed report generated for a source. While there are admittedly “warts” that occur (as in any automated system), this example demonstrates that even in moderately crowded fields the system performs well and continues to the point of fitting a power-law spectrum model to the source spectrum. Obviously, this opens a powerful possibility for virtual observatories, namely that searches could be performed on high-level quantites such as spectral slope. Even though human inspection would of course still be necessary for science-grade results, an automated system such as this could cull samples and produce useable results for many of the sources, both of which would be particularly useful prior to observing proposal deadlines (especially if the data of interest had only recently become publically available and had not been published yet). Queries are also submitted to HEASARC to find correlations of X-ray sources with 2$`\mu `$ass, USNO, FIRST, and other catalogs.

## 4. Pipelines

To demonstrate XAssist’s capabilities, pipelines333http://www.xassist.org/xassist/Pipelines.jsp; have been set up to run XAssist on ROSAT HRI and Chandra observations of RC3 galaxies. The ROSAT analysis is nearing completion, and crontab jobs are checking for the public release of Chandra observations that may containing RC3 galaxies in the FOV of view, and these data sets are downloaded (using the script cda.py also available for download). XMM-Newton has recently been added as a supported mission, and we are in the process of establishing a pipeline for public XMM data. The pipeline products are searchable and work is in progress to allow users to request fields to be added to the processing lists.

## 5. Current and Future Work

* Finishing off details (especially for ASCA and ROSAT PSPC analysis)
* Better installation and configuration support (including the creation of a GUI)
* Web access to XAssist
* Web services: limited status reporting is already available using SOAP and it will be possible to request full or partial processing of data via SOAP requests