text:
# DISCOVERY OF A 75 DAY ORBIT IN XTE J1543-568

## 1 Introduction

The transient X-ray source XTE J1543-568 was discovered during a slew of the Proportional Counter Array (PCA) on the Rossi X-ray Timing Explorer (RXTE) on Feb. 6.84, 2000, when it had a 2 to 10 keV flux of $`2\times 10^{10}`$ erg cm-2s-1 (Marshall, Takeshima & In ’t Zand 2000). A subsequent pointed PCA observation on Feb. 10.69 showed a similar intensity and revealed a pulsar with a period of 27.12$`\pm 0.02`$ s. Later, the pulsar was found in earlier data from the Burst and Transient Source Experiment (BATSE) on the Compton Gamma-Ray Observatory; the first detection of the pulsar was on Jan. 25 (Finger & Wilson 2000). The most accurate location was recently determined by serendipitous observations with the Wide Field Cameras on BeppoSAX: $`\alpha _{2000.0}=15^\mathrm{h}44^\mathrm{m}01^\mathrm{s}`$, $`\delta _{2000.0}=56^\mathrm{o}45\stackrel{}{\mathrm{.}}9`$ with a 99% confidence region delimited by a circle with a radius of 2$`\stackrel{}{\mathrm{.}}0`$ (Kaptein, In ’t Zand & Heise 2001). Thus far, no optical counterpart has been reported.

XTE J1543-568 continues to be active one year after discovery and is being monitored through regular pointed PCA observations. We here discuss the results from a timing analysis of data from the first year of monitoring, and report the discovery of the binary orbit of XTE J1543-568, present the light curve and discuss implications for the companion of the neutron star.

## 2 Observations and data analysis

A target-of-opportunity program was initiated after the discovery with RXTE to monitor the pulsar and to try to detect Doppler delays from a binary orbit. The frequency of the observations went from once every day for the first 7 days, to once every 2 days during the next 10 days, to once every 3 days during the next 75 days, and finally to once every week for the remainder of the time. Weekly observations still continued one year after discovery. We here report about the first year of observations (i.e., up to Feb. 3, 2001) for which the total accumulated net exposure time is 172 ks. The number of observations is 77. Most observations were carried out over two consecutive RXTE orbits with a total span of 1.9 hr, but the time span has been as high as 5.5 hr. We discuss measurements that were carried out with the PCA on RXTE. This instrument (for a detailed description, see Jahoda et al. 1996) consists of an array of 5 co-aligned Proportional Counter Units (PCUs) that are sensitive to photons of energy 2 to 60 keV with a total collecting area of 6500 cm2. The spectral resolution is 18% full-width at half maximum (FWHM) at 6 keV and the field of view is 1o FWHM. During the observations the number of active PCUs varied between 1 and 5.

### 2.1 Light curve

We extracted spectra for each RXTE orbit of observation, for all active PCUs combined, and for layer 1 only. All spectra were successfully modeled between 3 and 20 keV using a power-law function with low-energy absorption (following Morrison & McCammon 1983), multiplied with an exponential cutoff function and supplemented with a narrow emission line near 6.5 keV which may be identified with Fe-K line emission for neutral to C-like ionization states and which has a constant equivalent width of $`127\pm 5`$ keV. Figure 1 presents the best-fit spectral parameters and 3-20 keV flux as a function of time, with a time resolution of one observation. The light curve exhibits variability on a time scale of weeks. The peak was reached about 90 days after discovery and is nine times the lowest flux levels observed thus far. Outside this main peak, the flux occasionally reaches half the peak value. There is no obvious periodicity.

We also studied the light curve as measured with the All-Sky Monitor on RXTE since 1996. This light curve is quite noisy because of the small flux. There is no conclusive evidence for outbursts other than the one reported here, and there is no indication that XTE J1543-568 was ever brighter than the flux levels measured with the PCA.

### 2.2 Pulse period and binary orbit

For the study of the pulse period, we employed event data with 0.5 s time resolution and 64 readout channels for photon energy (i.e., configuration E\_500us\_64M\_0\_1s). We optimized the accuracy of the pulse period determination as a function of detector layers and bandpasses and find the best result if all layers are included and photons with energies higher than 20 keV excluded, although this optimum is only marginally better than if only layer 1 is included. We combined the data from all active PCUs, applied standard filters to eliminate bad data, and corrected the resulting light curves to the solar system barycenter. The latter includes a correction for the satellite orbit. For each of the observations, we measured the pulse period by epoch folding the data for a relevant range of test periods and calculating the $`\chi ^2`$ values for a constant flux hypothesis. The resulting curves show sinc-like trends. To obtain the best-fit period we fitted for each observation the peak curve with a 3rd degree polynomial and identified the period with the peak position of the cubic function. The 68% confidence error was estimated by determining the range of periods for which the $`\chi ^2`$ value is within 1 from the maximum. This procedure was verified through Monte Carlo simulations. The resulting measurements are presented in figure 2. They are supplemented with 78 BATSE measurements from the public domain111 made available through the WWW at URL http://www.batse.msfc.nasa.gov/batse/ and from Finger & Wilson (2000), pertaining to TJD 11572-11678 (TJD = JD$`2,440,000.5`$). The BATSE observations stopped prior to de-orbiting the satellite in June 2000. Figure 2 clearly shows that there is a 75 d periodicity which must be due to an expected binary orbit. Also, long time scale spin-up trends are visible which are expected to result from accretion torque. The spin-up is faster in the first $`100`$ days which is consistent with an accretion torque explanation because the flux and, therefore, the accretion rate is larger at that time.

In order to disentangle the orbital modulation from the spin up trend, we modeled in an initial step the orbital modulation through a simple sinusoidal variation. The spin up was modeled through various simple functions with similar success and we here present the results when the period derivative $`\dot{P}`$ is described by a Gaussian function on top of a flat, non-zero, level. We only modeled the data after TJD 11610 because the spin up seems to be more complex before that. For the determination of the orbit this is not detrimental because the excluded time interval represents only about 5% of the complete time span of the observations. The results of the modeling are presented in figure 2. In a second step, we subtracted the model for the spin up from the data and fitted a model of an eccentric binary orbit to the residuals. The best-fit $`\chi _\nu ^2`$ is rather high at 10.2 (for 118 degrees of freedom $`\nu `$), which we attribute to an incomplete spin-up model. In order to obtain estimates for the errors in the model parameters that take into account the systematic uncertainties, we forced $`\chi _\nu ^2`$ to 1 through a multiplication of all pulse period errors with a constant factor $`\sqrt{10.2}=3.2`$. The results for the orbital and spin-up parameters are given in table 1. In figure 2 (bottom panel) we plot the data after subtraction of the orbital modulation model, in figure 3 after subtraction of the spin-up model and folding with the orbital period. As a double check, we carried out the same procedure for the data after TJD 11670 which excludes the times when the Gaussian component of the spin-up model is important. Thus, possible strong systematic errors are avoided at the expense of accuracy. The result, given in an extra column in table 1, is consistent with the previous result.

The pulse profile was studied for 3 observations throughout the total observation span, see figure 4 for the profile at flux peak. The profile consists of two broad peaks at a phase difference of 0.5 and with different magnitudes. The pulsed fraction is consistently 65$`\pm 5`$%. This kind of profile and amplitude is not unusual for pulsars among high-mass X-ray binaries (c.f., Bildsten et al. 1997 and White et al. 1995). We determined the profiles in two energy bands (2 to 10 keV and 10 to 20 keV) and do not find strong morphological differences. The pulse fraction is about 10% higher in the upper energy band.

## 3 Discussion

Given the projected semi-major axis $`a_\mathrm{x}\mathrm{sin}i`$ and $`P_{\mathrm{orb}}`$, we can place constraints on the companion of the neutron star. The mass function is

$`f(M)`$ $`=`$ $`{\displaystyle \frac{4\pi ^2(a_\mathrm{x}\mathrm{sin}i)^3}{GP_{\mathrm{orb}}^2}}={\displaystyle \frac{(M_\mathrm{c}\mathrm{sin}i)^3}{(M_\mathrm{x}+M_\mathrm{c})^2}}=8.2\pm 0.5\mathrm{M}_{}`$ (1)
If $`M_\mathrm{x}=1.4`$ M, this implies a minimum mass for the companion star of 10.0 M. For a main sequence star, the implied spectral type is B3 or earlier (Zombeck 1990). The orbital period places XTE J1543-568 in the ’Corbet’ diagram of pulse period versus orbital period unambiguously on the branch of optically confirmed Be transients (Corbet et al. 1986, c.f. Bildsten et al. 1997).

In principle the distance is constrained by the spin-up rate and flux. The flux is a measure of the mass accretion rate which determines the torque on the neutron star and the spin-up rate. Given the steadiness of the spin-up rate, we consider it likely that the accretion occurs through a disk. In that case Ghosh & Lamb (1979) predict a spin-up rate of

$`\dot{P}`$ $`=`$ $`1.9\times 10^{12}\mu _{30}^{2/7}m^{3/7}R_6^{6/7}I_{45}^1P^2L_{37}^{6/7}\mathrm{s}\mathrm{s}^1`$ (2)
for a neutron star with magnetic dipole moment $`\mu _{30}`$ in units of 1030 G cm3, mass $`m`$ in units of 1.4 M, radius $`R_6`$ in units of 10 km, moment of inertia $`I_{45}`$ in units of 1045 g cm2, and luminosity $`L_{37}`$ in units of 1037 erg s-1. In figure 5 we present the period derivative as a function of 3-20 keV flux. Our data are rather noisy and do not cover a large dynamic range in flux, but they appear to be roughly consistent with a 6/7 power law relationship. Given the spectral shapes found, the overall bolometric correction is about a factor of 2. If we calibrate the above $`\dot{P}/L_{37}`$ relation at the $`\dot{P}=1\times 10^8`$ s s-1 level and assume all normalized parameters to be 1, the implied distance is 26 kpc. This would place XTE J1543-568 at an unlikely position. Acceptable distances would perhaps be at least about 50% smaller. It is hard to achieve that by adjusting the other parameters individually: $`\mu _{30}`$ would have to be at least 64 times larger, $`m`$ 16 times smaller, $`R_6`$ 4 times larger, or $`I_{45}`$ 3 times smaller. Therefore, given the correctness of the Ghosh & Lamb formulation, there is at least the suggestion that the distance is large.

Another constraint for the distance, that does not require $`\mu `$, may be obtained from an argument over the magnetospheric radius $`r_{\mathrm{mag}}`$ at which the neutron star magnetosphere begins to dictate the mass flow. This radius scales with $`\dot{M}^{2/7}`$ or $`F^{2/7}`$. If the flux $`F`$ varies between $`F_{\mathrm{min}}`$ and $`F_{\mathrm{max}}`$ at times when the mass accretion persists, $`r_{\mathrm{mag}}<r_{\mathrm{co}}(F_{\mathrm{min}}/F_{\mathrm{max}})^{2/7}`$ at peak flux, where $`r_{\mathrm{co}}`$ is the co-rotation radius at which the Keplerian period is equal to the spin period ($`r_{\mathrm{co}}=(GMP^2/4\pi ^2)^{1/3}`$). The torque $`2\pi I\dot{\nu }`$ is given by $`\dot{M}\sqrt{GMr_{\mathrm{mag}}}`$. Substition in the above condition reveals

$`\dot{M}`$ $`>`$ $`5.6\times 10^2{\displaystyle \frac{I_{45}\dot{P}}{m^{2/3}P^{7/3}(F_{\mathrm{min}}/F_{\mathrm{max}})^{1/7}}}\mathrm{M}_{}\mathrm{yr}^1`$ (3)
For standard parameters and $`F_{\mathrm{max}}=4.5F_{\mathrm{min}}`$ (after TJD 11670), we obtain $`\dot{M}>9.8\times 10^{10}`$ Myr-1. If all liberated potential energy is released in the form of radiation, the bolometric luminosity is constrained by $`L>1.2\times 10^{37}`$ erg s-1. For $`F_{\mathrm{max}}=4.5\times 10^{10}`$ erg s-1cm-2 and a bolometric correction of 2, a minimum distance of 10 kpc is implied. Again this is a fairly large distance.

These distance estimates suggest that searches for the optical counterpart may prove difficult, also because the large value of $`N_\mathrm{H}`$ points to a high extinction. For a 10.0 M main sequence star, $`M_\mathrm{V}=2`$. If $`N_\mathrm{H}=1.3\times 10^{22}`$ (see figure 1; this is consistent with the total Galactic column as inferred from maps by Dickey & Lockman 1990), $`A_\mathrm{V}=7.4`$ (according to the standard conversions introduced by Predehl & Schmitt 1995). For a distance of 26 kpc, the apparent visual magnitude would be $`V=23`$, for a distance of 10 kpc, $`V=21`$.

The most interesting result of our analysis is the low eccentricity. It is among the lowest for the 12 confirmed or suspected Be X-ray binaries whose orbits have now been determined (8 cases are reviewed by Bildsten et al. 1997 and since then orbits have been determined for 2S 1845-024 by Finger et al. 1999, for SAX J2103.5+4545 by Baykal, Stark & Swank 2000, and for X Per by Delgado-Martí et al. 2001). There is one other system apart from XTE J1543-568 which has an eccentricity that is consistent with 0 ($`e<0.09`$ at $`2\sigma `$ confidence for 2S 1553-54, see Kelley et al. 1983). However, one should note that Kelley et al. have a sparse coverage of less than one binary orbit and do not solve for spin-up. Finding low-eccentricity systems is remarkable because one would expect most long-period Be systems to have undergone supernovae with large kick velocities for the neutron stars and to have circularization time scales far longer than their age. XTE J1543-568 confirms the suspicion that there may be a substantial population of neutron stars formed with little or no kick (Delgado-Martí et al. 2001).

We thank the project scientist Jean Swank for continuing the RXTE observations of XTE J1543-568 through public TOOs after our TOO program time was spent, Evan Smith for his tireless efforts to schedule the many observations, Mike Tripicco for his assistance in the barycentric corrections for 2001 PCA data, Ron Remillard for his advice on the ASM light curve, and the referee Mark Finger for valuable suggestions. JZ acknowledges financial support from the Netherlands Organization for Scientific Research (NWO).


format: text

text:
# Critical slowing down in polynomial time algorithms

## Abstract

Combinatorial optimization algorithms which compute exact ground state configurations in disordered magnets are seen to exhibit critical slowing down at zero temperature phase transitions. Using arguments based on the physical picture of the model, including vanishing stiffness on scales beyond the correlation length and the ground state degeneracy, the number of operations carried out by one such algorithm, the push-relabel algorithm for the random field Ising model, can be estimated. Some scaling can also be predicted for the 2D spin glass.

There has long been a close link between the concepts of statistical physics and the algorithms used to simulate condensed matter systems. The fundamental connection is the mathematics of graphs, which is applied analytically, for example, to compute high temperature series. Using sophisticated connectivity algorithms borrowed from computer science, one can numerically calculate quantities in percolation to a high precision . The Fortuin-Kastelyn cluster representation  for magnets is the basis for the Swendsen-Wang algorithm . By implementing nonlocal dynamics, these algorithms can greatly reduce the time needed for simulations near a phase transition. The study of disordered systems, such as spin glasses, pinned vortices in superconductors, and random field magnets, lead to the introduction of graphs with nonuniform edges . Early in the study of disordered systems, it was realized that the study of such graphs is directly related to issues of computational complexity. In some cases , computing the ground state of a disordered material is computationally intractable, as the relevant optimization questions on a graph are NP-hard . Quite interestingly, some computationally intractable problems have been found to have phase transitions . For example, given an ensemble of logical expressions characterized by the number of Boolean variables $`N`$ and clauses $`M`$, the fraction that are satisfiable can exhibit finite size scaling about a critical value of $`M/N`$.

In this letter, I present results on the behavior of ground state algorithms near phase transitions in two models, the random field Ising magnet (RFIM) and the 2D spin glass (2DSG), which, in contrast with problems such as satisfiability, are always solvable in time polynomial in the size of the graphs. These phase transitions lead to singularities in the mean time to solve the graph optimization problems. For the RFIM, a close connection is made between the critical slowing down of the ground state algorithm, the correlation lengths, and the degeneracy of the ground state in the thermodynamic limit. Numerical evidence is presented that suggests that the dynamic critical exponent is $`z1`$ for the RFIM; scaling arguments suggest $`z1`$. The behavior far from the transition can also be explained. The algorithm for the 2DSG is more difficult to analyze, but the observed uniform time per spin in the ferromagnetic (FM) phase is seen to be natural.

Generally, the dynamics at low temperature $`T`$ is exceedingly slow in disordered magnets, as seen in experiment and in Monte Carlo simulations using local moves . The glassy dynamical behavior is attributed to the complex structure of the energy functional. The free energy barriers to equilibration at a length scale $`\mathrm{}`$ scale as $`\mathrm{}^\psi `$, so that the time to equilibrate a portion of the sample of size $`\mathrm{}^d`$ are expected to scale as $`\mathrm{exp}(\mathrm{}^\psi /T)`$. However, for some random magnets there are algorithms which take time polynomial in $`N`$ to find an exact ground state. The process of finding the ground state uses ”nonphysical” configurations or moves: a local search or simulated annealing that uses only physical configurations and local moves is hindered by the large energy barriers.

When using local Monte Carlo moves at finite $`T`$ to model uniform magnets, the relaxation or correlation time at continuous transitions scales as $`L^z`$, with $`z_{\mathrm{loc}}\gamma /\nu `$ . Nonlocal cluster moves, such as used in the Swendsen-Wang algorithm, can reduce the critical exponent $`z`$, with $`z_{\mathrm{cl}}\alpha /\nu `$ , with $`\alpha `$ and $`\nu `$ the exponents for the heat capacity and correlation length, respectively. As the algorithms used to find the ground states of disordered magnets do not use local moves and are not designed to utilize clusters related to a critical point, it is less clear how many operations are required. Polynomial bounds on the worst case behavior of these algorithms do exist. For graphs with $`n`$ vertices and $`m`$ edges, the highest level version of the push-relabel (PR) algorithm, used here for the RFIM, will use $`O(n^2\sqrt{m})`$ time. Useful algorithms for general matching , applied here to the 2DSG, use from $`O(n^{3/2})`$ to $`O(n^3)`$ operations, up to logarithmic corrections, assuming $`mn`$. In practice, however, for typical disorder realizations of physical interest in finite dimension $`d`$, these algorithms are much faster, with the average running time for many problems scaling roughly as $`N^q`$, with $`q`$ typically in the range $`1.11.3`$ . For systems with a single $`T=0`$ fixed point, such as the elastic medium in a random potential , this scaling for the time would not be expected to vary with parameters.

The $`T=0`$ RFIM has a phase transition that has been extensively studied using a mapping  of the ground state to the optimization problem max-flow . The Gaussian RFIM has Hamiltonian $`H_R=J_{<ij>}s_is_j_ih_is_i`$, where the allowed spin values on $`N=L^d`$ lattice sites $`i`$ are $`s_i=\pm 1`$, the ferromagnetic couplings $`J`$ are positive, and the random fields $`h_i`$ are Gaussian distributed, with mean $`0`$ and variance $`\mathrm{\Delta }^2J^2`$. In the paramagnetic phase ($`\mathrm{\Delta }>\mathrm{\Delta }_c2.27`$), the spins are pinned by the external fields $`h_i`$ and the magnitude of the net magnetization $`m=N^1s_i`$ vanishes as $`L0`$. In the ferromagnetic phase ($`\mathrm{\Delta }<\mathrm{\Delta }_c`$), the ferromagnetic coupling $`J`$ dominates and $`|m|`$ has a finite limit, with the sign of $`m`$ determined roughly by the total of $`h_i`$ (some spins will be reversed by fluctuations in $`h_i`$.) Results for the time for this algorithm to find the ground state in a 3D cubic lattice with periodic boundaries are plotted as a function of $`\mathrm{\Delta }`$ and linear system size $`L`$ in Fig. 1. Plots of the number of primitive operations executed are nearly identical in form. Near the phase transition separating the ferromagnetic from paramagnetic phase, the time to find the ground state shows a characteristic critical slowing down. This notable effect reflects the deep connection between the dynamics of the ground-state algorithm and the physics of the RFIM.

The application of max-flow algorithms to the RFIM is well established, but to make the connection between scaling in the RFIM and algorithm timing, it is useful to review the algorithm. The network flow algorithm generally used to solve the RFIM, because of its speed, is the PR algorithm of Tarjan and Goldberg . The algorithm described here uses a modification that removes the need for source and sink nodes , reducing memory usage and also clarifying the physical connection.

The modified PR algorithm starts by assigning an “excess” $`x_i`$ to each lattice site $`i`$, with $`x_i=h_i`$. Residual capacity variables $`r_{ij}`$ between neighboring sites are initially set to $`J`$. A height variable $`u_i`$ is then assigned to each node via a global update step. In this global update, the value of $`u_i`$ at each site in the set $`𝒯=\left\{j\right|x_j<0\}`$ of negative excess sites is set to zero. Sites with $`x_i0`$ have $`u_i`$ set to the length of the shortest path, via edges with positive capacity, from $`i`$ to $`𝒯`$.

The ground state is found by successively rearranging the excesses $`x_i`$, via “push” operations, and updating the heights, via “relabel” operations. When no more pushes or relabels are possible, a final global update determines the ground state: those sites which are path connected by bonds with $`r_{ij}>0`$ to $`𝒯`$ have $`s_i=1`$, while the sites which are disconnected from $`𝒯`$ have $`s_i=1`$. A push operation moves excess from a site $`i`$ to a lower height neighbor $`j`$, if possible, that is, whenever $`x_i>0`$, $`r_{ij}>0`$ and $`u_j=u_i1`$. In a push, the working variables are modified according to $`x_ix_i\delta `$, $`x_jx_j+\delta `$, $`r_{ij}r_{ij}\delta `$, and $`r_{ji}r_{ji}+\delta `$, with $`\delta =\mathrm{min}(x_i,r_{ij})`$. Push operations tend to move the positive excess towards sites in $`𝒯`$. When $`x_i>0`$ and no push is possible, the site is relabeled, with $`u_i`$ increased to $`1+\mathrm{max}_{\{j|r_{ij}>0\}}u_j`$. In addition, if a set of highest sites $`𝒰`$ become isolated, with $`u_i>u_j+1`$, for all $`i𝒰`$ and all $`j𝒰`$, the height $`u_i`$ for all $`i𝒰`$ is increased to its maximum value, $`N`$, as these sites will always be isolated from the negative excess nodes. A proof of the correctness of the PR flow algorithm can be found in standard textbooks  and its application to the RFIM is well known . Periodic global updates, here applied every $`N`$ relabels, are often crucial to the practical speed of the algorithm. The highest site heuristic is used here, which applies pushes and relabels where $`u_i`$ is maximal and $`x_i>0`$.

The PR algorithm is intuitively appealing: when the initial capacities within a region are large, the excess can be rearranged so that the positive excesses cancel the negative excesses as much as possible. The remaining excess values having the sign of the original total excess for the region. As $`r_{ij}=J`$ and $`x_i=h_i`$ initially, large capacities correspond to the ferromagnetic bonds being strong enough to favor alignments of the spins, with the spin direction given by the sign of the total $`h_i`$ for the region. If the initial capacities are not large enough (weaker $`J`$), the regions align independently, according to the local field, and the excesses do not cancel. The number of steps needed to move excess across the diameter of a region via push operations is bounded below by the linear size of the region. Note that the running time in the FM phase is somewhat dependent on the majority magnetization, due to the up/down asymmetry of the algorithm.

I now argue that these correspondences can be used to bound the running time of the algorithm, using the physical properties of the RFIM, particularly the behavior of the ground state degeneracy in the thermodynamic limit. Recent work  has given strong numerical evidence of insensitivity of the interior spins in the ground state to boundary conditions, when $`\mathrm{\Delta }>\mathrm{\Delta }_c`$. (This is in contrast with the scenario for, e.g., mean field  or highly disordered  spin glasses in $`d>8`$, where the entire solution is sensitive to the volume.) Most importantly, this implies that the ground state solution is determined by the $`h_i`$ within a volume typically of the size of the correlation volume $`\xi ^d`$. However, in the FM phase, the interior is sensitive to the boundaries and a finite fraction of the spins flip infinitely often as the sample volume is increased, as the net magnetization is given by a coarse-grained global sum of the $`h_i`$.

First, consider the case $`\mathrm{\Delta }2d`$, where, in the ground state, almost all spins satisfy $`s_i=\mathrm{sgn}(h_i)`$. In the algorithm, the initial positive excesses at sites $`i`$ where $`h_i>2dJ`$ remain positive, as the total capacity of bonds leaving the site is only $`2dJ`$, so that push operations can only reduce $`h_i`$ by $`2dJ`$. The sites with excess $`x_i<2dJ`$ also always have negative excess. There will be some small rearrangement, but only locally, and the sign of the excess will change only at a very few sites during the execution of the algorithm, which terminates after a number of operations $`N`$ (up to logarithmic corrections .)

A similar scenario holds, but at scales $`\xi `$, for $`\mathrm{\Delta }>\mathrm{\Delta }_c`$. The algorithm establishes the boundaries of correlation volumes by rearrangement of excess over distances of scale $`\xi `$. Further rearrangement is blocked by the effective decrease in residual capacity with scale (as the stiffness, corresponding to the scale dependent $`J`$, decreases rapidly on scales greater than $`\xi `$ .) The question to be answered, then, is how long the algorithm takes, per spin, to redistribute excess on scale $`\xi `$. The number of push and relabel operations in a volume $`\xi ^d`$ is bounded below by $`\xi ^{d+1}`$. For the excess to be pushed over a distance $`\xi `$, the relative heights must differ by at least $`\xi `$, so that at least $`O(\xi ^{d+1})`$ relabels $`R`$ must be performed for each correlation volume (global updates lead to height changes, but these do not appear to affect the scaling, empirically.) This gives the estimate $`RL^3\xi `$. This scaling is consistent with numerical results, up to logarithmic corrections, for $`d=1,3`$. The inset to Fig. 1 show a scaling plot for $`R`$, for example, with $`(RL^3)^{1/\nu }\delta `$ plotted as a function of the finite size scaling variable $`\delta L^{1/\nu }`$, $`\delta =(\mathrm{\Delta }\mathrm{\Delta }_c)/\mathrm{\Delta }_c`$, with the values $`\nu =1.37`$ and $`\mathrm{\Delta }_c=2.27`$ fixed parameters, determined independently . A fit in $`d=1`$, with $`\nu =2`$ and $`\mathrm{\Delta }_c=0`$ also fixed by known results, describes the data well for $`L5\times 10^6`$, with $`RL\xi \mathrm{ln}(L/\xi )`$ (without the global relabeling heuristic.)

One limit where the asymptotic time for the algorithm can be described more precisely is when $`\mathrm{\Delta }[L\mathrm{ln}(L)]^{1/2}`$.) Here, the capacities do not limit the rearrangement of excess: the final state is either $`x_i0`$ everywhere or $`x_i0`$ everywhere, corresponding to a uniform $`s_i=\pm 1`$ state, according to whether $`h_i0`$ or $`h_i0`$, respectively. The dynamics of the PR algorithm is set by the fluctuations in the total $`h_i`$ in regions at each length scale. At any time scale of the computation, the positive excess will be pushed towards the nearest negative excess region, with the distribution of excess negative or positive over a particular volume with $`\mathrm{}^d`$ sites, as the rearrangement of excess will have been completed over shorter lengths at an earlier time scale. Generally, but especially given sufficiently frequent global updates, the sites with greatest $`u_i`$ will be furthest from the set $`𝒯`$. As the highest sites are examined for pushing, the excess will be moved from these sites to the next lower level. These sites will then have their excess moved to the next lower level and the algorithm will “sweep” the excess through the volume of diameter $`2\mathrm{}`$. This will establish a distribution of excess that will have uniform sign over a region of size $`2\mathrm{}`$. In this fashion, the algorithm will solve for the sum of $`h_i`$ recursively. The number of steps at each stage will be $`L^d`$ and $`O(\mathrm{ln}L)`$ stages will be required, giving a total running time $`L^d\mathrm{ln}(L)`$. This result is consistent with the numerical timings for very small $`\mathrm{\Delta }`$, with the data for $`RL^d`$ linear in $`\mathrm{ln}(L)`$ over sample dimensions $`16L128`$ to within the $`1\%`$ numerical error for $`d=3`$ and over $`10^3L5\times 10^6`$ to within the same error for $`d=1`$. Coarsening of the height variables during the algorithm is displayed in Fig. 2.

To indicate the generality of critical slowing down and partial arguments for other systems, it is useful to compare the RFIM results with results for the 2DSG. The Hamiltonian is $`H_S=_{<ij>}J_{ij}s_is_j`$, where again the $`s_i`$ are Ising spins. The Gaussian distributed $`J_{ij}`$ have mean $`J_0`$ and variance $`1`$. The SG to FM transition  takes place as $`J_0`$ is increased through the critical value $`J_c0.96`$. The mapping from the 2DSG with free boundaries to a general matching problem is given in Ref. : energy is minimized by selecting a state with a minimum total weight for frustrated bonds. Timing results for the 2DSG are shown in Fig. 3. One notable difference from the RFIM is the apparent convergence to constant time per spin in the FM phase. As the algorithm used for the 2DSG uses a bond representation, the algorithm does not need to distinguish degenerate states related by global spin flips. If the frustration is low enough, the FM phase is obtained by local operations giving a solution with percolating unfrustrated bonds. For low $`J_0`$, in the SG phase, in contrast, though locally the ground state is insensitive to the boundaries, the global ground state is sensitive to the disorder and the operations (augmenting paths) must therefore be carried out over all scales, so that the time per spin may diverge as a power of $`L`$ (numerically, approximately as $`L^{0.78\pm 0.10}`$ for $`96L720`$ and $`J_0=0`$.) At the critical point, independent calculations  show that the domain wall fractal dimension increases slightly, so that bond updates on large scales are more expensive.

In summary, the time for a polynomial time algorithm to find ground states is examined near zero temperature critical points for two models, the random field Ising model and the 2D spin glass. At the critical points, the combinatorial optimization algorithms used, PR flow and general matching, while exponentially faster than, say, simulated annealing at finding the ground state, slows down. This slowing down of the algorithm is argued to be closely related to important physical ideas, namely, the uniqueness or two-fold degeneracy of the ground state in various phases and the divergence of the correlation length as the transition is approached. The time for the algorithm can be understood in detail for the RFIM with small random fields. It will certainly be of interest to consider such slowing down and scaling at other fixed points for other polynomial time algorithms, to expand our understanding of the physics and algorithms for these systems. These considerations are clearly relevant for developing efficient parallel algorithms. I would like to thank J. Machta for useful discussions. This work was supported in part by the NSF (DMR-9702242.)


text:
# How large are the spectral gaps?

## Proof of Theorem 1

Since $`E_\mathrm{\Lambda }`$ is a frame for $`L^2(D)`$ if and only if $`E_{\mathrm{\Lambda }\mu }`$ is also a frame for $`L^2(D)`$ (with the same frame constants) for any $`\mu ^n`$, and our estimates do not depend on the choice of $`\mathrm{\Lambda }`$, it is sufficient to consider the case $`\mu =(0,\mathrm{},0)`$.

By the frame property and Lemma 7 we get

$$A|D|\underset{\mathrm{\Lambda }}{}|F_D\chi _D(\lambda )|^2=\underset{\mathrm{\Lambda }}{}|\widehat{\chi }_D(\lambda )|^2=\underset{Q_R}{}|\widehat{\chi }_D(\lambda )|^2+\underset{\{\lambda Q_R\}}{}|\widehat{\chi }_D(\lambda )|^2.$$
$`25`$
Using Lemma 9 we see that if $`R=2^{k_0}`$,

$$\underset{\{\lambda Q_R\}}{}|\widehat{\chi }_D(\lambda )|^2=\underset{k=k_0}{\overset{\mathrm{}}{}}\underset{Q_{2^{k+1}}Q_{2^k}}{}|\widehat{\chi }_D(\lambda )|^2CB2^{k_0(n\alpha )}=\frac{BC}{R^{n\alpha }}.$$
$`26`$
So by $`(25)`$ and $`(26)`$

$$\underset{Q_R}{}|\widehat{\chi }_D(\lambda )|^2A|D|\frac{BC}{R^{n\alpha }}$$
$`27`$
which proves that if $`R>\left(\frac{BC}{A|D|}\right)^{\frac{1}{n\alpha }}`$, then

$$\mathrm{\Lambda }Q_R\mathrm{}.$$
$`28`$
Moreover, the above proof shows that $`CC^{}|D|_\alpha `$ where $`C^{}`$ depends only on $`n`$.

###### Remark Remark

In the proof above the key estimate is $`|\{DDh\}|C|h|^{n\alpha }|D|_\alpha `$. While this is the best possible estimate uniform in $`h`$, in the proof we are have a wide choice of $`h`$’s as long as $`|h|=2^k`$ and the estimates $`(18)`$, $`(19)`$, and $`(21)`$ are satisfied.

This observation can be used to handle the family of examples given by Example 5 above. For convenience we take $`\mathrm{\Lambda }=(\frac{1}{2},0)+^2`$. We can now take all $`h`$’s in the proof of Theorem 1 of the form $`h=(h_1,0)`$ and for this choice of $`h`$’s it is easy to check that $`|\{D_kD_kh\}|C|h|diameter(D_k)`$, where $`C`$ is a uniform constant, since the ”teeth” of $`D_k`$’s point in the $`y`$-direction. Since $`diameter(D_k)`$ is uniformly bounded above and below, the lack of sharpness of Theorem 1 exposed in Example 5 is resolved for this family of examples.

References


text:
# Suppression of Superconducting Critical Current Density by Small Flux Jumps in 𝑀𝑔𝐵₂ Thin Films

## I Introduction

For applications of superconductors, a high transition temperature $`T_c`$ and superconducting critical current density $`j_c`$ are desirable. The discovery of the new superconductor $`MgB_2`$  with a remarkably high transition temperature of 39 K and high critical current density $`j_c`$ provides a promising candidate for such applications. One big issue concerns the stability of the critical current and what process dominates the flux motion in the new superconductor $`MgB_2`$. For bulk samples the flux dynamics and the vortex phase diagram have been intensively investigated . Usually it is believed that the $`j_c`$ is controlled by either continuous flux creep ( in intermediate and high temperature region ) or sudden, large and discontinuous flux jumps ( in low temperature region ). The former can give rise to a continuously and parallelly marching flux front and thus a continuous magnetization vs. field M(H) curve, whereas the latter will generate many big blasts at the flux front leading to big discontinuous steps on the M(H) curves. Since there are generally more defects ( so more pinning centers ) in thin films, it is interesting to know whether the same flux dynamics is occurring in thin films as in bulk samples. In this paper, we present the experimental observation of many small flux jumps ( SFJ ) in low temperature and field region of $`MgB_2`$ thin films. It is further shown that the central magnetization peak of the magnetization-hysteresis-loop ( MHL ) is smeared out and the nominal $`j_c`$ is suppressed by this effect in comparison to that due to the continuous flux creep.

## II Experiment

The thin films of $`MgB_2`$ were fabricated on (1102) $`Al_2O_3`$ substrates by using pulsed laser deposition technique, which was described clearly in Ref.. They are typically 400 nm thick with predominant c-axis orientation ( c-axis is perpendicular to the film surface ). A rectangular sample of size $`2.1mm\times 4.9mm`$ was chosen for the magnetic measurements. The temperature dependence of the diamagnetic moment was carried out by a Quantum Design superconducting quantum interference device ( SQUID, MPMS 5.5 T ) and the magnetization-hysteresis-loops ( MHL ) were measured with a vibrating sample magnetometer ( VSM 8T, Oxford 3001 ) at temperatures ranging from 2 K to $`T_c`$ and external field up to 8 T along the c-axis. The M(H) curve was measured with a field sweep rate of 0.01 T/s and integration time of 60 ms. The pressure of helium gas in the sample chamber for thermal exchange was kept at 0.04 bar during the measurement.

## III Results

In the inset of Fig.1, we show the temperature dependence of the zero-field-cooled (ZFC) magnetization measursed by SQUID at $`\mu _0`$H = 0.001 T. It is clear that the superconducting transition temperature $`T_c`$ is about 38 K and the transition is rather sharp indicating a good quality of the film.

The main frame of Fig.1 shows a typical MHL measured at 2 K with the external field sweep rate of 0.01 T/s. One can see that the MHL is symmetric about the M = 0 axis showing the dominance of the bulk superconducting current here. Another interesting finding is that the MHL is almost closed at about 8 T, being very similar to the observation on bulk $`MgB_2`$. This relatively low irreversibility field $`H_{irr}`$ determined from the closing point of MHL was first attributed by Wen et al. to the existence of quantum vortex liquid in rather clean system of $`MgB_2`$. The nominal $`j_c`$ estimated using the Bean critical state model is about 1.2$`\times 10^7A/cm^2`$ at 2 K and zero field which is about one order of magnitude higher than that in high-pressure synthesized bulk samples. This further indicates a rather good quality of the film. Worthy of noting here is that multiple small and irregular instabilities appear in the magnetization in a low field region. These irregular instabilities have magnitudes typically in the order of $`10^3`$ emu. They are much higher than the noise background ( $`10^5`$ emu ) of our VSM. The magnitude of these instabilities are much larger in low fields than in high fields ( where the M(H) curve becomes smooth ) although the field sweep rate and the data acquisition speed are the same indicating that the instabilities are not due to the noise background of the VSM. These instabilities have been rechecked with caution in a sense that after one week and two weeks, the same MHLs were obtained during the remeasuring process on $`MgB_2`$ thin films, confirming the reproducibility of the instabilities in $`MgB_2`$ thin films.

In order to know whether these irregular instabilities are related to the flux jumps appearing in bulk $`MgB_2`$ samples, we have carried out detailed measurements at different temperatures. The results are shown in Fig.2. One can clearly see the following interesting findings: (i) The instabilities in the magnetization appear in a region below a certain value of field and temperature, e.g. below 1.3 T at 2 K; (ii) When these instabilities appear the central magnetization peak is strongly flattened out ( see e.g. the data of 2, 4, 6, and 8 K ) revealing a suppression of the nominal $`j_c`$; (iii) The magnetizations at 2, 4, 6 and 8 K nearly merge at low field, though at higher fields they are separated gradually showing the usual order of the $`j_c`$ vs. H and T, i.e., a higher $`j_c`$ at a lower field and temperature; (4) When the temperature is increased the small instabilities will evolve into some larger ones ( at 10 K ) and disappear completely at a higher temperature ( at 14 K ). The larger instabilities here look similar to those appearing in bulk samples, therefore it is tempting to regard these instabilities as small flux jumps.

In Fig.3 we show a part of MHL measured at 4 K. One can clearly see that when the field is swept through 0.7 T from above, the continuous MHL curve becomes discontinuous: small and irregular instabilities appear in low field region. Interestingly these two regions are separated by a clear kink point on the MHL. The same feature appears for other MHLs measured at 2, 6 and 8 K. This kink may be understood in the following way: the SFJ appearing in low field region suppresses the magnetic moment ( thus the nominal critical current density ) of the sample; while in high field region the M(H) curve will behave in another way since no SFJ occurs there. However we are not sure whether this kink point is corresponding to a phase transition of the vortex system.

## IV Discussion

Fig.4 shows the field dependence of the nominal $`j_c`$ determined using the Bean critical state model via $`j_c=20\mathrm{\Delta }M/Va(1a/3b)`$, where $`\mathrm{\Delta }M`$ is the width of the MHL, V, a and b are the volume, width and length ( $`a<b`$ ) of the sample, respectively. Although for system with flux jumps the Bean critical state model may be inapplicable, it can be used to estimate the nominal $`j_c`$ value for a qualitative comparison. At 14 K and near zero field, the estimated $`j_c`$ is as high as $`1.7\times 10^7A/cm^2`$, which is rather high. It is necessary to note that the $`j_c`$ value at 14 K is a real magnetic critical current density instead of a nominal one since here no flux jumps occur. However, at lower temperatures, the nominal $`j_c`$ is clearly suppressed by these SFJ in low field region. In the high field and temperature region the M(H) curve becomes continuous showing the gradual setting in of the normal flux creep.

To get a more comprehensive understanding to the SFJ observed in a $`MgB_2`$ thin film, it is worthwhile to compare it with that in bulk $`MgB_2`$. The MHLs of a high pressure synthesized bulk sample are presented in Fig.5. It is clear that after each flux jump, the optimal slope may not fully recovered, while there is no suppression of the central magnetization peak ( and thus the nominal $`j_c`$ ) near zero field at 2 and 4 K in contrast to what observed in $`MgB_2`$ films. In other words, although strong flux jumps occur at 2 and 4 K in bulk, the outlines of the MHLs at these temperatures are still wider than those at higher temperatures where no flux jumps are observed.

Now we have a close look at the part of MHLs where the flux jumps occur. In Fig.6(a) and Fig.6(b) the magnetization M and the derivative dM/dH are plotted against the field in the flux jump region for the bulk and the thin film respectively. Two types of flux jumps can be clearly seen here for these two different samples. The jumps in the bulk sample are sparse and relatively large. After each jump, the magnetization changes significantly and then it gradually comes back to the main branch of MHL. A nearly constant field interval for every jump ( about 0.07 T ) is observed. During multiple measurements, the jumps repeatably tend to occur at the same field with the same magnitude. All these can be qualitatively understood based on the Swartz and Bean’s adiabatic theory. The flux jump can be triggered when the gradient of magnetic flux profile inside the sample exceeds some critical value. After very short time ( usually in the order of $`ms`$  ) the magnetic profile drops to a new configuration with a lower gradient, then the flux only creeps slowly until the gradient of magnetic flux profile induced by varying the external field exceeds the critical value again. In each jump many vortices are involved in the thermomagnetic avalanche which normally expands to a large part of the sample volume. In thin films, however, the situation is completely different: many small local avalanches occur ( as shown in Fig.6(b)). Although the H-T region for flux jumps to appear doesn’t change in different round of measurement of in $`MgB_2`$ thin films under identical conditions, the specific positions and the magnitude of the SFJ are however completely unrepeatable. All these cannot be explained by the adiabatic theory. The different avalanches observed in bulk samples and thin films may be induced by the different structural details and thermal diffusibility. For example, in bulk samples, there are many large grain boundaries, which act as strong pinning centers. The gradient of the flux profile near these boundaries can be broken at a certain limit. Once a blast occurs in a bulk sample, the thermal energy induced by drastic flux motion cannot easily diffuse out and be carried away by the environment. Therefore this self-heating will lead to increase the region in which the vortex instability occurs leading to a large jump on the magnetization in a bulk superconductor. One can also understand from this picture that the number of flux jumps cannot be large in bulk samples. In thin films the situation can be very different. On one hand scanning-electron-microscope ( SEM ) data indicates a high density of small defects formed during the preparation process of the thin films leading to much stronger critical current densities. Therefore there are many places for the avalanche to occur. On the other hand thermal diffusion is much easier in thin film samples due to their very small thickness and large surface area exposed to the environment. Therefore in thin films each avalanche is small in magnitude but the number of avalanches can be huge. This picture may give an explanation to many small vortex avalanches observed in Nb-film and the YBCO film. However why some of these small avalanches will grow in a dendritic structure is still an open question. But clearly the very fine disorder structure and the relatively better thermal diffusion in thin films are two key factors to be considered here.

Finally we suggest a criterion to identify the specific region for these SFJ on the H-T vortex phase diagram. The boundary point for SFJ is defined as the clear kink point shown in Fig.3. Above this field no flux jumps could be observed above the noise background of the instrument. This is very helpful to illustrate the field and temperature region in which the application of the $`MgB_2`$ film can be hampered by the thermal instability. In Fig.7 a borderline $`H_{SFJ}(T)`$ for the SFJ is plotted together with the irreversibility line $`H_{irr}(T)`$ determined from the closing point of the MHL with a criterion of $`\mathrm{\Delta }M=10^4emu`$. The SFJ appears only in the low temperature and low field region. Beyond this region the M(H) curves become continuously showing normal flux creep with very low creep rate. Therefore it is safe to conclude that the application in large part of the vortex solid will not be influenced by the SFJ. However, it is important to point out that the SFJ observed in these films will prevent using SQUID device made from these thin films at low temperatures.

We have been aware of a recent result by Johansen et al.  who found the SFJ and some dendritic avalanches in $`MgB_2`$ thin films by doing the magneto-optical measurement. One can see from their data that many tiny avalanches occur first at the edge of the film and some of them will gradually grow into a dendritic structure. Each tiny jump on our MHL may correspond to a local avalanche or the growth on one branch of the dendritic structure. These tiny jumps cannot be observed from the MHL measured by SQUID as presented by Johansen et al. since the SQUID has a very low speed for data acquisition. They can be seen clearly from our data measured by VSM with a fast data reading capacity.

## V Concluding remarks

In conclusion, we have reported the observation of the suppression of the central magnetization peak and thus the nominal critical current density in the new superconductors $`MgB_2`$ film at low temperatures due to many small flux jumps. A comparison with a $`MgB_2`$ bulk is made. It is suggested that the small vortex avalanches in the thin film are closely related to the high density of small defects and the relatively easy thermal diffusion. A borderline for this effect to occur is determined on the H-T phase diagram. This gives helpful information for the application of $`MgB_2`$ films.

###### Acknowledgements.

This work is supported by the National Science Foundation of China (NSFC 19825111) and the Ministry of Science and Technology of China ( project: NKBRSF-G1999064602 ). The work at Pohang University was supported by the Ministry of Science and Technology of Korea through the Creative Research Initiative Program.


text:
# Late Light Curves of Type Ia Supernovae

## 1 Introduction

Type Ia supernovae (SNe) are an integral part of many current astrophysical investigations. Their luminosities briefly rival entire galaxies, making them useful as distance indicators to large $`z`$ (Riess et al. 1998; Perlmutter et al. 1999). Understanding of their contributions to galactic chemical evolution is necessary to sort out the star formation history of various other stellar populations in the Galaxy (Timmes et al. 1996). The escape of gamma rays from their ejecta may produce a significant portion of the diffuse gamma ray background in the 400 - 2000 keV energy range (Watanabe et al. 1999). The escape of positrons from their ejecta may be large enough to explain a majority of the 511 keV annihilation radiation measured by the CGRO/OSSE instrument (Chan & Lingenfelter 1993; Milne, The & Leising 1999, hereafter MTL). All of these investigations rely upon a basic understanding of the SN Ia event.

Type Ia SNe are widely-believed to be the thermonuclear explosion of a mass-accreting carbon-oxygen white dwarf (WD). The conversion of carbon and oxygen to iron-peak and intermediate-mass elements releases a large amount of energy initially, and subsequently the decay 56Ni $``$ 56Co $``$ 56Fe powers the supernova through gamma- and X-ray photons and positrons.222Although neutrinos are also created in these decays, the lower neutrino opacity makes neutrinos an insignificant contributor to the energy deposition rate during the epoch of interest in this work. See Nadyozhin (1994) for major features of this decay chain. The deposition and diffusion of this energy produces the optical emission that defines the Type Ia SN. The observational support for this paradigm is typified by the ability of a standard carbon-oxygen WD model, W7, to reproduce the light curves and spectra of the SN Ia 1981B (Branch et al. 1985).

During the earliest epoch after the SN explosion, the energy released in the 56Ni $``$ 56Co $``$ 56Fe decay is efficiently trapped by the SN ejecta. The deposition of this energy creates energetic electrons, which thermalize (primarily via ionization and excitation of atoms) and recombine. The optical photons created in the recombinations diffuse to the surface, and escape the SN ejecta. The emerging optical spectrum is dictated by the nature of the diffusion of the optical photons. The expansion of the SN ejecta lowers the column density to the surface, decreasing both the time for optical photons to diffuse outward and the efficiency of trapping of gamma-rays and energetic positrons and electrons. The SN thus transitions from an epoch during which the energy deposition is essentially complete and instantaneous and the emission depends upon the diffusion of optical photons (a “diffusion-dominated” epoch) to an epoch during which the diffusion time-scale is negligibly short and the emerging emission depends upon the transport of the decay products (a “deposition-dominated” epoch).

Observations of SNe during these two epochs probe different characteristics of the SN explosion. As will be described in the next section, observations of SNe Ia during both epochs have led to considerable advances in the understanding of the SN event. This work concentrates upon late emission. It follows a previous work (Milne, The & Leising 1999) which also investigated late emission, fitting model-generated energy deposition rates to photometry of 10 SNe Ia. That work treated the 50d -200d time period in an approximate fashion, concentrating upon the 200d \- 1000d emission and will be referenced for discussions relating to gamma-ray, X-ray and positron transport and the yield of escaping positrons. This work performs similar comparisons, but concentrates more on the 50d -200d period, which spans the transition from diffusion-dominated emission to deposition-dominated emission.333In this work, all SNe Ia will be assumed to rise to peak luminosity in 18 days. In addition, this paper treats the information that is available from the R and I band observations of SNe Ia, as well as from late spectra. We compare model-generated energy deposition rates with B,V,R, and I band photometry for a large collection of well-observed SNe Ia, before and after creating bolometric corrections from a collection of late SN Ia spectra. Through these comparisons, we address four specific questions: 1) Is there any order to the late light curves of SNe Ia? 2) Is there observational support of the suggestion that positrons escape the SN ejecta? 3) Do the late light curves afford any insight into the correct explosion scenario(s)? 4) Are the late light curves of SNe Ia suitably explained by the decay products of the 56Ni $``$ 56Co $``$ 56Fe without the inclusion of additional sources of energy deposition.

## 2 Physics of SN Ia Emission

### 2.1 Emission during the first 50 days

SN light curves peak in brightness and fade on time-scales related to the lifetimes of 56Ni ($`\tau `$=8.8d) and 56Co ($`\tau `$=111d). By 100 days after the explosion, SN Ia emission has faded to less than 4% of its maximum luminosity. Largely because of the dimming of the emission, SNe Ia are better studied during the first 50 days than during later epochs. As a result, most investigations address issues relating to the diffusion-dominated epoch. Once thought to be a homogeneous class, more and better observations have revealed that intrinsic differences exist among SNe Ia (Phillips et al. 1987, Filippenko et al. 1992, 1992b, Suntzeff 1996). These differences are evident in the time evolution of the spectra, the shape of the light curves, and in the absolute magnitude of the emission at peak luminosity. Inhomogeneity at the level observed has necessitated the enlargement of the paradigm from a single, standard model to families of models.

The level of understanding of this inhomogeneity has progressed on many fronts over the last decade. Observationally, “twins” have been found for both anomalously bright, and anomalously faint SNe Ia, solidifying the existence of a large range (and perhaps a continuum) of phenomena.444Li et al. 2000 concluded that only 64% of all SNe Ia are normally-luminous, with $``$16% sub-luminous and $``$20% super-luminous. The very existence of inhomogeneity has altered the investigations of explosion scenarios. Four explosion scenarios dominate the present paradigm. The first consists of a carbon-oxygen white dwarf near the Chandrasekhar mass accreting hydrogen or helium from a binary companion until it reaches a mass at which the core carbon ignites. If the resulting burning front accelerates to become a detonation in the outer layers of the WD, a “delayed detonation” results. If the burning front remains sub-sonic, the result is a “deflagration”. These Chandrasekhar mass (CM) scenarios account for inhomogeneity by variations in the propagation of the burning front due to density and/or compositional differences in the progenitor carbon-oxygen WD. The second scenario consists of a lower mass carbon-oxygen WD accreting a helium shell, which becomes thick enough to produce a helium shell detonation. This, in turn, triggers central carbon ignition. In this sub-Chandrasekhar (SC) mass scenario, inhomogeneity is due to the different nucleosynthesis that results when the progenitor mass varies from 0.65 - 1.1 M. The third scenario merges two carbon-oxygen white dwarfs, with the more massive white dwarf accreting the companion. As with the Chandrasekhar mass scenario, central carbon ignition results. In this case, the accreted envelope is carbon. The range of masses of these explosions have been suggested to vary from 1.2 - 1.8 M, with a roughly constant nickel yield. This scenario is also referred to as a “double degenerate” explosion. The fourth scenario is not a thermonuclear explosion at all, but an accretion-induced collapse (AIC) of a white dwarf (either CO or ONeMg). It has been argued that in some cases, electron capture within an accreting white dwarf may lead to a collapse rather than central carbon ignition (as with Types II/Ib/Ic SNe). These events eject less total mass and less nickel than the other scenarios, and have been suggested to explain sub-luminous SNe Ia. Families of SN models have been developed within these scenarios (including variations of these scenarios such as pulsed delayed detonations). Using these models, many authors have demonstrated that compositional and kinematic differences within families of these SN models can roughly simulate spectral variations in observed SNe. These differences are seen and simulated near peak luminosity (Mazzali et al. 1992, 1995, 1997, Hőflich et al. 1996, Jeffery et al. 1992, Baron et al. 1996), and during the later, nebular epoch (Ruiz- Lapuente et al. 1992, 1996, Liu et al. 1997a, 1997b, 1998, Bowers et al. 1997).

The existence of inhomogeneity affects the use of SNe Ia as distance indicators at the level currently employed. There is an apparent absolute luminosity/peak-width relationship inferred from distance estimates to the host galaxies of nearby SNe Ia, quantified variously as $`\mathrm{\Delta }`$m15(B) (Phillips 1993), MLCS (Riess et al. 1996), and stretch (s) (Perlmutter et al. 1997). The ability to theoretically explain the luminosity/peak-width relationship is critical for relating well-observed, nearby SNe to SNe at the cosmological distances at which SNe Ia are used as distance indicators. This task is difficult because calculating the transport of optical photons through SN ejecta is a very complicated procedure. The ejecta is constantly evolving in both density and temperature. A photon crossing this non-equilibrated ejecta is red- (or blue-) shifted relative to local matter, wreaking havoc with the line-dominated opacity. Further complicating matters are the uncertainties of the cross sections of many relevant interactions. Despite these difficulties, detailed studies have been performed. Hőflich (1995) fit the B,V,R,I band observations of the well-observed SN 1994D with a variety of Chandrasekhar mass models, demonstrating a high level of discrimination between models. Hőflich & Khokhlov (1996) then generated B-M band light curves for Chandrasekhar mass and sub-Chandrasekhar models, and fit B,V,R,I band observations of 26 SNe Ia. One conclusion from that study is that sub-Chandrasekhar models appear too blue at peak to explain observations of sub-luminous SNe Ia. The same conclusion was reached by Nugent et al. (1997). Pinto & Eastman (1996, 2000a,b,c) investigated the influence of progenitor mass, nickel mass, nickel distribution, explosion energy and opacity upon the bolometric light curves. They then concentrated upon the ability of Chandrasekhar mass models to reproduce the B,V,R luminosity/peak-width relation in the range 0.85 $``$ $`\mathrm{\Delta }`$m15(B) $``$ 1.75. A significant achievement of those efforts is the agreement between their calculated 18d -20d rise-times of the B and V peak luminosity, with the observed rise-times of SNe Ia (Riess et al. 1999a,b; Aldering, Knop & Nugent 2000).555The models shown in Hőflich (1995) and Hőflich & Khokhlov (1996) rise to V peak in 13d -15d. Recent simulations by Hőflich et al. (1998) show that the V peak can be delayed by $``$3 days depending upon the C/O ratio, but fits of these models to observed SNe (equivalent to the previous works) have not been published.

Collectively, these investigations have reached the loose consensus that the Chandrasekhar mass explosion scenario is the favored scenario to account for the range of optical observations. Whether any subset of SN Ia events occur as sub-Chandrasekhar, merger or AIC explosions remains unclear. In this work, we simulate light curves for models representative of all four scenarios, and demonstrate general results of late light curve studies.

### 2.2 Emission after 50 days

Neither Hőflich & Khokhlov (1996) nor Pinto & Eastman carried their photometric simulations beyond 120d, both groups suggesting that their simulations become inadequate at late epochs. As described in MTL, during the time interval of interest in this work, the energy deposition is dominated by interactions involving the decay products of the 56Co $``$ 56Fe decay; gamma-ray photons and positrons. The photons possess $``$30 times more energy per decay than the positrons, but are more penetrating. This leads to the photons initially dominating the energy deposition but transitioning at later times to positron dominance. The photon transport was performed with a Monte Carlo algorithm adopting the prescription of Podznyakov, Sobol, & Sunyaev (1983). A detailed description of the Monte Carlo algorithm is given in The et al. (1990), with its application to SNe Ia light curves in Burrows & The (1990).666The resulting gamma-ray light curves agree with Hőflich (1995) and Hőflich et al. (1996) to within 5%. The dominant interaction is Compton scattering, which produces energetic electrons. These secondary electrons have a mean energy typically 300 keV at 100d. MTL assigned these electrons zero lifetimes, depositing the energy in-situ. This work improves on that by including the transport of the secondary electrons with the same algorithm employed for the positron transport, but using Mőller scattering rather than Bhabha scattering. The effect of the secondary-electron transport is small but non-negligible; in the Chandrasekhar mass models with radial positron escape and low ionization, the light curves are fainter by $``$ 0.10m; for low-mass sub-Chandrasekhar models, the effect reaches 0.18m.

Positron transport depends upon the nature of the magnetic field. Three scenarios have been suggested to model the magnetic field (Ruiz-Lapuente & Spruit 1998, hereafter RLS). The first suggests that the field is too weak to confine positrons, and positrons follow straight line trajectories, with a fraction escaping the ejecta (Colgate, Petschek & Kriese 1980). The second suggests a stronger field that confines positrons, but with the field lines radially-combed by the homologous expansion. The positrons spiral along these radial field lines with their pitch-angles decreasing due to the field gradient (beaming), with a fraction escaping the ejecta (Chan & Lingenfelter 1993). The third scenario suggests a strong field that is turbulently disordered such that positrons mirror frequently with no net transport (Axelrod 1980). Positrons may survive non-thermally at late times, but none escape. Colgate, et al. (1980) argued that the first two situations are equivalent. Simulations by MTL found this to be approximately true for most SN models. The first two scenarios will herein be referred to as the “radial” scenario, the third will be referred to as the “trapping” scenario. The differing positron escape leads to differences in model energy deposition rates, which makes late observations of SNe Ia a probe of the photon and positron transport.

The positron and secondary electron transport was performed with a Monte Carlo algorithm, as explained in MTL. The dominant energy loss mechanisms are ionization and excitation of bound electrons for low levels of ionization, and plasma excitation for higher ionization. The model-generated energy deposition rate for the sub-Chandrasekhar model, HED8, is shown in Figure 1. The dashed line (D) assumes instantaneous deposition of all decay energy, an assumption that is invalid due to gamma-ray photon escape by 30d. The dotted line (G) uses the results of the gamma-ray/secondary electron deposition only and assumes no deposition of positron kinetic energy. By 100d, the deposition of positron kinetic energy is an important contributor to the total energy deposition rate, but the positron lifetimes are short enough that all positron curves approximate instantaneous, in-situ deposition of positron energy (In). At later times, densities in the ejecta are low enough to allow escape in the radial scenario (dark shading, R) and non-thermal survival in the trapping scenario (light shading, T). Both scenarios were calculated for a range of ionizations, from 1% of all nuclei being singly ionized, to all nuclei being triply ionized. Positron escape removes energy from the ejecta, leading to fainter light curves than in-situ deposition. Lower ionization permits more positron escape and thus is the fainter edge of the radial curve. A trapping field permits large non-thermal lifetimes, but no escape. This leads to storage of energy and a late light curve that is brighter than in-situ deposition. For this scenario, low ionization leads to a larger energy storage and the brighter edge of the trapping curve.777The normally-luminous sub-Chandrasekhar model shown allows a relatively large escape fraction for the radial scenario compared to normally-luminous Chandrasekhar mass models. Nonetheless, the basic features shown in Figure 1 are characteristic of all the SN Ia models simulated.

### 2.3 Model/Observation Comparisons

The ideal model/observation comparison would entail an explicit NLTE determination of photon and positron energy deposition rates and from these rates the generation of optical/IR spectra. These spectra would then be compared with a sequence of observed SN Ia spectra, comparing the evolution of the total flux from early to very late epochs. This ideal situation is far from being realized, both computationally and observationally. Falling short of that solution, an alternative approach has been employed in numerous works; fitting energy deposition rates to multi-band photometry or estimates of the uvoir bolometric luminosity. Two issues must be addressed for these comparisons to be meaningful, photon diffusion delay and color evolution. Photon diffusion delay refers to the portion of the time delay between the energy deposition (from the scattering of gamma- and X-ray photons and slowing of positrons) and the production of the optical light that can be seen by an external observer that is due to the diffusion of optical photons out of the ejecta. This diffusion delay determines the shape of the SN light curve peak at early times, but decreases to negligible values at late times. Thus, at late epochs, energy deposition is instantaneously recycled into optical emission. In previous studies of late light curves of SNe Ia, two approaches have been used to account for photon diffusion delay. Studies that attempt to fit observations continuously from early to late epochs derive this time delay as a function of explosion epoch. The alternative approach is to study only the late epochs and assume instantaneous recycling. Observationally, the late epoch begins when the emitted spectrum transitions from continuum emission to (forbidden) nebular lines. Color evolution refers to the changing fraction of the total optical flux that a given photometric band samples due to the evolution of the spectrum. Again, two different approaches have been employed. The first approach creates bolometric light curves from the photometry (applying weighting factors derived from spectra) and compares these light curves to model-generated energy deposition rates.888Contardo et al. (2000) discuss the issues related to generating bolometric light curves from band photometry. The second approach determines the epoch after which a given photometric band scales with the total optical flux (the epoch is determined from spectra), and compares this band photometry to model-generated energy deposition rates.

Two groups that carried their investigations to late epochs have treated photon diffusion in their comparisons. Colgate, Petschek and Kriese (1980) fit Monte Carlo model simulations of the energy deposition rate to B band data from SNe 1937C and 1972E. That work assumed a grey opacity for both photon and positron transport and a single-zone model. They concluded that positron escape was required to explain the late light curves. Cappellaro et al. (1998) (hereafter CAPP) and Salvo et al. (2001) transported gamma-ray photons, positrons and optical photons through W7-like models (scaled to various masses) with a Monte Carlo code and fit the energy deposition rates to V band data from SNe 1991T, 1994D, 1992A, 1993L, 1996X and 1991bg.999Both Colgate et al. (1980) and CAPP transported the positrons with routines that were produced for photon transport. The positrons were given a speed $`c`$. That work assumed a grey opacity for photons and varied the positron opacity to fit the late light curves. They concluded that positron escape was suggested for some, but not all SNe.

RLS transported gamma-ray photons and positrons through Chandrasekhar mass and sub-Chandrasekhar models without treating photon diffusion delays, and fit the energy deposition rates to bolometric light curves for the SNe 1972E, 1992A and 1991bg.101010The bolometric light curves were generated from photometry for SNe 1992A and 1991bg. For SN 1972E, the bolometric light curve was generated from a series of optical spectra. That work concentrated upon very late times ($``$ 100d) at which time the photon diffusion delay time-scale is assumed to be short. They concluded that normally-luminous SNe Ia transition from a positron-trapping phase to a positron-escape phase after 450d, and can be fit with Chandrasekhar mass models. They further concluded that sub-Chandrasekhar models best describe the sub-luminous SNe Ia and feature positron escape as early as 150d. MTL transported gamma-ray photons and positrons through Chandrasekhar mass and sub-Chandrasekhar models, and fit the energy deposition rates to bolometric, V band and B band light curves for 10 SNe, including the SNe studied by the previously mentioned groups. That work also concentrated on very late times and did not treat photon diffusion delays. They concluded that both Chandrasekhar mass and sub-Chandrasekhar models can explain normally- and super- luminous SNe Ia, with positron escape suggested to be consistent with seven of the eight SNe. They further concluded that none of the models tested could explain the sub-luminous SN 1991bg, a result in agreement with CAPP, but in disagreement with RLS.

This work does not simulate photon diffusion delays, comparing model-generated energy deposition rates to multi-band photometry as though the delays were negligible. This approach uses the comparisons to determine the onset of instantaneous recycling. Similarly, the initial comparisons of the model-generated energy deposition rates to each of the B,V,R and I band data-sets assume no color evolution. The differences between these comparisons and the later comparisons that incorporate spectrally-derived bolometric corrections afford a measure of the importance of color evolution at late epochs.

## 3 SN Ia Observations

To perform model-generated energy deposition rate/SN light curve comparisons, a large collection of SNe Ia observations have been compiled. The principle sources of B,V,R,I band photometry are the 29 SNe Ia observed by Hamuy et al. (1996), and the 22 SNe Ia observed by Riess et al. (1999). MTL fit the 10 SNe Ia best observed to very late times. Of those ten SNe, only six are included in this study (as explained below). These SNe have been supplemented with observations of SN 1995D by Sadakane et al. (1996), SN 1996X by Salvo et al. (2001), SN 1997cn by Turatto et al. (1998), SN 1998bu by Jha et al. (1999), Suntzeff et al. (1999), Garnavich et al. (2000), Cappellaro et al. (2001), and SN 1998de by Modjaz et al. (2000). The SNe Ia have been divided into luminosity sub-classes (normally-, sub- and super-luminous). Hereafter these categories will be denoted in figures by N, sb & SP. The categorization is based upon the $`\mathrm{\Delta }`$m15(B) (Phillips 1993) values of the SNe. Numerous empirical relations have been derived for SNe Ia, plotting a given variable versus the $`\mathrm{\Delta }`$m15(B) parameter. Shown in Figure 2 are the maximum B,V and I band absolute magnitudes for 21 of the 22 SNe Ia analyzed in this study combined with 18 others from Hamuy et al. 1996.111111This figure only included the SNe with well-defined peak shapes. Early studies assumed a linear relationship between these variables for all SNe Ia. A larger sampling has suggested that all these relations deviate from linearity at $`\mathrm{\Delta }`$m15(B) $``$ 1.6. It is unclear whether this deviation is best explained by a continuous quadratic function applicable to all SNe Ia (Phillips 1999), or by two different linear relations, one for each sub-class of SNe Ia. Regardless of the interpretation of this break, the normally-/sub-luminous cut has been set at 1.60. By that definition, SN 1993H is considered sub-luminous and SN 1992A is considered normally-luminous. The normally-/super-luminous cut ($`\mathrm{\Delta }`$m15(B)=0.95) was chosen such that the well-observed SN 1991T (which showed spectral deviations from normal SNe Ia) was considered super-luminous. It is important to note that the absence of a template SN with a wider peak than SN 1991T hinders the accurate estimation of $`\mathrm{\Delta }`$m15(B) for super-luminous SNe. As will be shown, the exact location of the normally-/super-luminous cut is unimportant for this work, but the normally-/sub-luminous cut does differentiate two distinct late light curve sub-classes.

From 64 observed SNe Ia, 22 were included in this study. To be included, a SN had to be observed at least once after 85d with at least one observation within $`\pm `$20 of the normalization epoch (explained below), chosen to be 65d post-explosion. Additionally, the SN had to be first observed no later than one week post-B maximum (a requirement relaxed for the SN 1992K due to the undersampling of the sb sub-class). The SNe Ia which meet this criteria are; N:90N,90O,92A,92al,93ag,94D,95D,95E,95al,96X,98bu, SP:91T,91ag,92bc,94ae,95bd, sb:86G,91bg,92K,93H,97cn,98de. Four SNe Ia from MTL were not included in this study. SNe 1937C and 1972E were not well-observed with multi-band photometry and were excluded. SN 1993L was discovered more than a week post-peak and was excluded. The late photometry of SN 1989B is being re-analyzed and thus that SN was excluded.121212The SN has been identified as a light echo candidate (Boffi et al. 1999). As published, the late light curve of SN 1989B is more than 1m brighter than the suggested models and is the singular exception to the trends seen in the other SNe. All three sub-classes are better sampled from 50d -150d than after 150d.

The phenomenon studied in this work is the evolution of the light curves at late epochs. For this reason, a relative magnitude scale has been employed rather than an absolute or apparent magnitude scale. All SNe have been normalized to have zero magnitude at 65d post-explosion (assuming an 18d rise-time for all sub-classes). The initial motivation for the 65d normalization date was based upon the suggestion in Pinto & Eastman (2000) that the photon diffusion time-scale becomes negligibly short at roughly that time. As the next section shows, a high level of homogeneity exists within the multi-band photometry of each sub-class after 65d. Additionally, the greater number of 45d -85d observations relative to later epochs makes the normalization algorithm perform better at early times. These two factors led us to retain the 65d normalization date for all data-sets, although the model-generated energy deposition rates do not fit the band photometry at 65d for any band/sub-class combination. The normalization has been performed by linear interpolation of all data between 45d and 85d. For two SNe, only a single observation existed in that interval. For those SNe, the latest pre-45d observation was used to perform the interpolation.

By treating all SNe within a given $`\mathrm{\Delta }`$m15(B) range as a single object, we are assuming homogeneity within each luminosity sub-class. The amount of scatter about a template gives some measure of the quality of that assumption, but numerous sources of systematic errors exist to cloud that interpretation. Each SN light curve is measured relative to the host galaxy’s background light, which is unique to that SN. If the SN is near the detectability limit for that observation, the influence of the background subtraction increases in importance. The irregularity of the sampling means that highly uncertain observations are mixed throughout every epoch of the late light curves. Similarly, different telescopes, filters and detectors were used to produce this data-set, introducing additional systematic differences between individual data-points. As was shown by Suntzeff (1999), observations of the same SN by the same observer at the same site with the same class telescope can yield different photometry. The differences were found to be in excess of the quoted uncertainties. Suntzeff showed that in this optimal case, the “instrumental uncertainty” was on the order of 0.06m for the V band during the 60d -75d epoch. For the case of this mixed data-set, the uncertainties are likely to be larger. Errors in this work were also introduced by the fitting algorithm, which fit the individual light curves to the estimated 65d magnitude.

All errors discussed above have been related to estimates of the relative magnitude. As SNe are always detected at some time after the explosion, the explosion date is also somewhat uncertain. Two factors may influence uncertainties in the explosion date; 1) uncertainty of a sub-class rise-time to peak, and 2) uncertainty of the peak date. For this work, all SNe were assumed to rise to peak B magnitude in 18d regardless of the luminosity sub-class. To minimize the effect of peak date uncertainty, only SNe first observed no later than a week after maximum light were used (with the exception of SN 1992K). The determination of the peak date depends upon template fitting, so peak date uncertainties are only relative to the template used. SNe discovered before tB(max) have peak date uncertainties typically of $`\pm `$1d. SNe discovered post-peak have larger uncertainties. Nonetheless, as the time-scales involved in this work are on the order of tens of days, uncertainties on the order of a few days will not significantly alter these results.

## 4 Results

### 4.1 V Band

The 22 SNe included in this study are shown fit (by eye) with the delayed-detonation, Chandrasekhar mass model DD23C in Figure 3. All 22 SNe have been normalized at 65d. Three tendencies are apparent in the data. First, within each sub-class the data show remarkable homogeneity. The sub-classes are DEFINED by differences in the early light curve shape, but after 65d, there is no evidence of variations. Second, the normally-luminous and super-luminous data appear to evolve similarly, while the sub-luminous sub-class continues to fall more steeply (recall that a steep decline from peak DEFINES this sub-class). Third, the separation between the normally-/super-luminous and the sub-luminous sub-classes far exceeds the scatter within the sub-classes. This argues against a continuous transition between the sub-luminous and normally-luminous sub-classes. Unfortunately, the sub-luminous sub-class is under-sampled, so many more sub-luminous SNe need to be observed before addressing whether there exists an absolute separation of sub-classes.

After $``$80d the normally-luminous and super-luminous sub-classes appear to be explainable by this model, assuming radial escape. By contrast, the sub-luminous sub-class cannot be explained by either magnetic field scenario. To improve the visualization of the 40d -120d normally-luminous and super-luminous data, in the inset of Figure 3 we show the data as residuals relative to the model curve. The similarities between the normally-luminous and super-luminous data-sets is apparent. Also apparent in the inset are the failure of the model light curves to fit the V band data in detail until after $``$120d. These results are not specific to the delayed-detonation, Chandrasekhar mass model, DD23C. Shown in Figure 4 are six models suggested to explain normally-luminous and super-luminous SNe, shown in the “delta magnitude” format.131313The delta magnitude format shows the residuals of the data and model-generated light curves to the instantaneous deposition approximation, in units of magnitude. The models are defined in Table 1. All model light curves show the same fundamental structure, with the models fitting the data after $``$170d and suggesting positron escape. The Chandrasekhar mass models fit the normally-luminous data better than does the sub-Chandrasekhar mass model, but the improvement is modest. As was discussed in MTL, the super-luminous model light curves are similar, both to other super-luminous SN models and to normally-luminous Chandrasekhar mass model light curves. We assert that the 16 normally- and super-luminous SNe do not need to be differentiated into distinct sub-classes after 60d. We further assert that if the V band emission scales with the bolometric emission to the level of a few tenths of a magnitude, then the late emission from all 16 SNe can be suitably explained with a single light curve featuring positron escape, and that Chandrasekhar mass, sub-Chandrasekhar mass, and merger explosions can all explain these data.

It is apparent from Figure 3 that the sub-luminous sub-class differs considerably from the normally-/super-luminous sub-classes. MTL concluded that none of the models they tested could explain the 60d -560d evolution of the light curves. One interpretation of that result is a rejection of both the Chandrasekhar mass and the sub-Chandrasekhar mass models tested. A different interpretation would be the failure of the assumption that the energy deposition rate scales with the V band during that epoch. Figure 5 shows five sub-luminous models fit to the sub-luminous sub-class data. The models are normalized to fit the data at a later epoch ($``$170d). It is clear from the figure that none of the models can explain the V band data before 150d.141414The fifth panel (lowest) is an accretion-induced collapse model, which ejects 0.2 M of material. The model was created by scaling ONeMg to 0.2 M. This model is included upon the suggestion of Fryer et al. (1996,1999) that AIC’s may explain sub-luminous SNe Ia. Even at this extremely low mass, the early escape of positrons is not enough to reproduce the steepness of the 65d -170d light curve. After 150d, the V band data can be suitably explained by all models if positron escape is assumed.151515 It is important to note that the range of ionizations allowed lead to fairly thick curves for low-mass models, which exaggerates the visual impression of their ability to explain the data. A treatment of the level of ionization would reduce the range to a single curve. This latter interpretation implies that the energy deposition rate scaling with the V band occurs later for the sub-luminous sub-class than for the other sub-classes. For this interpretation to be correct, there must be an explanation which is consistent with the spectral observations, which suggest that sub-luminous SNe Ia enter the nebular phase earlier than normally-luminous or super-luminous SNe Ia (Mazzali et al. 1997).

### 4.2 BVRI Bands

In the previous section, we fit model-generated energy deposition rates to V band data. That procedure relies upon the assertion that the V band traces the bolometric luminosity. For the normally-/super-luminous data-set, the data are approximately fit with the models, suggesting a validation of that assertion. By contrast, for the sub-luminous data-set the models could not fit the data until a later epoch. All SNe Ia used in this study were also observed in the B, R and I bands. By independently fitting these data-sets to the model-generated energy deposition rates, we further explore the nature of the late-time emission from SNe Ia.

Shown in Figure 6 are the BVRI data-sets for SNe Ia fit to the model, W7. Within each sub-class the data has been normalized to 65d. Each sub-class has then been independently fit to the model light curves. For all four photometric bands, the normally-luminous and super-luminous light curves are similar. The B band data for the normally-luminous and super-luminous SNe Ia deviate considerably from the model before 100d, but by $``$170d they roughly scale with the energy deposition rate. The B band data for the sub-luminous sub-class differ from the normally-luminous and super-luminous sub-classes early, but appear similar at late times. Perhaps the most striking results are in the R and I bands, where the three sub-classes are quite similar. We argue that the differences between the late evolution of the B and V band data do not carry over to the longer wavelengths. The model fits initially deviate from the R and I band data for all three sub-classes, but then converge upon the model curve by $``$170d. The convergence in the I band is less certain due to scatter.

It is interesting to note that the onset of the energy deposition scaling with the photometry ($``$170d) occurs later than the transition from photospheric continuum emission to nebular emission ($``$60d) for all three SN Ia sub-classes and all four photometric bands. To determine whether this delay is explainable as a consequence of color evolution, we produce spectral templates for normally-/super-luminous and sub-luminous SNe Ia. The templates have been created from a collection of SN Ia spectra; 1994D (Patat et al. 1996; Turatto 2000; Filippenko 1997), 1991T (Filippenko et al. 1992), 1996X (Salvo et al. 2001), 1991bg (Turatto et al. 1996; Ruiz -Lapuente et al. 1993), 1987L (Ruiz -Lapuente et al. 1993), 1994ae (Bowers et al. 1997), 1998de (Modjaz et al. 2000), 1998bu (Jha et al. 1999; Matheson 2000), 1995D (Matheson 2000), 1981B (Branch et al. 1983), 1984A (Branch 1987), 1991F (Gomez & Lopez 1995), 1992A (Kirshner et al. 1992). We include all emission within the 3500Å-9700Å wavelength range, and ignore all emission outside of this range. Spectra that do not span that range have been linearly interpolated from earlier and later spectra.161616The interpolation was performed in two steps. The shape of the interpolated portion was determined from the two complete spectra adjacent in time. That shape was then spliced to the existing portion. The spectral sequences for a sampling of each sub-class are shown in Figures 7 and 8. From these templates, we have determined the evolution of the “filter efficiency”, or the fraction of the energy emission in the 3500Å-9700Å wavelength range that is detected by each band. The filter efficiencies for normally-luminous and super-luminous SNe Ia, shown in Figure 9, match the shape of the residuals of the data-model fits shown in Figure 6 for the normally-/super-luminous data-set. The V band filter efficiency remains nearly constant, in agreement with previous assertions that the V band scales with the bolometric luminosity.171717Contardo et al. (2000) arrived at a similar result. Shown in Figure 10 are the B,V,R and I band data fit to the model W7, after deriving bolometric corrections from these filter efficiencies. The data for all four bands are suitably explained by the energy deposition rates of the model if positron escape is assumed. Before $``$60d, the B,V and I data vary dramatically from the energy deposition rates. This would be consistent with the decrease to negligible values of the photon diffusion time-scale, as suggested by Pinto & Eastman (2000a). The correction lessened, but did not completely remove the residuals from the B and V band data during the 60d -170d epoch. It is unclear whether the residuals remain due to a modest failure of this relatively crude bolometric correction calculation, or rather they remain because this treatment has ignored some physical phenomenon.

The application of these corrections essentially creates bolometric light curves from each photometric band. A single bolometric light curve for normally- and super-luminous SNe Ia was generated by combining the four curves shown in Figure 10 after weighting each curve according to the photometry errors. The resulting bolometric light curve can be compared with model-generated light curves. These comparisons, shown in Figure 11, exhibit the same tendencies that are seen in the V band; the energy deposition rates for all four models are consistent with the data, and positron escape is suggested. These curves differ slightly from the bolometric light curve derived by Contardo et al. (2000) from much of the same data. That study was primarily concerned with emission near the luminosity peak and assumed constant bolometric corrections at late times ($``$ 130d). As we have shown that some evolution occurs during the 60d \- 200d time period, we assert that this bolometric light curve is more accurate at late times than the light curve from Contardo et al. (2000). Although bolometric light curves include information from four photometric bands rather than simply from the V band, the V band is the best observed. Thus, we present the bolometric light curve primarily to demonstrate that conclusions derived from studies of the V band shape are valid.

There are few nebular spectra for sub-luminous SNe Ia, so for this sub-class, the filter efficiencies are very crude and are undersampled after 150d (Figure 12). This is problematic, because as shown in Figure 6, the 100d -200d time-span is a critical transitional epoch. Nonetheless, as shown in Figure 13, the B,V,R and I band data are also fairly sparse and can be roughly explained with the model PDD54, after correcting for filter efficiency. All four bands fit the energy deposition rates better with the correction, as seen through comparison with Figure 6. There are large residuals in the V and I bands before $``$50d. This is likely due to photon diffusion delays. It is evident in Figure 13 that the light curves of sub-luminous SNe Ia are consistent with model-generated energy deposition rates which allow positron escape. More spectral and photometric data of sub-luminous SNe Ia need to be collected to further the understanding of this epoch. Due to the crudeness of the bolometric corrections, a bolometric light curve was not generated for the sub-luminous sub-class.

The success of the bolometric corrections in reducing the residuals during the 60d -170d time-span suggests the existence of an epoch during which optical emission is instantaneously recycled, but the SN color continues to evolve. A possible explanation for this color evolution is the decay of 56Co to 56Fe. It is during the 65d -170d epoch that the daughter becomes the dominant species in that decay. As shown by both Liu et al. (1997a,c) and Bowers et al. (1998), cobalt emission dominates the wavelength ranges observed with the R and I bands. As 56Co decays to 56Fe, these bands would fade faster than the energy deposition rate due to color evolution until late-times at which time most cobalt has decayed and iron and stable nickel dominate the spectrum in that wavelength range. In the V band wavelength range, both Co and Fe lines are present and the bolometric correction would change only slightly during the transition. The B band may brighten due to the increased emission from the many \[FeII\] lines from 446nm -456nm. This explanation can be tested through comparisons of model-generated spectra with observations made during this epoch.

## 5 Discussion

Two tasks comprised this work. The first task was the compilation of a data-set of B,V,R,I photometry of type Ia SNe divided into sub-classes, and development of bolometric correction factors for each band. The second task was the generation of late-time light curves from SN Ia models and the fitting of these models to the compiled data-set. Certain findings can be derived from this study independent of the model-fitting, and are worthy of mention. We have demonstrated that there is order to the late light curves. The late B and V band light curves of sixteen normally- and super-luminous SNe follow a similar evolution, while those of six sub-luminous SNe follow a different evolution. By contrast, the late R and I band light curves of all 22 SNe Ia follow a similar evolution. The normally-/super-luminous sub-classes appear virtually identical although they span more than half of the $`\mathrm{\Delta }`$m15(B) range of SNe in this study. The B and V band light curves of the sub-luminous SNe appear to require a distinct sub-class. Whether there are in fact SNe Ia that connect them to normally-luminous objects is an important question. This question can only be answered by observational campaigns that focus upon monitoring SNe that span the range of $`\mathrm{\Delta }`$m15(B) values.

Comparisons of model-generated energy deposition rates with photometry in the B,V,R, and I bands reveal three epochs. There is an early epoch ($``$ 60d), where energy deposition/photometry comparisons are invalid presumably due to the time delay between the deposition of energy and the emergence of the resultant optical emission. That epoch is followed by an intermediate epoch (60d -170d) where the bolometric light curve is fit by the energy deposition rate, but color evolution must be addressed for individual photometric bands to be fit by the energy deposition rate. This epoch is followed by a late epoch ($``$ 170d) where the bolometric corrections are roughly constant and each individual photometric band can be fit by the energy deposition rate if positron escape is allowed for. These epochs exist for both the normally-/super-luminous sub-class, and the sub-luminous sub-class. The bolometric corrections during the intermediate epoch are better defined for the normally-/super-luminous SNe Ia than for the sub-luminous SNe Ia, the latter corrections are preliminary.

Collectively, the light curves of SNe Ia after 60d suggest that the interactions of the products of the 56Co $``$ 56Fe decay with the ejecta can explain the energy deposition without any additional energy deposition source. Positrons are seen to escape the ejecta in quantity for all SNe Ia during the late epoch. These findings hold for Chandrasekhar mass, sub-Chandrasekhar mass, merger and AIC models equally. The model independence of these results does not imply that the late emission from SNe Ia cannot probe the progenitor. The challenge for NLTE radiation transport calculations is to reproduce the spectra during this epoch.

There are two ramifications of positron escape from SN Ia ejecta that warrant further discussion. First, nebular spectra have been used to estimate the 56Ni production in SNe Ia (Ruiz-Lapuente & Filippenko 1996, Bowers et al. 1997) and to differentiate between Chandrasekhar mass and sub-Chandrasekhar mass models (Liu et al. 1997a,1997b,1998). In all of these studies, instantaneous positron energy loss was assumed. By 300d, the energy deposition is dominated by positron slowing, and a substantial fraction of those positrons are expected to escape the ejecta in the radial scenario. Even positrons that do not escape the ejecta diffuse from the location of their creation. In light of this and other recent studies, nebular spectrum studies should be calculated with realistic positron transport. Second, positrons that escape the SN Ia ejecta are thought to survive non-thermally on time-scales of 105 years, or more (Chan & Lingenfelter 1993, Guessom, Ramaty & Lingenfelter 1991). The collective positron contributions from SNe Ia, as inferred from extra-galactic SN rates, are sufficient to generate a large fraction of the 511 keV positron annihilation radiation observed by the CGRO/OSSE, SMM and TGRS gamma-ray telescopes (Milne et al. 1999,2000). The degree of dominance is enough that accurate characterization of positron annihilation radiation will trace recent SN Ia activity in the Galaxy.

Light echoes have been detected from two SNe Ia, 1991T and 1998bu. In addition, Boffi et al. (1998) have identified SN 1989B as a light echo candidate. Of all SNe Ia observed after 300d, these three are the only SNe to remain significantly brighter than those in this work. As light echoes are currently being used as geometric distance indicators, this work suggests that the templates shown in this work may be used to locate light echoes (Sparks et al. 1999).

There remain incomplete aspects to this study. More SNe Ia need to be monitored to reduce gaps in observations and determine the late behavior of more SNe Ia along the $`\mathrm{\Delta }`$m15(B) sequence. The exact nature of the transition from the onset of instantaneous recycling to the cessation of color evolution must be better observed spectrally and in each photometric band, particularly with sub-luminous SNe Ia. These unexplained aspects underscore the fact that much is to be learned about the physics of Type Ia SNe from late observations.

We thank P. Hőflich for ongoing access to SN models and for discussion relating to radiation transport through SN ejecta. We thank M. Turratto and S. Jha for access to SN photometry and spectra. We also thank T. Matheson, R. Lopez, and P. Meikle for SN spectra.