text:
# Untitled Document

Three-Dimensional Modeling of the Deflagration Stage

of a Type Ia Supernova Explosion

A. M. Khokhlov

Laboratory for Computational Physics and Fluid Dynamics,

Naval Research Laboratory, Washington, DC 20375

Abstract

The paper describes a physical model and numerical algorithm for modeling Type Ia supernova (SNIa) explosions in three-dimensions and presents first results of modeling a deflagration explosion in a non-rotating, Chandrasekhar-mass carbon-oxygen (CO) white dwarf. Simulations show that the turbulent flame speed grows exponentially, reaches $`30`$% of the speed of sound, and then declines as the large-scale turbulence is frozen by expansion. The freezing of turbulent motions appears to be a crucial physical mechanism regulating the rate of deflagration in SNIa. The energy of the explosion is comparable to that of a typical SNIa. However, the presence of the outer layer of unburned CO and the formation of intermediate mass elements and pockets of unburned CO near the center pose problems for the modeling of SNIa spectra. Delayed detonation is a way to alleviate these problems and to produce consistent spectra.

1. Introduction

Type Ia Supernovae (SNIa) are thought to be thermonuclear explosions of Chandrasekhar-mass carbon-oxygen white dwarfs (CO-WD). They play an important role in nucleosynthesis (Nomoto et al. 1984) and the balance of energy in interstellar medium, and they are important distance indicators used in cosmology (Perlmuter et al., 1999; Riess et al. 2000). Our understanding of SNIa explosions is far from complete. To predict the explosion outcome, one needs to know how thermonuclear burning propagates inside the exploding star. The mechanisms and the speed of thermonuclear burning in SNIa remain an unsolved theoretical problem. A recent review by Hillebrandt & Niemeyer (2000) contains a large list of relevant publications.

Burning can propagate in a supernova either as a flame caused by the electron heat conduction or as a detonation in which reactions are triggered by a shock. Both the detonation speed and the laminar flame speed $`S_l`$ in SNIa are known. However, in the gravitational field of the star, the laminar flame becomes unstable with respect to the Rayleigh-Taylor (RT) instability (Nomoto et al. 1976). This leads to a turbulent deflagration with some unknown turbulent flame speed $`S_t>S_l`$.

Observations indicate that burning of a Chandrasekhar-mass CO-WD must start as a deflagration. Spectra of SNIa at maximum light show the presence of large amounts of intermediate mass elements such as Si, S, Mg (Pskovskii 1977, Branch 1981). These elements can only be synthesized at densities less than $`10^7`$g/cc, where the temperature of burned products stays below $`(45)\times 10^9`$K. When burning starts as a deflagration, the initially high-density WD expands while it burns, its density decreases, and this creates the necessary conditions for the production of intermediate mass elements. In one-dimensional deflagration models, $`S_t`$ is a free parameter (Nomoto et al. 1976, Woosley & Weaver 1986). Delayed detonation models assume that a deflagration undergoes a transitions to a detonation at low densities (Khokhlov 1991). In these models, both $`S_t`$ and the moment of a deflagration-to-detonation transition (DDT) are free parameters. To fit observations, burning must be initially slow in order to expand the WD significantly. Then burning must accelerate quickly to incinerate expanding outer layers rapidly. This requires very carefully tuning $`S_t`$ as a function of time in deflagration models. Delayed-detonation models accomplish the incineration by switching to a supersonic detonation mode of burning.

Two-dimensional simulations of SNIa explosions show that burning leads to a significant expansion of a WD. However, the deflagration speed found in all of the simulations is too small to lead to a powerful explosion (Livne, 1993; Arnett & Livne, 1994; Reinecke et al., 2000). Whether this is an inherent difficulty of the deflagration model or simply a deficiency of two-dimensional simulations is a question. Two-dimensional simulations lack a key physical ingredient of turbulent burning – turbulent energy cascade that is directed in three dimensions from large to small spatial scales. Simulations of burning in a vertical column in a uniform gravitational field show that the turbulent flame speed grows faster in three dimensions and eventually reaches a turbulent steady state where $`S_t`$ becomes independent of $`S_l`$ and of the details of burning on small scales (Khokhlov, 1995; hereafter K95). In two dimensions, $`S_t`$ grows more slowly and depends on the value of $`S_l`$. However, turbulence in the exploding star is affected by expansion on spatial scales $``$, where the characteristic RT timescale $`\tau _{RT}(g)^{1/2}`$ is comparable or exceeds the expansion timescale $`\tau _eR/U_e`$, where $`R`$ is stellar radius and $`U_e`$ is the (time-dependent) expansion velocity. Expansion tends to freeze large-scale turbulent motions and thus decreases the turbulent flame speed (K95). The actual $`S_t`$ in a supernova must depend on the competition of the RT instability, turbulent energy cascade, and expansion. To predict $`S_t`$, these effects must be modeled in three dimensions. This paper describes a physical model and numerical algorithm for modeling SNIa explosions in three-dimensions and presents first results of modeling a deflagration explosion in a nonrotating, Chandrasekhar-mass CO-WD.

2. Physical and numerical model

2.1. Fluid dynamics

Stellar matter is described by the Euler equation for an inviscid fluid

$$\begin{array}{cc}\hfill \frac{\rho }{t}& =\left(\rho 𝐔\right),\hfill \\ \hfill \frac{\rho 𝐔}{t}& =\left(\rho \mathrm{𝐔𝐔}\right)P+\rho 𝐠,\hfill \\ \hfill \frac{E}{t}& =\left(𝐔\left(E+P\right)\right)+\rho 𝐔𝐠+\rho \dot{q},\hfill \end{array}$$
$`(1)`$
where $`\rho `$, $`E=E_i+\rho U^2/2`$, $`E_i`$, $`𝐔`$, $`𝐠`$ and $`\dot{q}`$ are the mass density, energy density, internal energy density, velocity of matter, gravitational acceleration, and nuclear energy release rate per unit mass, respectively. The acceleration of gravity is computed using an angle-averaged density distribution. The equation of state (EOS) of degenerate matter includes contributions from ideal Fermi-Dirac electrons and positrons, equilibrium Planck radiation, and ideal ions. Pressure $`P=P(\rho ,E_i,Y_e,Y_i)`$ and temperature $`T=T(\rho ,E_i,Y_e,Y_i)`$ are determined by the EOS as functions of $`\rho `$, $`E_i`$, electron mole fraction $`Y_e`$, and the mean mole fraction of ions $`Y_i`$.

2.2. Flame propagation

The flame is advanced using a flame-capturing technique which mimics a flame with a prescribed normal speed (K95). A scalar variable $`f`$, such that $`f=0`$ in the unburned matter and $`f=1`$ in the material which has passed through the flame, obeys a reaction-diffusion equation

$$\frac{f}{t}+𝐔f=K^2f+R,$$
$`(2)`$
with artificial reaction and diffusion coefficients

$$\begin{array}{cc}\hfill R& =\{\begin{array}{cc}C=\mathrm{const}.,\hfill & \text{if }f_0f1\text{;}\hfill \\ 0,\hfill & \text{otherwise,}\hfill \end{array}\hfill \\ \hfill K& =\mathrm{const},\hfill \end{array}$$
$`(3)`$
with some $`0<f_0<1`$. In this paper, $`f_0=0.3`$. Equation (2) has a solution $`f(xSt)`$ which describes a reaction front that propagates with the speed $`S=(KC/f_0)^{1/2}`$ and has a thickness $`\delta (K/C)^{1/2}`$ (see Appendix of (K95)). If $`K`$ and $`C`$ in (3) are set to

$$K=S(\beta \mathrm{\Delta }r)\sqrt{f_0},C=\frac{S}{(\beta \mathrm{\Delta }r)}\sqrt{f_0},$$
$`(4)`$
with $`\beta =\mathrm{constant}`$, the front propagates with a prescribed speed $`S`$ and spreads onto several computational cells of size $`\mathrm{\Delta }r`$. The choice of $`\beta =1.5`$ spreads the flame on $`34`$ cells. Making the front narrower is not practical, since the fluid dynamics algorithm spreads contact discontinuities on $`4`$ cells. The energy-release rate inside the front is defined as

$$\dot{q}=q_f\frac{df}{dt},$$
$`(5)`$
where $`q_f`$, defined below in Section 3.4, is the nuclear energy release inside the front. Due to the energy release, both density and velocity vary across the front, and fluid motions are generated which advect the flame relative to the computational mesh. The flame-capturing technique described above advances the reaction front relative to the fuel with the speed practically independent of the orientation of the front on the mesh, fluid motions, and resolution.

2.3. Burning on small scales

A turbulent flame in a vertical column of width $``$ subjected to an acceleration of gravity $`g`$ propagates in a steady-state with the speed $`S_t0.5\sqrt{Ag}`$ independent of the laminar speed $`S_l`$, where $`A=(\rho _0\rho _1)/(\rho _0+\rho _1)`$ is the Atwood number, and $`\rho _0`$, $`\rho _1`$ are the densities ahead and behind the flame front, respectively (K95). The time to reach a steady state is $`2/S_t`$, and the thickness of the flame brush is $`2`$. To describe burning on scales that are not resolved in the simulations, the value of $`S`$ in (4) is taken as

$$S=\mathrm{max}(S_l,S_{sub}),$$
$`(6)`$
with

$$S_{sub}=0.5\sqrt{A\alpha \mathrm{\Delta }r\left(s_1|𝐠|+s_2\mathrm{max}(0,𝐧𝐠)\right)}$$
$`(7)`$
based on the assumption that burning on small scales is driven by the RT instability at these scales and that at $`<<R_{WD}`$ burning adjusts to a local steady state (Livne 1993, Arnett & Livne 1994, K95). In (7), $`\alpha 1`$ determines the driving scale, $`\alpha \mathrm{\Delta }r`$, in terms of the computational cell size $`\mathrm{\Delta }r`$, and $`𝐧=f/f`$ is the unit vector normal to the flame front and directed towards the products of burning. The choice of $`s_1=1`$, $`s_2=0`$ gives a subgrid flame speed that is independent of the orientation of the flame front relative to the direction of gravity. The choice $`s_1=0`$, $`s_2=1`$ gives a directional-dependent speed that reduces to $`S_l`$ when the flame moves in the direction of gravity. Finally, $`s_1=s_2=0`$ gives the laminar flame speed. In this paper, only the limiting, angle-independent and laminar prescriptions for $`S`$ were used. The Atwood number $`A(\gamma 1)\overline{\rho }q_f/\mathrm{\hspace{0.17em}2}\gamma P`$ was calculated from conservation of enthalpy, assuming constant pressure across the flame, where $`\gamma `$ and $`\overline{\rho }=(\rho _1+\rho _0)/2`$ are average values of the adiabatic index and mass density. Under the conditions of interest here, $`A0.10.4`$. The laminar flame speed $`S_l`$, taken according to Timmes & Woosley (1992) and Khokhlov et al. (1997), is $`S_l(10^210^4)\times a_s`$ for $`\rho 10^910^7`$g/cc, where $`a_s5\times 10^8`$cm/s is the sound speed.

2.4. Nuclear kinetics

There are three distinct stages of carbon burning. First, the $`{}_{}{}^{12}C+^{12}C`$ reaction leads to the consumption of $`{}_{}{}^{12}C`$ and formation of mostly $`{}_{}{}^{20}Ne`$, $`{}_{}{}^{24}Mg`$, protons and $`\alpha `$-particles. Then the onset of the nuclear statistical quasi-equilibrium (NSQE) takes place, during which Si-group (intermediate mass) elements are formed. Finally, Si-group elements are converted into the Fe-group elements and the nuclear statistical equilibrium (NSE) sets in. The NSQE and NSE relaxation involves hundreds of species from carbon to zinc and thousands of reactions of these nuclei with protons, neutrons, and $`\alpha `$-particles. The timescales of all these stages strongly depend on temperature. The timescale of carbon consumption is much shorter than the quasi-equilibrium timescale $`\tau _{nsqe}`$ which, in turn, is much shorter than the NSE timescale $`\tau _{nse}`$.

A four-equation kinetic scheme is used here to describe the energy release, synthesis of Si- and Fe-group elements, and the neutronization of NSE matter. A similar scheme has been used in one-dimensional simulations of SNIa (Khokhlov 1991). The kinetic equation for the mole fraction of carbon $`Y_C`$

$$\frac{dY_C}{dt}=\rho A(T_9)\mathrm{exp}(Q/T_{9a}^{1/3})Y_C^2,$$
$`(8)`$
describes carbon consumption through the major reactions $`{}_{}{}^{12}C({}_{}{}^{12}C,p)^{23}Na(p,\gamma )^{24}Mg`$ and $`{}_{}{}^{12}C({}_{}{}^{12}C,{}_{}{}^{4}He)^{20}Ne`$ with the branching ratio $`1`$. Here $`Q=84.165`$, $`T_{9a}=T_9/(1+0.067T_9)`$, where $`A(T_9)`$ is a known function (Fowler et al. 1975).

Most of the nuclear energy is released during the carbon exhaustion stage and the subsequent synthesis of the Si-group nuclei. The energy release or consumption due to the transition from Si-group to Fe-group nuclei (NSE relaxation) is less than 10%. Therefore, the nuclear energy release rate is approximated as

$$\frac{dq_n}{dt}=Q_C\frac{dY_C}{dt}+\frac{q_{nse}q_n}{\tau _{nsqe}},$$
$`(9)`$
where $`q_n(t)`$ is the binding energy of nuclei per unit mass, $`Q_C=4.48\times 10^{18}`$ ergs g-1 mol-1 is the energy release due to carbon burning (8), and $`q_{nse}(\rho ,T,Y_e)`$ is the binding energy of matter in the state of NSE.

The equation

$$\frac{d\delta _{nse}}{dt}=\frac{1\delta _{nse}}{\tau _{nse}},$$
$`(10)`$
traces the onset of NSE, where $`\delta _{nse}=0`$ in the unburned matter and $`\delta _{nse}=1`$ in the NSE products. Intermediate mass elements are expected where $`Y_C(t=\mathrm{})0`$ and $`\delta _{nse}(t=\mathrm{})<1`$. Fe-group elements are expected where $`Y_C(t=\mathrm{})0`$ and $`\delta _{nse}(t=\mathrm{})1`$. The ‘e-folding’ NSQE and NSE timescales

$$\begin{array}{cc}\hfill \tau _{nsqe}& =\mathrm{exp}\left(149.7/T_939.15\right)\mathrm{s},\hfill \\ \hfill \tau _{nse}& =\mathrm{exp}\left(179.7/T_940.5\right)\mathrm{s},\hfill \end{array}$$
$`(11)`$
approximate the results of the detailed calculations of carbon burning. The original kinetic scheme (Khokhlov 1991) also contained the equation for the mean ion mole fraction $`Y_i`$. To simplify the scheme in this paper, variations of $`Y_i`$ due to nuclear reactions were neglected since ions make a small contribution to the equation of state. A constant $`Y_i=0.07`$ was used instead, which is an average of $`Y_i`$ in the unburned and typical NSE matter in the SNIa explosion conditions.

Neutronization is described by the equation for the electron mole fraction $`Y_e`$,

$$\frac{dY_e}{dt}=R_w(\rho ,T,Y_e).$$
$`(12)`$
The corresponding term describing neutrino energy losses, $`\dot{q}_w(\rho ,T,Y_e)`$, is added to the energy conservation equation in (1). Values of $`q_{nse}`$, $`\dot{q}_w`$, and $`R_w`$ were computed assuming NSE distribution of individual nuclei, as described in (Khokhlov 1991). Recently, there has been an important development in theoretical computations which shows significantly smaller electron capture and $`\beta `$-decay rates in stellar matter (Langanke et al. 1999 and references therein). Following Martinez-Pinedo et al. (2000), this effect was approximately taken into account by decreasing $`\dot{q}_w`$ and $`R_w`$ by a factor of $`5`$. This should be sufficient to account for changes in $`q_{nse}`$ and in the corresponding nuclear energy release $`q_{nse}q_n(t=0)`$ caused by variations of $`Y_e`$. In future simulations, we plan to use $`\dot{q}_w`$ and $`R_w`$ based on the new individual capture rates.

Despite its simplicity, the kinetic scheme reasonably well describes all major stages of carbon burning. In particular, it takes into account the important effect of energy release or absorption caused by changes in NSE composition when matter expands or contracts (these changes happen on a quasi-equilibrium rather than on the equilibrium timescale). When the WD expands, this effect adds $`50`$% to the energy initially released by burning at high densities. Using even the simplest 13-species $`\alpha `$-network to account for this effect would be prohibitively expensive in three-dimensional hydrodynamical simulations.

To couple the kinetic scheme and the front tracking algorithm (Section 2.2), $`q_f`$ in (5) is taken $`q_f=q_{nse}q_n(0)`$ at densities $`\rho >2\times 10^7`$ g/cc, $`q_f=Q_CY_C(t=0)`$ at $`\rho <5\times 10^6`$ g/cc, and is linearly interpolated between the two values in the density range $`5\times 10^62\times 10^7`$ g/cc. The carbon mole fraction $`Y_C`$ inside the front is decreased in proportion to the increase of $`f`$. Thus, at high densities, both carbon consumption and the NSQE relaxation take place inside the flame front. At low densities where the NSQE stage of burning is slow, NSQE relaxation takes place outside the flame front.

2.5. Numerical method

Fluid dynamics equations (1) for $`\rho `$, U, $`E`$, and advection parts of equations (2), (8)-(10) and 912) for chemical variables $`f`$, $`Y_C`$, $`q_n`$, $`\delta _{nse}`$ and $`Y_e`$ are solved using an explicit, second-order, Godunov-type, adaptive mesh refinement algorithm ALLA (Khokhlov, 1998). A Riemann solver for an arbitrary EOS (Colella & Glaz, 1985) is used to evaluate fluxes at cell interfaces. Left and right states are computed using a monotone piecewise-linear reconstruction (van Leer 1979). A small amount of an artificial diffusion flux (Lapidus, 1967) is added to Euler fluxes following Colella and Woodward (1984). This feature is used only in the vicinity of strong shocks. Multidimensions are treated using directional splitting. The diffusion term in (2) is treated using an explicit, second-order finite differencing.

ALLA uses a fully threaded tree (FTT) adaptive mesh refinement to refine the computational mesh on the level of individual cells (Khokhlov 1998; Khokhlov & Chtchelkanova, 1999). FTT is a tree with inverted thread pointers directed from children to neighbors of parent cells. As opposed to conventional tree structures, all operations on the FTT, including tree modifications (refinement and unrefinement) are $`O(N)`$, where N is a number of cells. All of these operations can be performed in parallel. FTT has a $`2\frac{1}{2}`$ integers-per-cell memory overhead, compared with 8 integers for an ordinary tree or 14 integers for a tree that stores connectivity information. Computations using the FTT are virtually as fast (per computational cell) as those on a regular grid. A massively parallel version of the code uses a space-filling curve approach to maintain data locality, achieve load balance, and minimize interprocessor communications. ALLA has been extensively tested and used in various combustion problems involving shocks, flames, turbulence, and their interactions (Gamezo et al., 1999,2000; Khokhlov & Oran, 1999; and references therein) and in astrophysics (Khokhlov et al., 1999).

Mesh refinement was used to refine shocks, contact discontinuities, and regions with large gradients of $`\rho `$, $`P`$, and $`f`$. The expressions for the corresponding shock, $`\xi _s`$, discontinuity, $`\xi _c`$, and gradient refinement indicators, $`\xi _P`$, $`\xi _\rho `$, $`\xi _f`$, are given by the equations (11), (12) and (13) of Khokhlov (1998), respectively. In addition, refinement was introduced for the gradients of tangential velocity using an indicator

$$\xi _u=\underset{j,k=1,\mathrm{},3}{\mathrm{max}}\left((1\delta _{jk})\left(\frac{U_k}{x_j}\right)\right)\frac{\mathrm{\Delta }x}{\alpha a_s}$$
$`(13)`$
with $`\alpha =0.05`$. Refinement was initiated in places where the maximum of all indicators,

$$\xi =\mathrm{max}(\xi _s,\xi _c,\xi _P,\xi _\rho ,\xi _f,\xi _u),$$
$`(14)`$
was larger than the threshold value $`\xi >0.5`$. The mesh was unrefined in places where $`\xi `$ was less than the threshold value $`\xi <0.05`$. The cell size $`\mathrm{\Delta }r`$ is related to the level $`l`$ of the cell in the tree as $`\mathrm{\Delta }r=L/2^l`$. The mesh can be characterized by the minimum and maximum levels of leaves (unsplit cells) in the tree, $`l_{min}`$ and $`l_{max}`$, which were predefined at the beginning of a computation. During a computation, mesh refinement was allowed for cells with $`l_{min}l<l_{max}`$.

Computations were carried out on a 128 $`250`$MHz R10000 SGI Origin 2000 at the Naval Research Laboratory (NRL). On this machine, ALLA advanced $`2\times 10^4`$ cells at the rate of one timestep per second per processor on pure fluid-dynamic, three-dimensional problems with $`\gamma `$ held constant. With the addition of a degenerate matter EOS, flame capturing, and nuclear kinetics, the speed decreased to $`6\times 10^3`$ cell-steps/s/processor.

2.6. Problem setup and initial conditions

Most of the computation were performed for one octant of the WD assuming mirror symmetry along the $`x=0`$, $`y=0`$ and $`z=0`$ planes passing through the WD center. The initial temperature was taken as $`T=10^5`$K everywhere inside the WD. The initial central density was taken as $`\rho _c=2\times 10^9`$ g/cc. Starting from the central pressure $`P(\rho _c)`$, the equations of hydrostatic equilibrium, $`dP/dr=GM\rho /r^2`$ and $`dM/dr=4\pi \rho r^2`$, were integrated outward until $`P=0`$ was reached. The resulting WD configuration was interpolated onto a three-dimensional mesh. Nuclear kinetic variables for the unburned matter and $`f=0`$ were defined everywhere. To initiate burning in the center, $`f=1`$ was then set inside a small sphere of radius $`r_i=3\times 10^6`$ cm and mass $`M_i8\times 10^5M_{WD}`$. For a while, burning was not allowed inside the sphere in order not to impulsively disturb the hydrostatic equilibrium. This procedure smoothly triggered the flame propagation process just outside the sphere. Due to the finite numerical resolution, the shape of the initial flame bubble was not perfectly spherical but contained a number of perturbations. A quadrupole mode was the largest mode of perturbations consistent with the symmetries of the problem. Smaller perturbations were mostly burned out by the flame, but the quadrupole mode partially survived and gave rise to the development of the RT instability.

3. Results

3.1. Hydrostatic equilibrium

Before calculating flame propagation, tests were made to ensure that without burning, the WD remains hydrostatic and keeps its symmetry. Since a Chandrasekhar-mass WD is close to a collapse threshold, its equilibrium is very sensitive to the discretization errors. Thus a rather fine mesh is required near the WD center to maintain equilibrium. The initial mesh was constructed with fine cells at level $`l=l_{max}=10`$ for radii $`r<0.4R_{WD}`$, at $`l=l_{max}1`$ for $`0.4R_{WD}r<0.8R_{WD}`$, at $`l=l_{max}2`$ for $`0.8R_{WD}r<1.1R_{WD}`$, and coarser cells at larger $`r`$. Computations of hydrostatic equilibrium with $`\mathrm{\Delta }r(l_{max})=5\times 10^5`$ cm showed numerical noise $`<3\times 10^5`$cm/s everywhere except in the very outer layers of the WD where the density changed from cell to cell by more than an order of magnitude. Near the center, the noise was less than the fluid velocity generated by the laminar flame, $`S_l(\rho _0\rho _1)/\rho _02\times 10^6`$ cm/s. The total kinetic energy $`E_k`$ of the WD remained $`E_k<10^5E_{th}`$, where $`E_{th}`$ is the total internal energy, during more than two sound crossing times of the WD ($`0.9`$ s of integration). Deviations from spherical symmetry were negligible.

3.2. Deflagration of a $`0.5C+0.5O`$ WD

Deflagration of a CO-WD with initial composition $`X_C=0.5`$ by mass ($`Y_C=1/24`$), central density $`\rho _c=2\times 10^9`$g/cc, and angle-independent small scale burning $`s_1=1`$, $`s_2=0`$ and $`\alpha =1`$ in (7), serves as the baseline case. The computational domain for this simulation had size $`L=5.35\times 10^8`$ cm with levels of refinement $`l_{max}=10`$, $`l_{min}=7`$, and with maximum and minimum resolution $`\mathrm{\Delta }r(l_{max})=5.22\times 10^5`$ and $`\mathrm{\Delta }r(l_{min})=6.18\times 10^6`$ cm, respectively. The simulation was carried out for $`1.8`$s of physical time using Courant number $`0.7`$, required $`9,000`$ timesteps, and took $`6,000`$ CPU hours. The initial fixed mesh described in Section 3.1 was used for approximately $`1`$s. After that, the mesh refinement was turned on. The number of cells used in the simulation grew from $`5,000,000`$ in the beginning to $`30,000,000`$ at the end of the simulation.

Figure 1 presents the evolution with time of the released nuclear energy $`E_n`$, kinetic energy $`E_k`$, binding energy $`E_{tot}`$, and burned mass $`M_b`$ inside the WD. The figure shows that the deflagration resulted in a SNIa explosion. At $`1.8`$s, about $`1.3\times 10^{51}`$ ergs has already been released and $`E_{tot}8\times 10^{50}`$ergs is positive (the WD is unbound). At this moment, the WD radius is $`R_{WD}4.7\times 10^8`$ cm, $`\rho _c7\times 10^7`$g/cc and the expansion velocity $`U_{exp}1.0\times 10^9`$ cm/s. Burning will probably release a few more units of $`10^{50}`$ergs until central density reaches $`10^6`$g/cc and burning quenches completely.

Figure 2 shows the time evolution of the flame surface. Initial perturbations lead to the formation of rising plumes of burned material. In the beginning, the flame surface and the energy release are small, and the flame growth takes place in essentially hydrostatic conditions. At $`t<1`$s, the kinetic energy remains less that $`10^{47}`$ ergs or $`0.01`$% of the binding energy of the WD (Figure 1b). The flame surface increases significantly and becomes more wrinkled at later times. The energy release rate increases, and the WD begins to expand.

Figures 3 and 4 show the flame variable $`f`$ and the distribution of radial velocity $`V_r`$ at $`1.44`$s in three orthogonal planes, $`x=3.2\times 10^7`$cm, $`y=3.2\times 10^7`$cm, and $`z=3.2\times 10^7`$cm, offset from the WD center. At this time the flame has approximately reached the half-radius of the WD. The figures illustrate several important features of burning. Outer parts of the flame consist of big plumes rising through the WD. The absolute velocity of these plumes is $`6\times 10^8`$cm/s or $`a_s`$. This is about two times larger than the overall expansion velocity of matter, $`U_{exp}3\times 10^8`$cm/s. The plume velocity relative to the expanding material, $`0.5a_s`$, is subsonic. Figure 4 also shows regions of unburned matter sinking to the center. There is a complex flow pattern close to the WD center where fresh fuel is burned by the flame and different parts of the flame collide with each other and annihilate. This is similar to turbulent burning in a vertical column in a uniform gravitational field (K95). In that case, the leading part of the flame is also dominated by large scale bubbles whose motions determine the amount of fuel passing through the flame. The fuel was burned inside the flame brush which is dominated by smaller-scale structures. Since the scale of the largest plumes in a column is limited by the column width, the flame is able to reach a steady state on all scales. In a spherical star, however, more space becomes available for the plumes as the flame moves to larger radii, and the rate of burning must increase with time.

Figure 5 shows the effective turbulent flame speed $`S_t`$ which was estimated as the speed of an equivalent spherical flame burning matter with the same rate,

$$S_t\frac{1}{4\pi \rho r_f}\frac{dM_b}{dt},$$
$`(15)`$
where $`r_f`$ and $`\rho _f`$ are the average radius and the density ahead of the flame brush. In the simulation, a small-scale flame speed $`S`$ determined by (7) was $`(12)\times 10^7`$cm/s, depending on the radius, local density, and time. Figures 1 and 5 show that burning rate increased at $`t<1.41.5`$s and noticeably exceeded the value of $`S`$ due to the growth of the flame surface. However, $`S_t`$ reached maximum at $`1.5`$s and started to decline at later times. The WD becomes unbound soon after maximum $`S_t`$ is reached (see Fig. 1).

Figures 6 and 7 show the flame and radial velocity at the end of the simulation, $`1.8`$s. The leading part of the flame is still dominated by large plumes, whereas inner parts are wrinkled and contain more complex structures. The outer edge of the flame is located at $`0.8R_{WD}`$. The absolute velocity of the plumes has increased and reached $`10^9`$cm/s. However, the average expansion velocity of the WD has increased significantly and is also $`U_{exp}10^9`$cm/s. The velocity of the flame bubbles is essentially zero relative to the expanding matter, and they have practically stopped rising. The expansion completely froze the RT-instability on large scales. Generation of new flame surfaces on large scales ceased. Freezing through expansion is the reason for decrease in burning rate at $`t>1.5`$s.

There was an earlier three-dimensional computation of a deflagration explosion in a CO-WD (K95) that used a nonuniform expanding $`185^3`$ grid, assumed instantaneous transition of burned matter to the state of NSE, and the same prescription for $`S`$ as in this paper. This computation was terminated at $`1.7`$s. At $`1.7`$s, the simulation showed plumes of burned material rising at $`5\times 10^8`$cm/s. Approximately $`35`$% of the matter was burned, in agreement with this paper (Note that $`M_b`$ shown in Figure 16 of K95 must be multiplied by $`8`$, as it erroneously shows the amount of burned matter as per one octant of the star). In the earlier calculation, the WD is still slightly bound.

To illustrate mesh refinement, Figure 8 shows the mesh in the X-Y plane $`z=3.2\times 10^7`$cm at $`1.8`$s. Figure 9 shows the mesh, density and velocity field in a small region in the X-Z plane $`y=3.2\times 10^7`$cm with coordinates $`x=01.2\times 10^8`$cm, $`z=09.4\times 10^7`$cm. Despite the complexity of the flame surface and velocity field inside the flame brush, the overall expansion of the WD is remarkably one-dimensional. This is due to the subsonic nature of burning which allows pressure to equilibrate inside the WD. Small deviations of the WD surface from spherical symmetry are only noticeable at late times when the flame gets very close to it.

3.3. Importance of small-scale burning

Figure 10 shows the time evolution of $`E_n`$, $`E_k`$, $`E_{tot}`$, and $`M_b`$ for the simulation of a $`0.5C+0.5O`$ WD with $`S_{sub}=0`$ in (6). All other parameters were the same as in the baseline case presented in Section 3.2. Flame development was similar in the beginning. However, in this case, burning slowed down and then practically stopped because the expansion decreased $`S_l`$ by $`3`$ orders of magnitude. There was no explosion in this case and the WD remained bound. This shows that some amount of burning on unresolved scales must be incorporated. Although $`S10^7`$cm/s in the baseline case was too small to produce an explosion by itself, it is necessary to allow nuclear energy release at a flame front whose area increases due to turbulence. A finite $`S`$ is also necessary to mimic the reduction in flame surface caused by merging and burnout of small-scale structures. This is crucial in order to reproduce numerically a self-similar regime of turbulent burning in which the burning rate becomes independent of the details of flame behavior on small scales (K95). Figure 11 compares the flame surface at $`1.3`$s for both cases in the $`XY`$ plane passing through the center of the WD. This figure illustrates the tendency of the flame to create more surface when $`S`$ is small and to create less surface when $`S`$ increases. Simulations with more resolution and different prescriptions for $`S`$ are required to determine if we are indeed beginning to see the self-similar behavior of the flame.

3.4. Nucleosynthesis

Figure 12 shows the distribution inside the WD of Fe-group, Si-group (intermediate mass elements), and unburned CO matter at $`1.8`$s. The location of Si-group elements is determined using the criterion $`\delta _{nse}<0.95`$ (equation (10)). The figure shows that unburned CO penetrates deep inside the flame brush and is present in significant amounts near the WD center. Incomplete burning and the formation of intermediate mass elements occur at densities less than $`10^7`$g/cc and take place at radii greater than $`0.3R_{WD}`$. Closer to the center, the density is still too high, and burning produces of Fe-group nuclei. The production of S-group nuclei will move inwards as the WD expands further. Pockets of completely unburned CO are likely to survive in central parts. Before burning quenches, one should also expect the formation of large amounts of magnesium in the central parts of the WD. In one-dimensional models magnesium only forms immediately before quenching of incomplete burning in the outer parts of the WD. here, outermost parts of the WD are expected to stay completely unburned since the penetration of the flame into the these layers has stopped (Section 3.2).

Figure 13 shows the distribution at $`1.8`$s of the neutron excess of matter $`\eta =1=2Y_e`$. The maximum neutron excess, $`\eta 0.013`$, is much smaller than $`0.1`$ typically obtained in one-dimensional SNIa simulations (Nomoto et al, 1984; Woosley & Weaver, 1986; Khokhlov, 1991). This is mostly due to a reduction of the theoretical electron capture rates (Section 2.4). In addition, burned matter spends less time near the center due to RT-mixing. Figure 13 shows that the high-$`\eta `$ products are present at radii up to $`0.8R_{WD}`$.

3.5. Sensitivity to initial conditions

Figure 14 shows the results of a numerical experiment in which the location of the initial flame bubble was offset from the WD center by $`\frac{1}{3}`$ of its radius ($`10^6`$cm). The computational domain was doubled, and a simplified flame propagation scheme with constant $`q_f=q_{nse}q_n(0)`$ equal to the instantaneous NSE energy release was used. Equations (8)-(11) were turned off to save time. A $`30`$% deviation of the ignition cite from the center was enough to cause a highly asymmetric flame.

In another numerical experiment the initial carbon mass fraction $`X_C`$ was decreased from $`X_C=0.5`$ (Case B) to $`X_C=0.2`$. This resulted in less energy release and less buoyancy of burned products. Comparison between the $`X_C=0.5`$ and $`X_C=0.2`$ runs (Figure 15) shows a significant delay and slower increase of burning.

4. Discussion and conclusions

This paper described a physical model and an adaptive mesh refinement numerical algorithm for three-dimensional modeling of Type Ia supernovae. It also presented the first results of three-dimensional simulations using this model of a deflagration in a nonrotating, Chandrasekhar-mass, carbon-oxygen white dwarfs.

The simulation of a $`0.5C+0.5O`$ CO-WD with central density $`\rho =2\times 10^9`$g/cc was carried out for $`1.8`$s. The leading edge of the flame was dominated by large-scale, rising plumes of burned matter. Between these plumes, unburned CO was sinking towards the WD center. Fine-structured flame surface developed closer to the center below the plumes. The simulation lasted long enough to demonstrate the effect of freezing of large-scale turbulence due to expansion predicted in (K95). The effective turbulent burning speed $`S_t`$ increased exponentially due to the growth of the RT-instability, reached maximum value $`S_t2\times 10^8`$ cm/s at $`1.5`$s, and then started to decline. The expansion velocity of the WD reached $`10,000`$km/s by the end of the simulation, and the motion of the buoyant plumes relative to unburned matter practically stopped.

The deflagration resulted in a rather healthy explosion. About 60% of the WD was burned by the end of the simulation, and this released $`1.3\times 10^{51}`$ergs of nuclear energy. The total energy, $`E_{tot}8\times 10^{50}`$ergs, was positive and the WD was unbound. By $`1.8`$s, the WD radius was $`4.7\times 10^8`$cm, and burning still continued in its central parts. It should be expected that a few more units of $`10^{50}`$ergs would be released before burning would be quenched completely, so that the total kinetic energy at infinity should be in the range $`(11.3)\times 10^{51}`$ergs. This is $`5080`$% of the kinetic energy required to produce a typical SNIa explosion. The total burned mass in this simulation is expected to be around $`1M`$. Detailed nucleosynthesis postprocessing is required to predict the amount of $`{}_{}{}^{56}Ni`$ synthesized during the explosion. By analogy with one-dimensional simulations, we expect up to $`50`$% of the burned mass, or $`0.5M_{}`$, to be $`{}_{}{}^{56}Ni`$. This is nearly enough to power a SNIa light curve.

There are several features of the computed explosion, however, that lead us to conclude that it is an incomplete model of a SNIa. One distinct feature is the presence of a massive, $`0.4M_{}`$, outer layer of unburned CO surrounding burned matter. The presence of such a layer was a typical feature of all one-dimensional deflagration models with constant $`S_t`$ or with $`S_t`$ decreasing at the end of the explosion (Woosley & Weaver, Khokhlov 1991). The exception is the W7 deflagration model (Nomoto et al., 1976) which created intermediate mass elements at velocities up to $`14,000`$km/s by maintaining a high deflagration speed and simultaneously avoiding rapid expansion of the outer layers. The simulations presented here show the decrease of the burning rate in the outer layers and are not consistent with the W7-like behavior. The presence of the unburned outer layer significantly limits the maximum expansion velocities of the intermediate mass elements, and this presents a difficulty for the spectra of SNIa at maximum and before maximum light. The observations and spectral analyses indicate that the line-forming region for the intermediate mass elements extends up to $`15,000`$ and in some cases up to $`30,000`$km/s in velocity space (Branch, 1981; Benetti et al. ,1991; Wells et al., 1994; Höflich, 1995; Höflich & Khokhlov, 1996; Lentz et al. 2000).

The simulations also show that the formation of intermediate mass elements in a three-dimensional deflagration explosion occurs throughout almost the entire WD, and that some intermediate mass elements may form very close to the WD center. This presents another difficulty for spectral modeling. The spectra of normal bright SNIa show that the minimum velocity of the line-forming region of $`Si`$ is $`8,00011,000`$km/s (Höflich & Khokhlov, 1996), $`Ca`$, $`4,0006,000`$km/s (Fisher et al., 1997), $`Mg`$, $`13,00014,000`$km/s (Meikle et al., 1996; Höflich, 1997; Wheeler et al., 1998; Gerardy et al., 2000). There is also an indication that the minimum velocity of carbon in SN1990N (CII at 6580 A) is larger than 26,000 km/s (Fisher et al. 1997), which means that practically all of the material in this supernova has been burned. According to the simulations, unburned CO is also likely to remain in the central regions of the deflagration supernova where it will coexist with $`{}_{}{}^{56}Ni`$. Whether this can lead to the excitation of oxygen lines in the late-time spectra when the positron energy deposition becomes nonlocal (Milne et al. 1999)is an open question which requires further investigation (Höflich, private communication).

The formation of a massive outer layer of unburned material is a large-scale effect caused by the competition between the of buoyancy of large-scales structures (rising plumes) and the global expansion of the WD. In this competition, buoyancy eventually loses. During the explosion, the expansion velocity of the star increases gradually, and at $`1.8`$, outer layers have velocities comparable or even exceeding the velocity of the large-scale plumes of burned products. Even if burning on small scales is underestimated in the simulation, increasing the energy released by small scales would increase the expansion velocity, but would hardly affect the buoyancy of the plumes. As the WD continues to expand into the vacuum, the rarefaction will accelerate the outer layers and this will further increase their velocity relative to the burned material.

The difficulties associated with the deflagration explosion can be avoided if burning turns into a detonation during the explosion (delayed detonation). Detonation would easily overcome the expansion and would incinerate the outer unburned layers of CO. This was one of the initial reasons for the introduction of the delayed detonation model of SNIa. It is also clear that the delayed detonation will help with the potential problems associated with three-dimensional effects by wiping out composition inhomogeneities in the inner layers.

These are the first simulations with the model presented in the paper. As such, they show new effects but also rais many questions. Computations with more resolution and with different subgrid burning models are required to determine if the resolved scales are small enough that burning on these scales is in a self-similar regime. In particular, at maximum burning ($`1.5`$s), the simulations show large-scale motions with relative velocities of $`5\times 10^8`$cm/s. Thus there should be shear instabilities and a strong turbulent cascade on the sides of rising plumes. This should be incorporated into the model. The question of the deflagration-to-detonation transition (DDT) must also be addressed. A rough estimate can be made as follows. Assume that at $`1.5`$s, the turbulent velocity is $`10^8`$cm/s on the driving scale $`L10^8`$cm. Then, from Kolmogorov cascade, the turbulent velocity $`U_\lambda 10^5\lambda `$ cm/s on a scale $`\lambda `$, where $`\lambda `$ is in centimeters. At $`\rho 10^7`$g/cc, the thickness of the laminar flame is $`\delta _l10^2`$cm, and on the scale $`\lambda \delta _l`$ the velocity $`U_\lambda 10^7`$cm is much greater than $`S_l`$. Thus, the turbulent velocity on this scale may be strong enough to move the flame from the flamelet into the distributed regime, create hot spots of nonuniformly preheated unburned material, and thus trigger a detonation (Khokhlov et al., 1997). To trigger a detonation, the size of a nonuniform region must exceed a critical value which depends on density. At $`\rho 10^7`$g/cc, the critical size $`10^6`$cm is much less than $`R_{WD}`$ and the size of large scale structures obtained in the simulations. Thus the transition to detonation seems to be possible. Significant work requires to test these ideas. This paper shows that high-resolution three-dimensional modeling of SNIa is now a reality.

The calculations presented here indicate that the result of the explosion may be sensitive to varying the position of the ignition site and the initial carbon fraction $`X_C`$ in a WD. Decreasing of $`X_C`$ leads to less energy release and less buoyant products of burning. This results in much slower development of the RT-instability, so that the freezing due to expansion may be more effective. On the other hand, smaller energy release leads to slower expansion as well. It is most likely that decreasing $`X_C`$ would result in a weaker explosion. Sensitivity to the position of the ignition site indicates that the results may also be sensitive to stellar rotation, which may lead to a preferential growth or suppression of the initial perturbations along the rotational axis. This makes both $`X_C`$ and rotation potentially good candidates for being the ‘hidden parameters’ to describe the observed diversity of SNIa (Phillips, 1993).

Acknowledgments. I am grateful to Rob Rosenberg (NRL) and Almadena Chtchelkanova (Berkeley Research Associates) who developed FTT-based tools for the surface rendering and visualization used in this work. This paper benefited greatly from discussions with Peter Höflich (UT Austin), Elaine Oran (NRL), Craig Wheeler (UT Austin), Vadim Gamezo (Berkeley Research Associates), and Stirling Colgate (LANL). Work was supported in part by the Naval Research Laboratory through the Office of Naval Research and by the NASA ATP program.

References

Arnett, W.D. & Livne, E. 1994. ApJ 427, 314.

Benetti S., Cappellaro E., Turatto M. 1991, AA 247 410

Branch, D. 1981. ApJ, 248, 1076.

Colella, P. & Glaz, H.M. 1985. J. Comput. Phys., 59, 264.

Colella, P. & Woodward, P.R. 1984. J. Comput. Phys., 54, 174.

Fisher A., Branch D., Höflich P., Khokhlov A. 1995. ApJ 447L, 73.

Fisher A., Branch D., Nugent P., Baron E. 1997, ApJ 481L, 89.

Fowler, W.A., Caughlan, G.R. & Zimmerman, B.A. 1975. ARA&A, 13, 69.

Gamezo, V.N., Khokhlov, A.M. & Oran, S.E., Proceedings of the 17th International Colloquium on the Dynamics of Explosions and Reactive Systems, Universität Heidelberg, IWR, 1999.

Gamezo, V.N., Vasil’ev, A.A., Khokhlov, A.M. & Oran, S.E., Proceedings of the 28th International Symposium on Combustion, vol.28, 2000 (in press).

Gerardy C., Höflich P., Fesen R., Wheeler J.C. 2000, AJ, to be submitted

Hillebrandt, W. & Niemeyer, J.C 2000. astro-ph/0006305

Höflich, P. 1995. ApJ 443, 89.

Höflich P. 1997, in: Second Oak Ridge Symposium on Atomic & Nuclear Astrophysics, ed. A. Mezzacappa, IOP Publishing, p. 693

Höflich & Khokhlov 1996, ApJ 457, 500

Khokhlov, A.M. 1991. A&A, 245, 114.

Khokhlov, A.M. 1995. ApJ, 449, 695.

Khokhlov, A.M., Oran, E.S., Wheeler, J.C. 1997. ApJ, 478, 678.

Khokhlov, A.M. 1998. J. Comput. Phys., 143, 519.

Khokhlov, A.M. & Chtchelkanova, A.Yu. 1999. Proceedings of the 9-th SIAM Conference on Parallel Processing for Scientific Computing, San Antonio, TX, March 22-24, 1999.

Khokhlov, A.M., Höflich, P.A., Oran, E.S., Wheeler, J.C., Wang, L. & Chtchelkanova, A.Yu. 1999. ApJ, 524, L107.

Khokhlov, A.M. & Oran, S.E. 1999. Combustion & Flame, 119, 400.

Langanke, K. & Martinez-Pinedo, G. 1999. Phys. Lett. B., 453, 187.

Lapidus, A. 1967. J. Comput. Phys. 2, 154.

Lentz, E.J., Baron, E., Branch, D., Hauschieldt, P.H. 2000. ApJ, submitted; astro-ph/0007320.

Livne, E. 1993. ApJ 415, L17.

Martinez-Pinedo, G., Langanke, K. & Dean, D.J. 2000. ApJ, 126:493.

Meikle W.P.S., Cumming, , R.J., Geballe, T.R., et al. 1996. MNRAS 281, 263.

Milne P. A., The, L.-S.; Leising, M. D. 1999, ApJS 124, 503

Nomoto, K., Sugimoto, D. & Neo, S. 1976. AP&SS, 39, L37.

Nomoto, K., Thielemann, F.-K., Yokoi, K. 1984. ApJ, 286, 644.

Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999. ApJ, 517:565.

Phillips, M.M. 1993. ApJ, 413, L105.

Pskovskii, I. p. 1977. Soviet Astronomy, 21, 675.

Reinecke, M., Hillebrandt, W. & Niemeyer, J.C. 2000. A&A, in press.

Riess, A.G., Filippenko, A.V., Lui, M.C., et al. 2000. ApJ, 563:62.

Timmes, F. & Woosley, S.E. 1992. ApJ, 396, 649.

Wells L.A., Phillips M.M., Suntzeff B. et al. 1994. AJ 108, 2233.

Wheeler J.C., Höflich P., Harkness R., Spyromillo J. 1998, ApJ 496, 908

Woosley, S.E. & Weaver, T. 1986. ARA&A, 24, 205.

Figure captions

Figure 1a. - Figure 1a. - Total ($`E_{tot}`$), kinetic ($`E_k`$), released nuclear ($`E_n`$) energy and burned mass $`M_b`$ of a $`\frac{1}{2}C+\frac{1}{2}O`$ WD during a deflagration explosion (Section 3.2).

Figure 1b. - Solid line - kinetic energy ($`E_k`$) of a $`\frac{1}{2}C+\frac{1}{2}O`$ WD during a deflagration explosion at earlier times (Section 3.2). Dashed line - kinetic energy with burning turned off (Section 3.1).

Figure 2a. - Time evolution of the flame surface during a deflagration explosion (Section 3.2).

Figure 2b. - Same as Fig. 2a. Details of the flame structure at $`1.8`$s.

Figure 3. - Flame variable $`f`$ (Equation (2)) for a deflagration explosion Section 3.2 at $`1.44`$s shown in three orthogonal planes offset by $`3.2\times 10^7`$cm from the WD center.

Figure 4. - Same as Figure 3 but shows the radial velocity.

Figure 5. - Estimated turbulent flame speed $`S_t`$ (Equation (15)) during a deflagration explosion Section 3.2.

Figure 6. - Flame variable $`f`$ (Equation (2)) for a deflagration explosion Section 3.2 at the end of the simulation $`1.8`$s shown in three orthogonal planes offset by $`3.2\times 10^7`$cm from the WD center.

Figure 7. - Same as Figure 6 but shows the radial velocity.

Figure 8. - Mesh and flame variable $`f`$ in the X-Y plane passing through the WD center for a deflagration explosion Section 3.2 at the end of the simulation $`t=1.8`$s.

Figure 9a. - Density and velocity in the X-Z plane passing through the WD center for a deflagration explosion Section 3.2 at the end of the simulation $`1.8`$s. The region $`x=1.2\times 10^8`$cm and $`z=09.4\times 10^7`$cm is shown.

Figure 9b. - Mesh in the X-Z plane passing through the WD center for a deflagration explosion Section 3.2 at the end of the simulation $`1.8`$s. The region $`x=1.2\times 10^8`$cm and $`z=09.4\times 10^7`$cm is shown.

Figure 10. - Time evolution of the total ($`E_{tot}`$), kinetic ($`E_k`$), released nuclear ($`E_n`$) energy and burned mass $`M_b`$ of a $`\frac{1}{2}C+\frac{1}{2}O`$ WD during a deflagration explosion with the laminar flame speed $`S=S_l`$ (Section 3.3).

Figure 11. - Flame surface in the X-Y plane at $`1.3`$s for the simulation with (right panel, Section 3.2) and without (left panel, Section 3.3) turbulent small-scale burning.

Figure 12. - Distribution of the Fe-group, Si-group, and CO matter for a deflagration explosion Section 3.2 at the end of the simulation $`t=1.8`$s shown in three orthogonal planes offset by $`3.2\times 10^7`$cm from the WD center.

Figure 13. - Same as Figure 12 but shows the neutron excess $`\eta =12Y_e`$.

Figure 14. - Deflagration starting with the offset ignition point (Section 3.5). Shows the flame variable and velocity at $`1`$s. Black square shows the ignition point.

Figure 15. - Evolution of the kinetic energy during deflagration explosion of a WD with the initial $`0.5C+0.5O`$ and $`0.2C+0.8O`$ compositions (Section 3.5).


format: text

text:
# Non–standard Construction of Hamiltonian Structures and of the Hamilton–Jacobi equation.

## I Introduction

The standard construction of Hamiltonian theories starts from a Lagrangian. The procedure is well known and it is the subject of very many textbooks. Nevertheless it is interesting to construct hamiltonian structures for classical systems of differential equations, without using a Lagrangian, which may even fail to exist, starting from the equation of motion only. These approaches usually are not systematic, conducting to non unique hamiltonian structures and deal with non canonical variables (in some systems these variables are more natural than canonical ones).

One systematic, but complicated procedure to construct hamiltonian structures using symmetries of classical equations was intoduced by E. Sudarshan and N. Mukunda in. They also comment the possible relevance of alternative hamiltonian of classical systems under the quantization procedure. Another systematic approach due also to Sudarshan and Mukunda , but developed by J. Marsden, A. Weistein, D. Holm and B. Kupershmidt in  is based on definition of Lie – Poisson structures on duals of Lie algebras. More recently S. Hojman, et. al  presented a newly devised method which constitutes a general technique for construction of hamiltonian structures using symmetries and constants of motion of dynamical systems.

Non standard introduction of hamiltonian structures enrich the possibilities we have in the analysis of both classical and quantum systems and for this reason one must be careful interpreting the results we obtained. In the section III of the present work we use this fact in the context of analysis of stability for dynamical systems in the vicinity of a critical point, demostrating that the criterion of stability used by Bridges in  is inadequate.

On the other hand the Hamilton – Jacobi equation (H–J) is one of the cornerstones of the theoretical physics. Its role in the theoretical setting of optics and both classical and quantum mechanic is well established. Its use as a tool for solving problems is also well known.

The usual way of deriving the H–J equation is based on the lagrangian or hamiltonian formalisms. In the section IV of this paper we present one way of constructing the H–J equation starting from equation of motion only, to do this we use the approaches due to S. Hojman et. al. for non standard construction of hamiltonian structures. In some sense this complement the analysis carry out in the section III, because we only concern here with the study of regular points of dynamical systems and not concerns singular points of it, as in the above mentioned section.

With the help of alternative hamiltonian structures introduced by Hojman’s recipe, we define the Poisson brackets on the phase space M for dynamical system. Using a Poisson map $`\phi `$, M is reduced to a symplectic submanifold N. On N we can define the $`action`$ and the respective H–J equation.

The hamiltonian structures on M are so on N. Furthermore the solutions of dynamical system with initial conditions on N stay all the time on it. For these reasons the solution of H–J equation is a particular solution of dynamical system constrained to N.

There it is not unique way to define hamiltonian structure on M, neither the Poisson map used for reduction of the phase space to submanifold N is unique. This means that we can define the action in different ways and obtain different H–J equations. In this paper we show that solutions obtained from different H–J equations, from alternative hamiltonian structures on M yield correct solutions of dynamical systems. This allow us to speak about some kind of equivalence between different H–J equations obtained this way.

In the section II, for illustration purpose of using the method introduced by S. Hojman, we present one example borrowed from D. Armbruster, J. Guckenheimer, and P. Holmes in . This example was employed in the analysis of normal forms with O(2) symmetry. Using the complete set of symmetries of this dynamical system, it is introduced two alternative hamiltonian structures. With the help of O(2) symmetry we can reduce the four equations dynamical system in cartesian coordinates to three equations in polar coordinates and the dynamical system obtained this way we will use as one of the example presented in the section II. The physics motivations for the study of these systems can be found in the above listed references.

The section III is dedicated to demonstrated that the criterion of Bridges in , that relates the slope of a Casimir function parametrized by the Lagrange multiplier to critical point type depends of the used hamiltonian structures and it is inadequate for this reason.

The section IV is devoted to illustrated the construction of H–J in the way described above. We employ the examples used in the section III and present a complete analysis of the Euler top.

## II One example of non–standard construction of hamiltonian structure.

The contents of this Section illustrates partially the results of .

Consider a dynamical system for the complex variables $`Z_1`$ and $`Z_2`$ defined by the following equations:

$`\dot{Z}_1`$ $`=`$ $`\overline{Z}_1Z_2`$ (1)
$`\dot{Z}_2`$ $`=`$ $`Z_1^2`$ (3)
The system (LABEL:equ) can be writing in cartesian coordinates:

$`\dot{X}_1`$ $`=`$ $`X_1X_2+Y_1Y_2`$ (4)
$`\dot{Y}_1`$ $`=`$ $`X_1Y_2Y_1X_2`$ (5)
$`\dot{X}_2`$ $`=`$ $`X_1^2+Y_1^2`$ (6)
$`\dot{Y}_2`$ $`=`$ $`2X_1Y_1`$ (7)
Where:

$`Z_1`$ $`=`$ $`X_1+iY_1`$ (8)
$`Z_2`$ $`=`$ $`X_2+iY_2`$ (9)
A hamiltonian structure for (LABEL:equ or 6) consist of an antisymmetric matrix, $`J^{ab}`$ and a Hamiltonian $`H`$ such that $`J^{ab}`$ is the Poisson bracket for the dynamical variables, which are in general non–canonical. In addition to its antisymmetry, the matrix $`J^{ab}`$ is required to satisfy Jacobi identity and to reproduce, in conjunction with the Hamiltonian $`H`$ the dynamical equation (LABEL:equ or 6) i. e.,

$$J_{}^{ab}{}_{,d}{}^{}J^{dc}+J_{}^{bc}{}_{,d}{}^{}J^{da}+J_{}^{ca}{}_{,d}{}^{}J^{db}0$$
(10)
and,

$$J^{ab}\frac{H}{x^b}=f^a.$$
(11)
Where:

$$f=\left(\begin{array}{cc}X_1X_2+Y_1Y_2& \\ X_1Y_2Y_1X_2& \\ X_1^2+Y_1^2& \\ 2X_1Y_1& \end{array}\right)$$
(12)
It has been proved  that one solution to the problem of finding a Hamiltonian structure for a given dynamical system is provided by one constant of the motion which may be used as the Hamiltonian $`H`$, and a symmetry vector $`\eta ^a`$ which allows for the construction of a Poisson matrix $`J^{ab}`$. The constant of the motion and the symmetry vector satisfy,

$`_fH`$ $`=`$ $`0,`$ (13)
$`(_t+_f)\eta ^a`$ $`=`$ $`0,`$ (14)
respectively, where $`_f`$ is the Lie derivative along $`f`$ (for a definition, see , for instance). In addition, it is required that the deformation $`K`$ of $`H`$ along $`\eta ^a`$,

$$K\frac{H}{x^a}\eta ^a=_\eta H,$$
(15)
be non-vanishing.

The Poisson matrix $`J^{ab}`$ is constructed as the antisymmetrized product of the flow vector $`f^a`$ and the “normalized” symmetry vector $`\eta ^b/K`$,

$$J^{ab}=\frac{1}{K}(f^a\eta ^bf^b\eta ^a).$$
(16)
The Poisson matrix so constructed has rank 2 and it is, therefore, singular. Adding together two Poisson matrices constructed according to (16) will not increase its rank. It will just redefine the symmetry vector used to construct it. One method to increase the rank of such a Poisson matrix is presented in .

The system (LABEL:equ) is invariant under the transformation:

$`Z_1`$ $``$ $`Z_1\mathrm{exp}(i\mathrm{\Phi })`$ (17)
$`Z_2`$ $``$ $`Z_2\mathrm{exp}(2i\mathrm{\Phi })`$ (19)
In cartesian coordinates the transformation (LABEL:trf) is writing as:

$`\left(\begin{array}{cc}\stackrel{~}{𝒳}_1& \\ \stackrel{~}{𝒳}_2& \end{array}\right)=\left(\begin{array}{cc}\stackrel{~}{X}_1& \\ \stackrel{~}{Y}_1& \\ \stackrel{~}{X}_2& \\ \stackrel{~}{Y}_2& \end{array}\right)=\left(\begin{array}{cc}(\mathrm{\Phi })& 0\\ 0& (2\mathrm{\Phi })\end{array}\right)`$

Where:

$`(\mathrm{\Phi })=\left(\begin{array}{cc}\mathrm{cos}\mathrm{\Phi }& \mathrm{sin}\mathrm{\Phi }\\ \mathrm{sin}\mathrm{\Phi }& \mathrm{cos}\mathrm{\Phi }\end{array}\right)`$

From this, we have the symmetry vector:

$$\eta ^{(1)}=\left(\begin{array}{cc}Y_1& \\ X_1& \\ 2Y_2& \\ 2X_2& \end{array}\right)$$
(20)
Another symmetry transformation, involving space coordinates and time of (6) is:

$`\left(\begin{array}{cc}\stackrel{~}{𝒳}_1& \\ \stackrel{~}{𝒳}_2& \\ \stackrel{~}{t}& \end{array}\right)=\left(\begin{array}{cc}𝒳_1& \\ 𝒳_2& \\ \lambda ^1t& \end{array}\right)`$

The infinitesimal version of this transformation is:

$$\stackrel{~}{\eta }^{(2)}=\left(\begin{array}{cc}𝒳_1& \\ 𝒳_2& \\ t& \end{array}\right)$$
(21)
This vector we can equivalently write in the form: (see)

$$\eta ^{(2)}=\left(\begin{array}{cc}X_1+t(X_1X_2+Y_1Y_2)& \\ Y_1+t(X_1Y_2Y_1X_2)& \\ X_2+t(X_1^2+Y_1^2)& \\ Y_2t(2X_1Y_1)& \end{array}\right)$$
(22)
The symmetry vectors $`\eta ^{(1)}`$, $`\eta ^{(2)}`$ and $`f`$ satisfy the symmetry equation (14) and the Lie algebra:

$`[\eta ^{(2)},f]`$ $`=`$ $`f`$ (23)
$`[\eta ^{(1)},\eta ^{(2)}]`$ $`=`$ $`0`$ (24)
$`[\eta ^{(1)},f]`$ $`=`$ $`0`$ (25)
For the system (LABEL:equ) we have the following constants of motion:

$`C_1`$ $`=`$ $`|Z_1|^2+|Z_2|^2`$ (26)
$`C_2`$ $`=`$ $`{\displaystyle \frac{1}{2i}}(Z_1^2\overline{Z}_2\overline{Z}_1^2Z_2)`$ (27)
Or in cartesian coordinates:

$`C_1`$ $`=`$ $`X_1^2+Y_1^2+X^2\mathrm{\_}2+Y_2^2`$ (28)
$`C_2`$ $`=`$ $`2X_1Y_1X_2Y_2(X_1^2Y_1^2)`$ (29)
We can now construct some alternative Poisson structures taking as Hamiltonian $`H`$ = $`C_1`$ and the Poisson matrix:

$$J^{ab}=\frac{1}{K}(f^a\eta ^{(2)b}f^b\eta ^{(2)a}).$$
(30)
Where $`K`$ is the deformation of $`H`$ = $`C_1`$ along $`\eta ^{(2)}`$,

$$K\frac{H}{x^a}\eta ^a=_\eta H=2C_1,$$
(31)
Or taking $`C_2`$ as a Hamiltonian $`H`$, and Poisson matrix:

$$J^{ab}=\frac{1}{Q}(f^a\eta ^{(2)b}f^b\eta ^{(2)a}).$$
(32)
Where $`Q`$ is the deformation of $`H`$ = $`C_2`$ along $`\eta ^{(2)}`$,

$$Q\frac{H}{x^a}\eta ^a=_\eta H=3C_2,$$
(33)
The deformations of both $`C_1`$ and $`C_2`$ along $`\eta ^{(1)}`$ are vanishing and for this reason we can not use this vector in the construction of Poisson matrix with the above procedure.

Finally we note that introducing polar coordinates and using the O(2) symmetry, the system (LABEL:equ) can be writing as:

$`\dot{r}_1`$ $`=`$ $`r_1r_2\mathrm{cos}\mathrm{\Theta }`$ (34)
$`\dot{r}_2`$ $`=`$ $`\pm r_1^2\mathrm{cos}\mathrm{\Theta }`$ (35)
$`\dot{\mathrm{\Theta }}`$ $`=`$ $`\left(2r_2\pm {\displaystyle \frac{r_1^2}{r_2}}\right)\mathrm{sin}\mathrm{\Theta }`$ (36)
Where: $`Z_j=r_j\mathrm{exp}\mathrm{\Phi }_j`$, $`\mathrm{\Theta }=2\mathrm{\Phi }_1\mathrm{\Phi }_2`$

## III Hamiltonian structures and linear stability criterion.

In this section we use another way to construct Poisson structures based on the knowledge of constants of motion for the dynamical system under consideration (See  for detail) to demonstrate that stability criterion used by Bridges in  strongly depends of the used Poisson structures. Let us consider as the first example the one we presented in the preceding section.

Example I:

$`\dot{X}_1`$ $`=`$ $`X_1X_2\mathrm{cos}\mathrm{\Phi }`$ (37)
$`\dot{X}_2`$ $`=`$ $`\pm X_1^2\mathrm{cos}\mathrm{\Phi }`$ (38)
$`\dot{\mathrm{\Phi }}`$ $`=`$ $`\left(2X_2\pm {\displaystyle \frac{X_1^2}{X_2}}\right)\mathrm{sin}\mathrm{\Theta }`$ (39)
The phase space of (38) is:

$`M=\{\mathrm{\Phi },X_1,X_2:(X_1,X_2)^{2+},0\mathrm{\Phi }2\pi \}`$

The integrals of motion for this system are:

$`C_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}(X_1^2X_1^2)`$ (40)
$`C_2`$ $`=`$ $`X_1^2X_2\mathrm{sin}\mathrm{\Phi }`$ (41)
Taking as a Hamiltonian $`H=C_1`$ and the function $`\mathrm{\Psi }=C_2`$ we have the following Poisson structure:

I Poisson Structure:

$`J^{(1)ab}={\displaystyle \frac{1}{X_1X_2}}\epsilon ^{abc}{\displaystyle \frac{\mathrm{\Psi }}{X^c}}`$

$$J^{(1)ab}=\left(\begin{array}{ccc}0\hfill & X_1\mathrm{cos}\mathrm{\Phi }& \hfill \frac{X_1}{X_2}\mathrm{sin}\mathrm{\Phi }\\ X_1\mathrm{cos}\mathrm{\Phi }\hfill & 0& \hfill 2\mathrm{sin}\mathrm{\Phi }\\ \frac{X_1}{X_2}\mathrm{sin}\mathrm{\Phi }\hfill & 2\mathrm{sin}\mathrm{\Phi }& \hfill 0\end{array}\right)$$
(42)
It is not difficult to see that (42) together with $`H=C_1`$ provide a Hamiltonian formulation for (38), and that $`\mathrm{\Psi }`$ is the Casimir function ( i. e. function which has vanishing Poisson bracket relations with any other dynamical quantity. For detail, see  for instance) for this Poisson structure.

In other words we see that,

$`\dot{𝒳}=\{𝒳,H\}_M`$

Where: $`𝒳=(X_1,X_2,\mathrm{\Phi })`$

Let $`\lambda `$ be a Lagrange multiplier, then a necessary condition for $`𝒳`$ to be a critical point of (38) given a Hamiltonian $`H`$ and a Casimir $`\mathrm{\Psi }`$ is that:

$`H(𝒳)\lambda \mathrm{\Psi }(𝒳)=0`$

Or

$`\left(\begin{array}{cc}X_1& \\ X_2& \\ 0& \end{array}\right)\lambda \left(\begin{array}{cc}2X_1X_2\mathrm{sin}\mathrm{\Phi }& \\ X_1^2\mathrm{sin}\mathrm{\Phi }& \\ X_1^2X_2\mathrm{cos}\mathrm{\Phi }& \end{array}\right)=0`$

From this equations we have the critical point:

$$P_c=(\frac{1}{\sqrt{2}\lambda },\frac{1}{2\lambda },\frac{\pi }{2})$$
(43)
On the other hand, the derivative of $`\mathrm{\Psi }_{P_c}`$ and $`H_{P_c}`$ with respect to Lagrange multiplier are:

$`{\displaystyle \frac{d\mathrm{\Psi }_{P_c}}{d\lambda }}`$ $`=`$ $`{\displaystyle \frac{3}{4\lambda ^4}}`$ (44)
$`{\displaystyle \frac{dH_{P_c}}{d\lambda }}`$ $`=`$ $`{\displaystyle \frac{3}{4\lambda ^3}}`$ (45)
Liniarizing the Poisson system (38) about a critical point $`P_c`$: $`𝒳^i=P_c+\xi ^i`$, satisfying $`H=\lambda \mathrm{\Psi }`$ we have:

$$\dot{\xi }^i=\frac{}{X^m}\left[\frac{1}{X_1X_2}(\times H)^i\right]_{Pc}\xi ^m=\left[\frac{1}{X_1X_2}\mathrm{\Psi }\times (D^2H\lambda D^2\mathrm{\Psi })_{P_c}\right]_m^i\xi ^m$$
(46)
Where:

$$\frac{1}{X_1X_2}\mathrm{\Psi }\times (D^2H\lambda D^2\mathrm{\Psi })_{P_c}=\left(\begin{array}{ccc}0\hfill & 0& \hfill \frac{1}{2\sqrt{2}\lambda ^2}\\ 0\hfill & 0& \hfill \frac{1}{2\lambda ^2}\\ 2\sqrt{2}\hfill & 4& \hfill 0\end{array}\right)$$
(47)
One eigenvalue of (47) is zero and the other two are given by:

$$\mu ^2=\frac{3}{\lambda ^2}=4\lambda ^2\frac{d\mathrm{\Psi }_{P_c}}{d\lambda }$$
(48)
Where we have used (44).

This prove that for the I Poisson structure of dynamical system (38) the critical point is elliptic (purely imaginary eigenvalues) when $`d\mathrm{\Psi }/d\lambda <0`$ and hyperbolic when $`d\mathrm{\Psi }/d\lambda >0`$. This statement is agree with the stability criterion used by Bridges in . But this criterion depends of the used Poisson structure as we will see, taking another Poisson structure the criterion change.

Another Poisson structure for (38) is obtained taking $`C_2`$ as a Hamiltonian $`H`$ and $`C_1`$ as a Casimir $`\mathrm{\Psi }`$.

II Poisson Structure:

$`J^{(2)ab}={\displaystyle \frac{1}{X_1X_2}}\epsilon ^{abc}{\displaystyle \frac{\mathrm{\Psi }}{X^c}}`$

$$J^{(2)ab}=\left(\begin{array}{ccc}0\hfill & 0& \hfill \frac{1}{X_1}\\ 0\hfill & 0& \hfill \pm \frac{1}{X_2}\\ \frac{1}{X_1}\hfill & \frac{1}{X_2}& \hfill 0\end{array}\right)$$
(49)
The critical point of (38) using the second Poisson structure is obtained from the equation:

$`H(𝒳)\beta \mathrm{\Psi }(𝒳)=0`$

From this, we have:

$$P_c=(\frac{\beta }{\sqrt{2}},\frac{\beta }{2},\frac{\pi }{2})$$
(50)
Note that directly from the equation of motion (38) the critical point must satisfy the relation: $`X_1^2=2X_2^2,\mathrm{\Phi }=\frac{\pi }{2}`$.

The derivative of $`\mathrm{\Psi }_{P_c}`$ and $`H_{P_c}`$ with respect to Lagrange multiplier are:

$`{\displaystyle \frac{d\mathrm{\Psi }_{P_c}}{d\beta }}`$ $`=`$ $`{\displaystyle \frac{3}{4}}\beta `$ (51)
$`{\displaystyle \frac{dH_{P_c}}{d\beta }}`$ $`=`$ $`{\displaystyle \frac{3}{4}}\beta ^2`$ (52)
Liniarizing the Poisson system (38) about a critical point $`P_c`$ we have the matrix:

$$\frac{1}{X_1X_2}\mathrm{\Psi }\times (D^2H\beta D^2\mathrm{\Psi })_{P_c}=\left(\begin{array}{ccc}0\hfill & 0& \hfill \frac{\beta ^2}{2\sqrt{2}}\\ 0\hfill & 0& \hfill \frac{\beta ^2}{2}\\ 2\sqrt{2}\hfill & 4& \hfill 0\end{array}\right)$$
(53)
The eigenvalues of (53) are zero and:

$$\mu ^2=3\beta ^2=4\beta \frac{d\mathrm{\Psi }_{P_c}}{d\beta }$$
(54)
Where we have used (51). This prove that for the II Poisson structure of dynamical system (38) the critical point is elliptic (purely imaginary eigenvalues) when $`d\mathrm{\Psi }/d\beta >0`$ and hyperbolic when $`d\mathrm{\Psi }/d\beta <0`$. This statement is completly opposite to the mentioned above stability criterion used by Bridges in . This fact confirms that this criterion depends of the used Poisson structure.

As the second example illustrating this fact we use the system presented by T. Bridges in .

Example II:

$`\dot{X}_1`$ $`=`$ $`2X_3`$ (55)
$`\dot{X}_2`$ $`=`$ $`2Q(X_1)X_3`$ (56)
$`\dot{X}_3`$ $`=`$ $`X_2+Q(X_1)X_1`$ (57)
Where $`Q(X_1)`$ is a real polynomial in $`X_1`$. The phase space of (56) is: $`M=^{3+}`$

The integrals of motion for this system are:

$`C_1`$ $`=`$ $`X_3^2X_1X_2`$ (58)
$`C_2`$ $`=`$ $`X_2{\displaystyle _0^{X_1}}Q(s)𝑑s`$ (59)
I Poisson Structure $`H=C_1`$, and $`\mathrm{\Psi }=C_2`$.

$`J^{(1)ab}=\epsilon ^{abc}{\displaystyle \frac{\mathrm{\Psi }}{X^c}}`$

$$J^{(1)ab}=\left(\begin{array}{ccc}0\hfill & 0& \hfill 1\\ 0\hfill & 0& \hfill Q(X_1)\\ 1\hfill & Q(X_1)& \hfill 0\end{array}\right)$$
(60)
With this Poisson structure the equation of motion (56) take the simple form: $`\dot{\stackrel{}{X}}=\mathrm{\Psi }\times H`$. When $`\mathrm{\Psi }`$ and $`H`$ are parallel the system (56) has a critical point.

$`H(𝒳)+\lambda \mathrm{\Psi }(𝒳)=0`$

From this, we have:

$$P_c=(\lambda ,\lambda Q(\lambda ),0)$$
(61)
The derivative of $`\mathrm{\Psi }_{P_c}`$ and $`H_{P_c}`$ with respect to Lagrange multiplier are:

$`{\displaystyle \frac{dH_{P_c}}{d\lambda }}`$ $`=`$ $`2\lambda Q(\lambda )+\lambda ^2{\displaystyle \frac{dQ}{d\lambda }}`$ (62)
$`{\displaystyle \frac{d\mathrm{\Phi }_{P_c}}{d\lambda }}`$ $`=`$ $`\left[Q(\lambda )+\lambda {\displaystyle \frac{dQ}{d\lambda }}\right]`$ (63)
Liniarizing the Poisson system (56) about a critical point $`P_c`$: $`X^i=P_c+\xi ^i`$, we have:

$$\dot{\stackrel{}{\xi }}=[\mathrm{\Psi }\times (D^2H+\lambda D^2\mathrm{\Psi })_{P_c}]\stackrel{}{\xi }$$
(64)
Where:

$$\mathrm{\Psi }\times (D^2H\beta D^2\mathrm{\Psi })_{P_c}=\left(\begin{array}{ccc}0\hfill & 0& \hfill 2\\ 0\hfill & 0& \hfill 2Q(\lambda )\\ \lambda \frac{dQ}{d\lambda }+Q(\lambda )\hfill & 1& \hfill 0\end{array}\right)$$
(65)
The eigenvalues of (65) are zero and:

$$\mu ^2=2\left(2Q+\lambda \frac{dQ}{d\lambda }\right)=2\frac{d\mathrm{\Psi }_{P_c}}{d\lambda }$$
(66)
Where we have used (63). We have proved that for the I Poisson Structure of dynamical system (56) the critical point is elliptic when $`d\mathrm{\Psi }/d\lambda >0`$ and hyperbolic when $`d\mathrm{\Psi }/d\beta <0`$.

II Poisson Structure $`H=C_2`$, and $`\mathrm{\Psi }=C_1`$.

$`J^{(2)ab}=\epsilon ^{abc}{\displaystyle \frac{\mathrm{\Psi }}{X^c}}`$

$$J^{(2)ab}=\left(\begin{array}{ccc}0\hfill & 2X_3& \hfill X_1\\ 2X_3\hfill & 0& \hfill X_2\\ X_1\hfill & X_2& \hfill 0\end{array}\right)$$
(67)
With the help of the Poisson structure (67) T. Bridges in  obtained the opposite criterion to the one proved above. This means that we can not use the slope of a Casimir function parametrized by the Lagrange multiplier to establish the critical point type, independently of the used hamiltonian structure.

## IV The Hamilton–Jacobi Equations

Consider the following equation of motion

$$\frac{dx^a}{dt}=f^a(x^b)a,b=1\mathrm{}N$$
(68)
If we have a hamiltonian structure on the phase space $`M`$ of (68) with a Casimir function $`\mathrm{\Psi }`$, then the Casimir determine the geometry of the phase space. $`\mathrm{\Psi }`$ is, of course, a constant of the motion of dynamical system (68). The equation $`\mathrm{\Psi }(x^b)=C`$ for different values of $`C`$ define a foliation of the phase space $`M`$. Each leaf of the foliation is an invariant submanifold, furthermore it is a symplectic submanifold $`N`$ of the phase space $`M`$. Using a Poisson map $`\phi `$, i. e. a transformation $`\phi :MN`$ that conserve the Poisson bracket: $`\{F(\phi ),H(\phi )\}_M=\{F,G\}_N(\phi )`$, for any functions $`F,H:N`$, we can define canonical variables, and an $`action`$ on $`N`$. The Hamilton–Jacobi (H–J) equation for this action yields a particular solution of dynamical system constrained to N. Let’s illustrate the above procedure using the first example (38) of the previous section with the II Poisson structure (49). For this srtucture we have the following Poisson brackets:

$`\{X_1,X_2\}_M=0`$ (69)
$`\{X_1,\mathrm{\Phi }\}_M={\displaystyle \frac{1}{X_1}}`$ (70)
$`\{X_2,\mathrm{\Phi }\}_M=\pm {\displaystyle \frac{1}{X_2}}`$ (71)
From (70) we have:

$$X_1\{X_1,\mathrm{\Phi }\}_M=\{X_1,X_1\mathrm{\Phi }\}=1$$
(72)
This suggest the Poisson map $`\phi :MN`$ , given by the inverse of:

$$X_1=q,\mathrm{\Phi }=\frac{p}{q},X_2=\sqrt{2C_1\pm q^2}$$
(73)
Where $`N`$ is the submanifold defined by: $`\mathrm{\Psi }(X_1,X_2,\mathrm{\Phi })=C_1`$

We can show that $`\phi `$ is a Poisson map, moreover it define canonical brackets on $`N`$:

$`\{X_1,X_1\mathrm{\Phi }\}_M=\{q,{\displaystyle \frac{pq}{q}}\}_N=\{q,p\}_N=1`$ (74)
$`\{X_1,X_2\}_M=\{q,\sqrt{2C_1\pm q^2}\}_N=0`$ (75)
$`\{X_2,\mathrm{\Phi }\}_M=\{\sqrt{2C_1\pm q^2},{\displaystyle \frac{p}{q}}\}_N=\pm {\displaystyle \frac{1}{\sqrt{2C_1\pm q^2}}}=\pm {\displaystyle \frac{1}{X_2}}`$ (76)
On $`N`$ in $`(p,q)`$ coordinates the Hamiltonian $`H`$ take the form:

$$H=q^2\sqrt{2C_1\pm q^2}\mathrm{sin}\frac{p}{q}$$
(77)
Using this Hamiltonian and introducing the action by the transformation: $`qq`$ and $`p\frac{S}{q}`$. We have the H–J:

$$\frac{S}{t}+q^2\sqrt{2C_1\pm q^2}\mathrm{sin}\frac{1}{q}\frac{S}{q}=0$$
(78)
Taking the solution of (78) in the form:

$`S=Q(q)Et`$

Replacing in the H–J equation we find:

$$tt_0=\frac{Q}{E}=_0^q\frac{\stackrel{´}{q}d\stackrel{´}{q}}{\sqrt{\stackrel{´}{q}^4(2C_1\pm \stackrel{´}{q}^2)E^2}}$$
(79)
The integral (79) is an elliptic integral and has solution in terms of special functions (The cnoidal and snoidal functions). The solutions $`q(t)`$ and from the Poisson map (73), the solutions of dynamical system (38) lie on concentric closed curves constrained to the submanifold $`N`$. We can deduce immediately that all solutions of (38) are periodic except for the stationary point (50).

Let’s take now the second example (56) of the previous section, with the Poisson structure (60). We have the following Poisson bracket on the phase space $`M`$:

$`\{X_1,X_3\}_M=1`$ (80)
$`\{X_2,X_3\}_M=Q(X_1)`$ (81)
$`\{X_1,X_2\}_M=0`$ (82)
Using the Poisson map $`\phi :MN`$ , given by the inverse of:

$$X_1=q,X_3=p,X_2=C_1+_0^qQ(s)𝑑s$$
(83)
where $`N`$ is the submanifold defined by: $`\mathrm{\Psi }(X_1,X_2,X_3)=C_1`$, we obtain, the canonical brackets on $`N`$ with the Hamiltonian:

$$H=p^2q\left(C_1+_0^qQ(s)𝑑s\right)$$
(84)
Introducing the action by the transformation: $`qq`$ and $`p\frac{S}{q}`$. We have the H–J equation:

$$\frac{S}{t}+\left(\frac{S}{q}\right)^2q\left(C_1+_0^qQ(s)𝑑s\right)=0$$
(85)
Taking the solution of (85) in the form:

$`S=K(q)Et`$

Replacing in the H–J equation we find:

$$2(tt_0)=\frac{K}{E}=_0^q\frac{d\stackrel{´}{q}}{\sqrt{\stackrel{´}{q}(C_1+_0^{\stackrel{´}{q}}Q(s)ds))+E}}$$
(86)
Taking as an example $`Q(s)=1`$ we find the following solution for (56):

$`X_1`$ $`=`$ $`{\displaystyle \frac{\left[Ae^{2(tt_0)}\frac{C_1}{2}\right]^2E}{Ae^{2(tt_0)}}}`$ (87)
$`X_2`$ $`=`$ $`{\displaystyle \frac{\left[Ae^{2(tt_0)}+\frac{C_1}{2}\right]^2E}{Ae^{2(tt_0)}}}`$ (88)
$`X_3`$ $`=`$ $`{\displaystyle \frac{\left[Ae^{2(tt_0)}\frac{C_1}{2}\right]\left[Ae^{2(tt_0)}+\frac{C_1}{2}\right]}{2Ae^{2(tt_0)}}}`$ (89)
Where the solution is obtained for: $`E\frac{C_1^2}{4}0`$ and $`A=\sqrt{E}+\frac{C_1}{2}`$ For other expressions of the polynomial $`Q(s)`$, the solutions of (56) are expressed in terms of elliptic functions.

Alternatively we can use, with this example the Poisson structure (67). For this structure we have the following Poisson brackets on $`M`$:

$`\{X_1,X_2\}_M=2X_3`$ (90)
$`\{X_1,X_3\}_M=X_1`$ (91)
$`\{X_3,X_2\}_M=X_3`$ (92)
Using the Poisson map $`\phi :MN`$ , given by the inverse of:

$$X_1=e^q,X_3=p,X_2=e^qp^2$$
(93)
where $`N`$ is the submanifold defined by: $`\mathrm{\Psi }(X_1,X_2,X_3)=C_2`$, we obtain the canonical brackets on $`N`$ with the Hamiltonian:

$$H=qp^2_0^{e^q}Q(s)𝑑s$$
(94)
The obtained H–J equation is:

$$\frac{S}{t}+e^q\left(\frac{S}{q}\right)^2_0^{e^q}Q(s)𝑑s=0$$
(95)
From (95) we find an implicit equation for $`q(t)`$:

$$2(t\stackrel{~}{t}_0)=_0^q\frac{d\stackrel{´}{q}}{\sqrt{e^{\stackrel{´}{q}}\left(_0^{e^{\stackrel{´}{q}}}Q(s)𝑑s+E\right)}}$$
(96)
and from $`q(t)`$ we obtain a solution of (56).

Finally we present a similar analysis for the Euler top. The free body equations of motion are:

$$\dot{\stackrel{}{L}}=\stackrel{}{\mathrm{\Omega }}\times \stackrel{}{L},\mathrm{\Omega }=\frac{\stackrel{}{L}}{I}$$
(97)
Where $`I=`$diag$`(I_1,I_2,I_3)`$ is the diagonalized moment of inertia tensor, $`I_1,I_2,I_3>0`$. $`\stackrel{}{L}`$ is the angular momentum, and $`\stackrel{}{\mathrm{\Omega }}`$ is the angular velocity.

The dynamical system (97) has the following conserved quantities:

$`C_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{L_1^2}{I_1}}+{\displaystyle \frac{L_2^2}{I_2}}+{\displaystyle \frac{L_3^2}{I_3}}\right)`$ (98)
$`C_2`$ $`=`$ $`L_1^2+L_2^2+L_3^2`$ (99)
This system is Hamiltonian with $`H=C_1`$, $`\mathrm{\Psi }=C_2`$ and the following Poisson matrix:

$`J^{(1)ab}={\displaystyle \frac{1}{2}}\epsilon ^{abc}{\displaystyle \frac{\mathrm{\Psi }}{X^c}}`$

$$J^{(1)ab}=\left(\begin{array}{ccc}0\hfill & L_3& \hfill L_2\\ L_3\hfill & 0& \hfill L_1\\ L_2\hfill & L_1& \hfill 0\end{array}\right)$$
(100)
The phase space of (97) is $`M=^3`$. $`M`$ is a Poisson manifold with the brackets:

$`\{L_1,L_2\}_M=L_3`$ (101)
$`\{L_2,L_3\}_M=L_1`$ (102)
$`\{L_3,L_1\}_M=L_2`$ (103)
Mapping $`\phi :MN`$, with the Poisson map defined by:

$`L_1`$ $`=`$ $`\sqrt{\lambda p^2}\mathrm{cos}q`$ (104)
$`L_2`$ $`=`$ $`\sqrt{\lambda p^2}\mathrm{sin}q`$ (105)
$`L_3`$ $`=`$ $`p`$ (106)
and $`N`$ being a submanifold defined by the equation: $`\mathrm{\Psi }(L_1,L_2,L_3)=C_2`$. We obtain the canonical brackets from (102), for the variables $`p`$ and $`q`$ on the symplectic submanifold $`N`$.

Introducing the action by the transformation: $`qq`$ and $`p\frac{S}{q}`$. We have the H–J equation:

$$\frac{S}{t}+\frac{1}{2I_1}\left[\lambda \left(\frac{S}{q}\right)^2\right]\mathrm{cos}^2q+\frac{1}{2I_2}\left[\lambda \left(\frac{S}{q}\right)^2\right]\mathrm{sin}^2q+\frac{1}{2I_3}\left(\frac{S}{q}\right)^2=0$$
(107)
Taking the solution of (107) in the form:

$`S=Q(q)Et`$

Replacing in the H–J equation and considering the conditions:

$`I_3>I_2>I_1,\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}EI_1<\lambda <2I_3E`$

We find the solution of (97):

$`L_1`$ $`=`$ $`I_1^{1/2}\left({\displaystyle \frac{2EI_3\lambda }{I_3I_2}}\right)cn\tau `$ (108)
$`L_2`$ $`=`$ $`I_2^{1/2}\left({\displaystyle \frac{2EI_3\lambda }{I_3I_2}}\right)sn\tau `$ (109)
$`L_3`$ $`=`$ $`I_3^{1/2}\left({\displaystyle \frac{\lambda 2EI_1}{I_3I_1}}\right)dn\tau `$ (110)
$`\tau `$ $`=`$ $`\sqrt{{\displaystyle \frac{(\lambda 2EI_1)(I_3I_2)}{I_1I_2I_3}}}(tt_0)`$ (111)
Where $`sn\tau `$, $`cn\tau `$, and $`dn\tau `$ are the Jacobi functions. Note that (109) is agree with the solution obtained by L. D. Landau in .

Another Poisson structure can be obtained taking a Hamiltonian $`H=\frac{1}{2}C_2`$, a Casimir $`\mathrm{\Psi }=C_1`$ and a Poisson matrix:

$`J^{(2)ab}=\epsilon ^{abc}{\displaystyle \frac{\mathrm{\Psi }}{X^c}}`$

$$J^{(2)ab}=\left(\begin{array}{ccc}0\hfill & \frac{L_3}{I_3}& \hfill \frac{L_2}{I_2}\\ \frac{L_3}{I_3}\hfill & 0& \hfill \frac{L_1}{I_1}\\ \frac{L_2}{I_2}\hfill & \frac{L_1}{I_1}& \hfill 0\end{array}\right)$$
(112)
The phase space of (97) with this structure, is a Poisson manifold with the brackets:

$`\{L_1,L_2\}_M`$ $`=`$ $`{\displaystyle \frac{L_3}{I_3}}`$ (113)
$`\{L_2,L_3\}_M`$ $`=`$ $`{\displaystyle \frac{L_1}{I_1}}`$ (114)
$`\{L_3,L_1\}_M`$ $`=`$ $`{\displaystyle \frac{L_2}{I_2}}`$ (115)
Using the Poisson map $`\phi :MN`$, defined by:

$`L_1`$ $`=`$ $`\sqrt{{\displaystyle \frac{\lambda p^2}{I_2I_3}}}\mathrm{cos}q`$ (116)
$`L_2`$ $`=`$ $`\sqrt{{\displaystyle \frac{\lambda p^2}{I_1I_2}}}\mathrm{sin}q`$ (117)
$`L_3`$ $`=`$ $`{\displaystyle \frac{p}{\sqrt{I_1I_2}}}`$ (118)
the phase space $`M`$ is reduced to the submanifold $`N`$, defined by the equation: $`\mathrm{\Psi }(L_1,L_2,L_3)=C_1`$. $`N`$ is symplectic in the variables $`p`$ and $`q`$. With the Poisson map (117) we obtain the canonical brackets for $`p`$ and $`q`$ from (114).

Introducing the action by the transformation: $`qq`$ and $`p\frac{S}{q}`$. We have the H–J equation:

$$\frac{S}{t}+\frac{1}{2I_2I_3}\left[\lambda \left(\frac{S}{q}\right)^2\right]\mathrm{cos}^2q+\frac{1}{2I_1I_3}\left[\lambda \left(\frac{S}{q}\right)^2\right]\mathrm{sin}^2q+\frac{1}{2I_1I_2}\left(\frac{S}{q}\right)^2=0$$
(119)
Taking the solution of (119) in the form:

$`S=Q(q)Et`$

Note that in this additive separation of the variables, the constant parameter $`E`$ is not the usual energy.

Replacing in the H–J equation (119) and considering the conditions:

$`I_3>I_2>I_1,\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}EI_1I_2<\lambda <2EI_3I_2,\lambda =2I_1I_2I_3C_1`$

We find the solution of (97):

$`L_1`$ $`=`$ $`{\displaystyle \frac{1}{I_2^{1/2}}}\left({\displaystyle \frac{\lambda 2EI_1I_2}{I_3I_1}}\right)cn\stackrel{~}{\tau }`$ (120)
$`L_2`$ $`=`$ $`{\displaystyle \frac{1}{I_1^{1/2}}}\left({\displaystyle \frac{\lambda 2EI_1I_2}{I_3I_2}}\right)sn\stackrel{~}{\tau }`$ (121)
$`L_3`$ $`=`$ $`{\displaystyle \frac{1}{I_2^{1/2}}}\left({\displaystyle \frac{2EI_2I_3\lambda }{I_3I_1}}\right)dn\stackrel{~}{\tau }`$ (122)
$`\stackrel{~}{\tau }`$ $`=`$ $`{\displaystyle \frac{1}{I_1I_2I_3}}\sqrt{I_1(I_3I_2)(2EI_2I_3\lambda )}(t\stackrel{~}{t_0})`$ (123)
It is a straightforward matter to realize that (121) is a solution of (97) and become to (109) taking:

$`\lambda =2I_1I_2I_3C_1,\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}E=C_2`$

## V Conclusions

In section II we have been able to construct two different structures based on one symmetry vector $`\eta `$ and two conserved quantities $`C_1`$ and $`C_2`$ for one dynamical system. The Poisson structures $`J^{(1)}`$ and $`J^{(2)}`$ are built as the antisymmetric product of the evolution vector $`f`$ and the symmetry vector $`\eta `$ normalized using the deformation $`K`$ and $`Q`$ of $`C_1`$ and $`C_2`$ respectively.

In section III we showed that the stability criterion involving the derivative of a Casimir function parametrized by a Lagrange multiplier strongly depends of the used Poisson structure. The situation is even worse, because we can construct infinitely many Poisson structures, using the Poisson matrix defined by:

$`J^{ab}=\mu \epsilon ^{abc}{\displaystyle \frac{D}{X^c}}`$

Where $`D=D(C_1,C_2)`$ is an arbitrary function of $`C_1`$ and $`C_2`$, and a Hamiltonian $`H=H(C_1,C_2)`$, which is a function independent of $`D`$. We have many Poisson structures depending on the different choises of arbritary functions $`D`$ and $`H`$.(See for details ).

In the section IV, the Hamilton–Jacobi equations for a number of simple classical systems were obtained using non–standard Hamiltonian, constructed with the Hojman’s techniques, developed in .

It is important to note that in the reduction of the phase space $`M`$ of a dynamical system, the different values of the Casimir functions define a foliation of $`M`$. Each $`sheet`$, specified by the values of the casimir contains entire trajectories. i. e. solution of the equation of motion, and cannot be connected to any other sheet, only states within a sheet are transformable into one another. In other words each selection of the Poisson structure and then each sheet yield a different physical system, though they all share the same equation of motion. Further, since classical H–J theory forms a link with quantum mechanics it is of particular importance to examine alternative appoaches of constructing H–J equation for one classical system, to see if they lead to a razonable quantum mechanical systems.

Alternative H–J equations for Euler top are explicitly exhibited and solved. We demonstrated that both obtained H–J equation yield the correct solutions for the equation of motion, thought the $`actions`$ must be different.

More examples will be discussed in forthcoming articles.

## VI Acknowledgments

M.H. is deeply grateful to Instituto de Cooperación Iberoamericana for the MUTIS Doctoral Fellowship, and to Universidad de Chile for partial financial support (Proyecto 073).