The long-term evolution and appearance of Type Iax postgenitor stars
DDraft version March 5, 2019
Typeset using L A TEX twocolumn style in AASTeX62
The long-term evolution and appearance of Type Iax postgenitor stars
Michael Zhang, Jim Fuller, Josiah Schwab, ∗ and Ryan J. Foley Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA Department of Astronomy and Astrophysics, University of California Santa Cruz, Santa Cruz, CA 95064, USA
ABSTRACTType Iax supernovae may arise from failed explosions of white dwarfs that leave behind a boundremnant (i.e., a “postgenitor” star) that could be identified in wide field surveys. To understandtheir observational signatures, we simulate these white dwarf (WD) postgenitors from shortly afterexplosion until they move back down the WD cooling track, and we consider several possible WDmasses and explosion energies. To predict the peculiar surface abundances of the WD postgenitors,our models take into account gravitational settling and radiative levitation. We find that radiativelevitation is significant at temperatures above a mass-dependent critical temperature, typically in therange T eff ≈ − × K, significantly increasing surface abundances of iron-group elements.Due to enhanced iron group opacity compared to normal WDs, the postgenitor peak luminosity andcooling timescale depend sensitively on mass, with more massive WDs becoming brighter but coolingmuch faster. We discuss our results in light of recently discovered hypervelocity white dwarfs withpeculiar surface compositions, finding that our low-mass postgenitor models match many of theirobservational characteristics. Finally, we explore the effects of thermohaline diffusion, tentativelyfinding that it strongly suppresses abundance enhancements created by radiative levitation, but morerealistic modeling is required to reach a firm conclusion.
Keywords: supernovae:general, ISM:supernova remnants, stars: evolution, white dwarfs INTRODUCTIONIn 2002, a very peculiar supernova (SN) was discov-ered at Palomar Observatory (Wood-Vasey et al. 2002).This supernova, soon to be named 2002cx, did not fitneatly into any category. It had the pre-maximum spec-trum similar to that of a Ia, with absorption lines fromintermediate-mass elements and iron. However, it alsohad a very low luminosity, half the typical expansionvelocity, and atypically red colors. The discovery paper(Li et al. 2003), aptly entitled ‘SN 2002cx: The MostPeculiar Known Type Ia Supernova’, concluded that noexisting model can explain the supernova.Since then, dozens of supernovae like 2002cx have beendiscovered. Foley et al. (2013) grouped these events intoa distinct class, called Type Iax, and estimate its eventrate as roughly 1/3 that of type Ia. SNe Iax are the mostcommon class of peculiar Ia-like supernovae (Jha 2017).Like their prototype, type Iax supernovae are character-
Corresponding author: Michael [email protected] ∗ Hubble Fellow ized by a Ia-like pre-maximum spectrum, but with peakluminosities typically a few magnitudes fainter and ex-pansion velocities a few times slower. SNe Iax are moreinherently diverse than SNe Ia, with no strong width-luminosity relation, peak luminosity varying over 4-5magnitudes, and velocities varying by a factor of 4. SNeIax also have other interesting properties. For example,they have never been observed in elliptical galaxies (al-though there is one example in an S0 galaxy: Foley et al.2010) and prefer star-forming spirals, yet there is alsono sign of ongoing star formation at the site of any Iax(Foley et al. 2013).Ia/Iax SNe result from the explosion of a white dwarfwith a binary companion, but neither the nature of thecompanion nor the mechanism that triggers the explo-sion is known. In the single degenerate Ia/Iax scenario,the companion is a non-degenerate star; in the doubledegenerate Ia/Iax scenario, it is another white dwarf.Wang et al. (2013) argue for a single degenerate sce-nario for SNe Iax, arguing that they could be the prod-uct of accretion from a helium star onto a CO WD, amechanism which reproduces the long delay times andluminosity diversity of SNe Iax. The detection of a a r X i v : . [ a s t r o - ph . S R ] M a r luminous blue progenitor of a SN Iax (McCully et al.2014) supports this scenario. It has also been proposedthat SNe Iax could result from hybrid CONe progenitors(Meng & Podsiadlowski 2014). Bravo et al. (2016) inves-tigates explosions of hybrid CONe WDs created by off-center carbon burning in intermediate mass stars, find-ing such explosions could leave behind bound remnantWDs. Kashyap et al. (2018) explain Ia/Iax SNe as be-ing due to the merger of a CO WD with a ONe WD.While the CO mixture burns easily, the ONe mixturedoes not, creating a low-luminosity transient with smallejecta mass.Despite these unknowns, the leading explanation forSNe Iax is that they are white dwarf (WD) deflagra-tions that do not lead to detonations (Branch et al. 2004;Jha et al. 2006; Phillips et al. 2007), which explains theIa-like spectrum, the low luminosity, and the low ex-pansion velocity. Deflagrations also tend to produce awide range of explosion energies, explaining the diversityof type Iax events, whereas Chandrasekhar-mass deto-nations are more uniform in their properties. Kromeret al. (2013) performed a 3D deflagration simulation ofa Iax explosion, successfully reproducing characteristicobservational features of SN 2005hk in the optical andnear-infrared. The asymmetric mass ejection may alsoimpart a kick of several hundred km/s to the boundremnant star (Jordan et al. 2012).If SNe Iax supernovae are truly partial deflagrationsof CO WDs or deflagrations and delayed detonations ofhybrid CONe cores, they will not be energetic enoughto unbind the white dwarf. Shen & Schwab (2017) sim-ulate bound WD postgenitors using the MESA stellarevolution code (Paxton et al. 2011). They take into ac-count delayed Ni decay and super-Eddington winds topredict the light curve from days after the explosion to1000 years afterwards. They report high uncertainties,but reasonably good matches for the late-time luminos-ity, temperature, and velocity of SN 2005hk, SN 2008A,and SN 2008ha.Recently, Vennes et al. (2017) discovered a hyperveloc-ity white dwarf (LP 40-365) with an oxygen-neon atmo-sphere and abundant intermediate-mass elements; theyconsider LP 40-365 a candidate for a Iax postgenitor.Using Gaia data, Raddi et al. (2018) confirm the hyper-velocity nature of the object while reporting a radius of0 . ± . R (cid:12) and a mass of 0 . +0 . − . M (cid:12) . Shen et al.(2018) discovered a few more peculiar hypervelocity sub-dwarf/white dwarf stars with similar properties. Theseobjects also have no detectable hydrogen or helium, butdo have strong carbon, oxygen, iron, magnesium, andcalcium features. The large space velocities, peculiarsurface compositions, and unusual masses/radii suggest these stars could be Iax postgenitors, although thoseauthors posit that they are the degenerate donor starcompanions of white dwarfs that exploded as SN Ia. Fi-nally, Kepler et al. (2016) reports an enigmatic WD witha nearly pure oxygen atmosphere. It is unclear how thisobject formed, but we speculate that it could be a ONeor CONe WD that deflagrated long ago.In this paper, we extend the postgenitor simulationsof Shen & Schwab (2017) to very late times to deter-mine the properties of these bound remnant stars. Wediscuss the setup of our simulations in Section 2, showthe most salient characteristics in Section 3, and analyzethe results while comparing them to hypervelocity starslike LP 40-365 in Section 4. MODELSTo simulate WD postgenitors, we perform stellar evo-lution calculations using MESA (Paxton et al. 2011,2013, 2015, 2018) version 10000. We set up a gridof models with different initial conditions, then evolvethem hydrostatically to predict the evolution of observ-ables such as temperature, luminosity, radius, and sur-face abundances.We do not attempt to simulate the supernova itself,nor do we take into account binary interaction, winds,radioactive decay, or detailed radiative transfer. Oursimulation is meant to start at late times (more than ∼
100 yr after explosion), when all significant radioactivenuclides have decayed, all winds have died down, and thesupernova remnant has long since become optically thin.The simulation assumes hydrostatic quasi-equilibriumat all times, taking into account convection, radiativetransport, element diffusion, radiative acceleration, andneutrino cooling to evolve the white dwarf until it is fardown its cooling track.2.1.
Initial conditions
To set up somewhat realistic initial conditions, we takea packaged white dwarf model from MESA and adjustits properties. We relax its composition to one appro-priate for carbon-oxygen Iax postgenitors (Kromer et al.2013), then relax the outer portions of the white dwarfto a constant entropy. These outer portions, which willhenceforth be called the “envelope”, represent a combi-nation of the nuclear ashes and the fall-back ejecta. Ahomogenous and constant entropy envelope is expectedin regions well mixed by convective burning, which en-forces a nearly isentropic structure, and is likely to bethe case in the thermally supported envelope of the WDsoon after deflagation.To account for the pollution with nucleosyntheticburning products of the deflagation, we relax the el-emental abundances in the stellar envelope to thoseshown in Table 1. These numbers are taken fromKromer et al. (2013), which describes 3D deflagrationsimulations of Iax supernovae and predicts final elemen-tal abundances of both the postgenitor and the ejecta.The precise numbers are not very important for our pur-poses, as we are more interested in the evolution of sur-face abundances over time due to diffusion and radiativelevitation. However, we note that the large iron-groupabundance, roughly 30x solar, is important for the evo-lution of the WD and its spectroscopic appearance.To account both for the inherent diversity of Iax su-pernovae and for uncertainty in the explosion processand its outcome, we use a grid of 24 initial conditions.The grid contains 4 postgenitor masses (0.15, 0.3, 0.6, 1 M (cid:12) ), 3 envelope fractions (10, 50, 90% by mass), and 2envelope specific entropies (3 × , 5 × erg g -1 K -1 ).The postgenitor masses and envelope fractions were cho-sen to encompass a large range of possibilities for theexplosion process. If the explosion is violent and ejectsa large amount of mass, for example, one might expecta small postgenitor mass. If significant burning takesplace but does not result in large amounts of ejecta,one might expect a large envelope fraction. The finalparameter, the specific envelope entropy, was chosen sothat the higher entropy corresponds to unbound or veryloosely bound white dwarfs, while the lower entropy al-ways corresponds to bound objects.2.2. Input physics
After white dwarfs are set up with the appropriateinitial conditions, we evolve them forward in time withMESA. Our inlist is provided in Appendix A, but herewe describe important settings. We use Type 2 opaci-ties derived from the Opacity Project, OP (Seaton 1995,2005), and enable both diffusion and radiative levitationof all elements being simulated. Furthermore, we restrictthe network of isotopes and reactions to include onlythose isotopes we simulate and no reactions, because theisotopes we added are not expected to undergo nuclearreactions. This prevents numerical errors from creatingspurious elements, which levitation or diffusion mightthen concentrate–a phenomenon we had previously seenin our models.In addition to diffusion, radiative acceleration, andconvection, we introduce an additional source of mixingwith min D mix=1.0 . This minimum diffusion coefficientof 1 cm s − ameliorates numerical problems, such asunphysically sharp composition gradients and unrealis-tically rapid composition fluctuations in the outer layersof the envelope. In physical terms it may correspond tosources of mixing not accounted for in our model, suchas rotational mixing. min D mix ∼ mz/custom mesa 10000.tar.gz.The algorithm for calculating radiative accelerationsis given in Seaton (2005). To summarize, the radiativeacceleration for element k is: g rad ,k = Fc µµ k κ R γ k (1) γ k = (cid:90) σ mtak ( u ) σ ( u ) F ( u )1 − e − u du (2) σ mtak = σ k ( u )(1 − e − u ) − a k ( u ) (3)Here: • u = hνk B T , • σ k ( u ) is the monochromatic cross section for ele-ment k at the wavelength corresponding to u , • σ ( u ) is the cross section for the mixture, • σ mtak ( u ) is the cross section for momentum transferto the atom, • a k ( u ) is a momentum transfer correction factor, Table 1.
Composition of SN Iax postgenitor envelopesIsotope Mass abundance (%) C 42 O 48 Ne 5.3 Mg 0.4 Si 1.5 S 0.4 Ca 0.03 Fe 3.6 Ni 0.3 • k R is the Rosseland mean opacity, • µ is the mean atomic weight of the mixture, • F ( u ) is the flux (assuming a blackbody), • γ k is an integrated cross section ratio.The cross sections and correction factors are both pro-vided by the Opacity Project for 10,000 wavelengths.For each element and each wavelength, OP providesthese quantities in a non-rectangular grid whose dimen-sions are temperature and electron density ( n e ).However, electron density is not a useful variable be-cause MESA uses temperature and density to define theconditions in a zone, not temperature and electron den-sity. Electron density can be converted to physical den-sity through another grid provided by OP, namely thatof electrons per atom for each element at each T and n e : ρ = n e e avg µ (4)= n e µ (cid:80) k f k e k ( T, n e ) (5)where e avg is the number of free electrons per atom forthe mixture, in which each element k has abundance f k .To calculate the radiative acceleration of a zone, wefirst calculate the physical density of every point on theT- n e grid, given the composition of the zone. We thenfind the 16 grid points closest to the zone’s temperatureand density, and calculate the radiative acceleration forthose 16 points. Cubic interpolation is then used tocalculate the acceleration at the zone’s temperature anddensity: g rad ( T, ρ ) = (cid:88) i =1 4 (cid:88) j =1 (log T ) i (log ρ ) j c ij , (6)where the 16 c ij constant coefficients are derived by fit-ting to the 16 grid points. We use cubic interpolationinstead of linear interpolation in order to preserve the continuity of the derivatives of opacity with respect totemperature and pressure, which MESA requires for itshydrostatic solver. RESULTS3.1.
Overview
Our grid of 24 models resulted in a diversity of out-comes, listed in Table 2. In the majority of cases, thesimulation followed a canonical pattern, explained inmore detail in Section 3.2. The WD would initiallycool and dim. After some time, there is a rapid re-brightening and reheating event. Depending on the peaktemperature, radiative levitation could become impor-tant at this stage, creating an atmosphere dominated bynickel and iron. Afterwards, the WD cools, radiative lev-itation fails, and gravitational settling takes over. TheWD enters onto the cooling track and follows it there-after. 3.2.
Canonical case
Figure 1 illustrates the behavior of a typical postgeni-tor in our simulations. The initial cooling and dimmingis due to the outer layers, which have a short thermaltimescale, radiating away heat. During this phase theouter envelope is convective–the constant entropy en-velope exactly fulfills the Schwartzchild criterion, andpreferential cooling of the outer layers only increases thetemperature gradient. Although there is abundant heatburied deeper in the envelope, this heat has not yet hadtime to diffuse out. When the heat does diffuse out,it results in the reheating and rebrightening event seenin all three diagrams. The WD becomes very hot andbright. At some point it begins to cool again, followingthe normal cooling track for WDs.A few features of Figure 1 are worth noting. First,the final cooling track is to the left of the initial cool-ing track, and also to the left of the rebrightening track.Since L = 4 πR σT by definition, this leftward shift in-dicates a substantial decrease in radius. Second, higher Table 2.
Simulation outcomesMass ( M (cid:12) ) Envelope entropy (erg g -1 K -1 ) Envelope fraction Outcome0.15 3 × × × × × × × × × × × × × × × × × × × × × × × × masses lead to higher temperatures and higher luminos-ity.Figure 2 and Figure 3 show the time evolution of tem-perature and luminosity for all low entropy models. Ascan be seen, all scenarios follow the canonical patternof dimming, rapid re-brightening, and re-dimming. Thetimescales, however, are drastically different. Higher-mass WDs evolve much faster, as do WDs with low en-velope fraction.Figure 4 shows the surface abundances of all elementsover time in one specific model, namely the one withWD mass 0.6 M (cid:12) , 10% envelope fraction and lower en-velope entropy. During the initial cooling stage, surfaceabundances are constant due to the convective zone inthe envelope. In fact, the convective zone ensures thatthe entire WD has near-uniform composition through-out this stage. After the re-brightening event at ∼
70 yr,the envelope becomes radiative, and the the high surfacetemperatures cause radiative acceleration of iron andnickel toward the surface. These two elements are pref-erentially levitated because they have a large number oflines, thus fulfilling g rad > g . As the postgenitor cools,radiative levitation eventually fails to hold heavy ele-ments aloft, and they fall out of the photosphere. Grav- itational settling then takes over, interacting with theartificially injected mixing (min D mix = 1) to createstable surface abundances after 1 Myr.3.3. Abnormal cases
The abnormal cases are simulations that fail to reachthe cooling track. In some cases this is because themodel has positive total energy, and is therefore notgravitationally bound. This occurs in the low massWDs with high entropies and high envelope fractions.In other cases the model is gravitationally bound, buthas enough energy in the envelope that the luminosityexceeds Eddington luminosity, the envelope expands totens or hundreds of solar radii, and MESA stalls. Thisoccurs preferentially in the high mass WDs with thinenvelopes, as high mass WDs have enough gravity tokeep the envelope bound. In both cases significant massloss is expected, though we do not attempt to calculatesuch mass loss here. One particularly interesting exam-ple of the second case is shown in Figure 5, where thered dwarf dims, rebrightens, dims again, and rebrightensagain for the final time while expanding into a red gi-ant. Not surprisingly, all of these abnormal simulationsoccur when the envelope entropy is high. l o g L ( L ) M l o g L ( L ) M l o g L ( L ) M log T eff (K) l o g L ( L ) M yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr10 yr Figure 1.
HR diagrams for different masses at the sameenvelope fraction (50%) and envelope entropy (3 × ergg -1 K -1 ). These pathological cases are interesting in their ownright. There is no reason, for example, why a survivingWD cannot have an ultra-hot envelope, or why the en-velope cannot puff up and drive a wind. Indeed, Foleyet al. (2016) find P-Cygni features on permitted linesat late times (t >
200 d) in a Type Ia SN, implyingan expanding envelope that the authors attribute to asuper-Eddington wind. The implied velocity of 410 km -1 is consistent with the escape velocity of a R (cid:12) post-genitor. Furthermore, the narrow forbidden lines have asimilar velocity as the expanding photosphere, implyingthey are also due to the wind. l o g T e ff ( K ) M l o g T e ff ( K ) M l o g T e ff ( K ) M log Age (yr) l o g T e ff ( K ) M f env = 0.9 f env = 0.5 f env = 0.1 Figure 2.
Effective temperature as a function of age for alllow entropy scenarios.
Unfortunately, it is difficult to simulate these highlyinflated mass-losing objects, because MESA runs intonumerical difficulties when mass becomes unbound orwhen the luminosity becomes super-Eddington nearthe photosphere. Although it is possible to intro-duce a wind, there is no guarantee that existing windprescriptions–developed for RGB and AGB stars–willbe suitable for these peculiar objects. There is also noguarantee that a wind would help with the convergenceproblems. Lau et al. (2012) model a similar instabilityassociated with super-Eddington winds in AGB stars,but their simulation also crashed due to numerical prob-lems. Thus, while these hyper-inflated carbon/oxygenare interesting, our MESA models are not trustworthyrepresentations, so we set aside these pathological casesand focus on the ‘Normal’ outcomes in Table 2. DISCUSSION l o g L ( L ) M l o g L ( L ) M l o g L ( L ) M log Age (yr) l o g L ( L ) M f env = 0.9 f env = 0.5 f env = 0.1 Figure 3.
Luminosity as a function of age for all low entropyscenarios.
In the previous section, we presented the salient char-acteristics of our simulations. It is worth discussingwhich aspects of our simulations are believable andwhich should be taken with a grain of salt. After thisdiscussion, we will compare our simulation results to theobservations of LP 40-365, a candidate Iax postgenitor.4.1.
Decline-rise-decline pattern
The most prominent characteristic of all the mod-els that do not become unbound or swell into red gi-ants is that they have a dimming phase, followed by arapid re-brightening, followed by a protracted coolingphase akin to those of ordinary WDs. This decline-rise-decline pattern is robust across a wide range of postgen-itor masses, envelope fractions, envelope entropies, andcompositions. However, its existence may depend on theconstant entropy assumption, which creates a tempera-ture profile that rises sharply with density ( T ∝ ρ / ), log Age (yr) M a ss a b un d a n c e c12o16 ne20si28 mg24s32 ca40fe56 ni58 Figure 4.
Surface abundances as a function of age. Plotting0.6 M (cid:12) , 10% envelope fraction, low entropy scenario. At veryearly times ( <
100 yr), the white dwarf composition doesnot change because there is a convection zone extending tothe surface. Radiative levitation is significant from 100 to10,000 yr, after which the white dwarf has cooled enough forgravitational settling to take over. which buries heat deep in the envelope. Since the ther-mal timescale increases rapidly with depth, the outerlayers cool before the inner layers can react, as shownin Figure 6. Eventually, heat from the interior diffusesand heats the envelope from the inside out. When theheat reaches the surface, the postgenitor is near peak lu-minosity. After a thermal time near the star’s core, theentire stars cools and descends the WD cooling track.This behavior is expected so long as we believe theconstant entropy assumption. This is not an obviouslybad assumption, as one might expect vigorous mixing inthe aftermath of a supernova, which would flatten outthe entropy gradient. Nevertheless, the large unknownsin the explosion mechanism prevents us from proving theassumption is accurate, so the initial decline and subse-quent rebrightening should be embraced cautiously.To test the sensitivity of our models to the constantentropy assumption, we created a model with a 0.6 M (cid:12) WD, 10% envelope fraction, and the following entropyprofile: d ln( s ) dr = 0 . H , (7)where S is the entropy per unit mass, H is the local scaleheight, and s=2 × erg g -1 K -1 at the bottom of theenvelope. This model has a qualitatively similar evolu-tion to the constant entropy model. There is a similardecline-rise-decline pattern, with a uniform compositionprofile during the decline phase, a photosphere domi-nated by heavy elements at peak luminosity, and grav-itational settling taking over during the cooling phase. T eff (K)101234 l o g L ( L ) yr 10 yr 10 yr 10 yr 10 yr (a) HR evolution log Age (yr) l o g R ( R ) (b) Radius vs age Figure 5.
A particularly interesting atypical result. Thisis the 1 M (cid:12) , 50% envelope fraction, high envelope entropymodel. We conclude that our results are not sensitive to theexact shape of the entropy profile.4.2.
Postgenitor luminosity and evolutionary timescales
The evolution of the postgenitor is largely controlledby the radiative diffusion of heat out of the deep in-terior, which is determined by the opacity structure ofthe star. A unique feature of our postgenitor modelsis their relatively large abundances of iron group ele-ments in their outer layers. The high opacities createdby these elements, coupled with the unusual initial con-ditions (constant envelope entropy, high iron content)of our models, creates the characteristic dimming andbrightening evolution described above.The initial dimming phase is easy to understand.Heat is transported outward by convection on a ther- mal timescale t therm , causing the outer layers of the starto cool. As the cooling front moves inward, t therm atthe base of the cool envelope increases (see red curvein Figure 6), and so the emerging luminosity decreaseswith time. This behavior continues until the coolingfront reaches a point in the star where t therm has a lo-cal maximum. This maximum can be easily seen atlog(1 − q ) ∼ − t therm initially has a local maximum fairlydeep in the star near the iron opacity peak at temper-atures of T ∼ K, which is especially important dueto the high iron abundance of our models.Eventually, heat diffusing into this layer from belowraises its temperature substantially, thereby decreasingits opacity which scales approximately as κ ∝ T − . , al-lowing more heat to diffuse from below. The layer heatsmore, further decreasing its opacity, causing a runawayprocess so that a heating wave runs through the enve-lope toward the surface of the star. The photospherictemperature and luminosity increase suddenly, as shownin Figures 2 and 3. The luminosity remains large forroughly one thermal time at the base of the high en-tropy envelope, after which the star steadily descendsthe WD cooling track.The timescale of the WD rebrightening can be esti-mated via the thermal time in the layers below the ironopacity peak. This timescale is t therm ∼ H K ∼ H κρ c P σ B T , (8)where K is the thermal diffusivity or thermal diffusioncoefficient, H ∼ P/ ( ρg ) is the local scale height, κ is theopacity, c P is specific heat at constant pressure (com-puted by MESA from the EOS), and the other variableshave their usual meaning. For bound-free and bound-bound opacity created by iron group elements, the opac-ity is approximately (Hansen et al. 2004) κ ∼ κ ρT − . , (9)where κ ∼ × Z cm K . g − and Z is the metal-licity. Deep in the star, H ∼ r , where r is the localradial coordinate, ρ ∼ M/ πr , and the temperaturecan be approximated from the virial relation, T ∼ GM µm p k B r (10)where µ is the mean molecular weight. This virial rela-tion holds because the high entropy envelope of our WD log (1-q) l o g T ( K ) log (1-q) l o g R ( c g s ) log (1-q) T h e r m a l t i m e s c a l e ( y r ) log (1-q) l o g S ( K ) Figure 6.
Temperature, opacity, thermal timescale, and entropy as a function of depth for selected ages of the 0.6 M (cid:12) , lowentropy, 10% envelope fraction model. The mass coordinate of the x -axis is log(1 − q ) = log(1 − M r /M ), such that the star’ssurface lies to the right. Colored lines indicate different ages since explosion. From early to late times: at the beginning, theouter layers cool while the interior stays at the same temperature. Around 71 yr, heat from the interior begins diffusing out,but has not yet reached the surface. The next profile shows the situation near peak luminosity. After that, the postgenitorcools. By 56 Myr, which the last profile depicts, the entire envelope has had time to cool, and the postgenitor is well on its waydown the cooling track. Note that the core-envelope transition is at log(1-q) = -1. models is non-degenerate and well-approximated by anideal gas.The relevant brightening time is set by the minimal t therm in layers below the iron opacity bump. Figure 6demonstrates the peak in t therm near the iron opacitybump, which prevents heat from these layers from dif-fusing outward. However, at larger depths the opacity islower and the t therm is shorter, such that heat diffusingfrom deeper in the star warms the gas in the iron opac-ity bump. Combining the above relations, this happens on the thermal timescale: t therm ∼ π k B κ acG µ m p M − T / ∼ × yr (cid:18) µ . (cid:19) − (cid:18) M . M (cid:12) (cid:19) − (cid:18) T K (cid:19) / . (11)Equation 11 provides a crude estimate of the diffusiontimescale that corresponds to the age at which the lu-minosity increases in Figures 1 and 3. The iron opacitypeak is around T Fe ∼ K at the densities present inour mid-envelope, so we should evaluate equation 12 atsomewhat warmer temperatures of ∼ K. The bright-ening age is smaller in higher mass postgenitors, largelybecause the internal temperatures are larger such thatthe iron opacity peak lies closer to the surface where0the density is lower and the diffusion time is smaller.In equation 11, the appropriate mass is the core mass,(1 − f env ) M , such that smaller envelope fractions havefaster evolution timescales.The timescale of the peak luminosity in Figure 3 isgiven by the photon diffusion time near the base of thehigh entropy envelope, t dif ∼ r ρκc (12)Replacing the density, temperature, and opacity as doneabove, this equates to t dif ∼ π k / B κ ( Gµm p ) / c M − / r − / . (13)Because the luminosity of the postgenitor is poweredby gravitational energy release as it contracts into a WD,we must evaluate equation 12 where the gravitationalenergy release is largest, i.e., where P gas ∼ P deg . At thislocation, we find the usual WD scaling relation r ∼ h G ( µ e m p ) / m e M − / . (14)Combining equations 13 and 14, we have t peak ∼ t dif (15) ∼ π k / B ( µ e m p ) / m / e κ G µ / m / p hc M − / (16) ∼ (cid:18) M . M (cid:12) (cid:19) − / . (17)The corresponding peak luminosity is simply the post-genitor’s gravitational binding energy divided by the dif-fusion timescale, L peak ∼ E bind t dif ∼ GM Rt dif ∼ L (cid:12) (cid:18) M . M (cid:12) (cid:19) / . (18)Although crude, these estimates approximately predictthe timescale and luminosity of the peaks in Figures 2and 3, and more importantly, they largely explain thesteep scalings with postgenitor mass, which are due tothe larger binding energies and lower opacities in moremassive postgenitors.4.3. Radiative levitation
During the hot and bright phase of our models, radia-tive levitation becomes strong enough to drive iron and nickel towards the surface, making them the most abun-dant elements at the surface. These elements are pref-erentially levitated because they have the most abun-dant absorption lines, and thus the highest momentumtransfer cross sections. Strong radiative levitation ofother rare elements (e.g., strontium or tellurium) is alsoprobable, but we do not include these elements in ourgrid because the Opacity Project does not provide opac-ities for them. The transition between a heavy elementphotosphere and a light element photosphere is abrupt,as can be seen in Figure 4. For a 0.6 M (cid:12) WD, thetransition occurs around T eff = 100,000 K. This criticaltemperature drops to 50,000 K for a 0.3 M (cid:12) WD, andclimbs to 250,000 K for a 1 M (cid:12) WD. (The 0.15 M (cid:12) WDs in our grid do not become hot enough for levita-tion.) The transition temperature has a simple physicalexplanation: it is the point at which g rad = g at thephotosphere for a given element.Figure 7 shows radiative acceleration as a function ofposition in the star for the hottest model in Figure 4. Itcan be seen that at this point in time, iron and nickelhave g rad > g while the other elements do not. Notsurprisingly, these two are by far the most abundantelements in the photosphere. As time passes, iron andnickel increase in abundance until g rad = g , which occurs ∼
20 yr after the profile shown. Around peak brightness(2600 yr), g rad for the other elements also approach g (Figure 8), and their surface abundances reach an equi-librium (Figure 9). Then the WD cools, and all radia-tive accelerations drop below gravitational accelerationat ∼ log (1-q) l o g g r a d ( c g s ) log_go16 ne20mg24 si28s32 ca40fe56 ni58 Figure 7.
The onset of radiative levitation. This plot shows g rad as a function of depth for all elements in the 0.6 M (cid:12) ,10% envelope fraction, low entropy model. This snapshotwas taken at 218 yr, when T eff = 150,000 K. log (1-q) l o g g r a d ( c g s ) log_go16 ne20mg24 si28s32 ca40fe56 ni58 Figure 8.
Same as Figure 7, but at peak brightness. Noticethat radiative acceleration for all elements is close to g atthe photosphere. This snapshot was taken at 2600 yr, when T eff = 200,000 K. possible that clumps of metal over-densities will form,with most of the flux escaping through gaps betweenthe clumps, and the clumps themselves shielding mostatoms inside from radiative levitation.An additional complication is that when the luminos-ity approaches or exceeds the Eddington limit–whichoccurs for most of our higher mass models–the atmo-sphere develops an inhomogenous porous structure andthe effective opacity is greatly reduced. This effect hasbeen suggested for many super Eddington scenarios, in-cluding nova outbursts (Kato & Hachisu 2005) and su-permassive stars (Shaviv 2000). Three-dimensional hy-drodynamic simulations of radiation-dominated massivestar envelopes (Jiang et al. 2015, 2017) reveal a com- A b un d a n c e c12o16 ne20si28 mg24s32 ca40fe56 ni58 Figure 9.
Abundance profile at 2640 yr, corresponding tothe same model and timestep as Figure 8. Note the highsurface abundance of iron and nickel due to strong radiativelevitation. A b un d a n c e c12o16 ne20si28 mg24s32 ca40fe56 ni58 Figure 10.
Abundance profile for the model from Figure 9,at very late times (56 Myr) when the WD has cooled. plex set of phenomena in super-Eddington atmospheresincluding shocks, porous atmospheres, and oscillations.MESA uses 1D models and cannot accurately model thisporosity, which, by reducing the effective opacity, mayalso reduce the radiative acceleration. Nor do our MESAmodels take winds into account. It is known that massloss strongly hampers the effects of diffusion and radia-tive levitation (Unglaub & Bues 1998; Matrozis & Stan-cliffe 2016) by removing levitated elements and pushingthe convective zone deeper into the star.The observational evidence for radiative accelerationindicates that extreme over-abundances of heavy ele-ments are possible, but a photosphere dominated byheavy elements is not. Werner et al. (2017) took UVspectra of two extremely hot DO WDs ( T eff = 115 , K and T eff = 125 ,
000 K) with moderate surface gravity2(log g = 7 ± . T eff = 65 , − ,
000 K) andfound light metals with subsolar abundances and iron-group elements with 1-100x solar abundances, whichthey interpret as the result of gravitational settling andradiative levitation. Hoyer et al. (2018) searched fortrans-iron elements in hot DO WDs and found very highabundances, indicating that radiative levitation is act-ing. The most extreme example, PG 0109+111, hasa tellurium abundance six orders of magnitude greaterthan solar. At a mass abundance of 6.2 × − , it isthe most abundant metal in the photosphere. The factthat hot WDs have been detected with extreme trans-Fe over-abundances but sub-solar intermediate mass el-ement abundances indicates that radiative levitation isnot completely overpowered by thermohaline diffusion.Radiative levitation is certain to be an important ef-fect in the luminous phase of our WDs, and we expectover-abundances of heavy elements. However, heavy el-ements are unlikely to become the dominant componentof the atmosphere.4.4. Candidate Iax postgenitor stars
Several peculiar WDs have been recently discoveredthat could be Iax postgenitors. Kepler et al. (2016) dis-covered a WD with an oxygen-dominated photosphere(SDSS J124043.01+671034.68) with no trace of carbon.However, since carbon burning is required to producea Iax deflagration, but carbon burning is incompletein such failed explosions, Iax postgenitors are likely tohave substantial carbon abundances. Radiative levita-tion and gravitational settling are unlikely to eliminatecarbon from the photosphere, so we find it unlikely thatJ1240+6710 is a Iax postgenitor. We speculate thatit could be the remnant of an oxygen deflagation aris-ing from an accreting ONe WD that nears the Chan-drasekhar mass. This scenario is similar to the CO de-flagation model we have considered, but beginning withan ONe WD, and would naturally explain the lack of car-bon. Other possibilities include a deflagation in a hybridCONe WD (Bravo et al. 2016) or a CO-ONe WD merger(Kashyap et al. 2018), though it seems likely such eventswould leave some carbon in the bound remnant.Shen et al. (2018) used
Gaia data to discover three hy-pervelocity WD stars. These stars are broadly similarto LP40-365 (see below) in temperature and luminosity,though very different in composition–LP 40-365 is rich inoxygen/neon with little or no carbon, whereas the three Shen objects have carbon in their atmospheres. TheShen objects are possible Iax postgenitors, although theauthors suggest they are instead the companions to Iaprogenitors. In fact, one of them appears to have origi-nated within a supernova remnant, lending credence toan explosive origin. We note here that Ia companionsand Iax postgenitors may look very much alike–they areboth expected to begin as hypervelocity objects withhigh entropy envelopes, inflated radii, and large abun-dances of iron group elements. Even though the goal ofour paper is to model SNe Iax, our models may turn outto be applicable to Ia companions as well.LP 40-365 is a peculiar hypervelocity WD (galacto-centric velocity = 852 km s -1 ) with peculiar abundances,originally discovered by Vennes et al. (2017). The mostabundant photospheric elements are oxygen and neon,followed by intermediate mass elements, while iron andnickel are detected at a number fraction of ∼ − . Theauthors propose that LP 40-365 is the postgenitor of anexploding carbon-oxygen-neon core. Using Gaia data,Raddi et al. (2018) measured the properties summarizedin Table 3. They confirmed the hypervelocity nature ofthe object, measured an abnormally large radius of 0.18 R (cid:12) , and found that it crossed the Galactic disk 5 . ± . Table 3.
LP 40-365 properties, from Raddi et al. (2018)Property Value T eff ± g . ± . R ( R (cid:12) ) 0 . ± . M ( M (cid:12) ) 0 . +0 . − . L ( L (cid:12) ) 0 . ± . ±
10 km s -1 Age > Curious about whether we could explain some of theobserved properties with our models, we looked throughour grid to find the model that most closely matchesLP 40-365. This turned out to be the 0.15 M (cid:12) modelwith 50% envelope fraction and 3 × erg g -1 K -1 enve-lope specific entropy, shown in Figure 11. After a longdimming phase lasting millions of years, this object ex-periences a broad peak in temperature and luminositythat places it close to LP 40-365 on the HR diagram.3 l o g T e ff ( K ) l o g L ( L ) log Age (yr) l o g R ( R ) Figure 11.
Evolution of the 0.15 M (cid:12) model with 50% enve-lope fraction and 3 × erg g -1 K -1 envelope specific entropy.The red crosses indicate the observed properties of LP 40-365at an assumed age of 23 Myr. The actual age is unknown,but probably between 5 and 100 Myr (indicated by the errorbars). The error bars on the measurements are plotted, butare too small to be seen In fact, our model comes strikingly close to match-ing the observed properties in Table 3, as can be seenfrom the color-magnitude diagram in Figure 12. At anage of 23 Myr, our model has a luminosity of 0.21 L (cid:12) ,a radius of 0.19 R (cid:12) , and a temperature of 8977 K. 23Myr ago, LP 40-365 would have been ∼ ∼
25 kpc from the Galactic center, rea-sonable for a SN Iax. We note that 23 Myr is close tothe peak in luminosity and temperature for this model,where evolution is slow, and the WD can linger in thisregion of the HR diagram for tens of millions of years.If LP 40-365 actually originated from the disk and is5 Myr old, a model with slightly higher mass (e.g. 0.2 M (cid:12) ) would be needed to match its properties, as a moremassive WD experiences faster evolution. Despite thisclose match, we do not claim that our model exactlyexplains LP 40-365. A key feature of LP 40-365 thatdisfavors the Iax postgenitor scenario is the low abun-dance of C in the photosphere, whereas in our modelsat this temperature, C is the most abundant element.This could indicate that LP 40-365 is better explainedas the donor star in a binary with an accreting WD G BP G RP G + l o g (( " )) + Figure 12.
Color-magnitude diagram with Gaia BP-RPcolors and G magnitudes. Green: Iax postgenitor candi-date LP 40-365. Red: the three candidates from Shen et al.(2018). Orange: our 0.15 M (cid:12) model with 50% envelopefraction and 3 × erg g -1 K -1 envelope specific entropy, con-verted to Gaia quantities by assuming a blackbody spectrum.Blue: Gaia objects within 100 pc, for reference. Note thatour model stays slightly above the main white dwarf cool-ing track at old age, due to its exceptionally low mass andconsequently high radius. (Shen et al. 2018), or perhaps a partially burnt O-NeWD (Jones et al. 2016). Nevertheless, we have shownthat it is natural for Iax postgenitor models to matchthe temperature, luminosity, and radius of LP 40-365 ata reasonable age.Aside from our CO WD simulations, we also per-formed an exploratory simulation with an oxygen/neoncomposition, matching that of LP 40-365. It encoun-tered numerical problems, but broadly matched thebehavior of the C/O white dwarfs in terms of theirdecline-rise-decline pattern and extreme element over-abundances caused by radiative levitation.4.4.1. Number of Detectable Postgenitors
We can estimate the number of detectable postgeni-tors assuming each type Iax SN produces a high velocityremnant star. The SNe Iax occurrence rate is ∼ ∼
300 years, a SN Iax wouldoccur every ∼ ∼ ∼ ∼ r = 2 kpc de-tection volume is N ≈ πr hn ≈ × , where h=3504pc is the galactic disk scale height, and n=0.14 pc -3 isthe local stellar density. Assuming the Milky Way has2 . × stars, roughly 0 .
2% are within the detectionvolume. Then we expect ∼ × yr as a SNremnant lifetime, we expect only 1 in ∼
70 of detectablepostgenitors could be traced back to a SN remnant, suchthat we expect to see only ∼ .
06 postgenitors associatedwith a SN remnant, potentially in tension with the oneobject traced back to a remnant by Shen et al. (2018).4.5.
Thermohaline mixing
One important effect we have not yet consideredis thermohaline mixing. Thermohaline mixing occurswhen a radiative region (defined by the Ledoux crite-rion) exhibits an inverse composition gradient, i.e., ithas layers of high molecular weight on top of layers oflow molecular weight. If a blob of high molecular weightmaterial is displaced downwards and no heat exchangeoccurs, the blob would be less dense than its surround-ings and float back up. However, if substantial heat ex-change occurs, the blob cools and becomes denser thanits surroundings, thereby continuing to sink. In Earth’soceans, thermohaline mixing gives rise to “salt fingers”–so called because sinking blobs create very salty tendrils,sticking deep into less salty subsurface layers.In stars and WDs, thermohaline mixing has the ef-fect of introducing mixing into radiative regions wheremixing would otherwise be negligible. This mixing is im-portant in scenarios like planetesimal accretion (Bauer& Bildsten 2018) and carbon-enhanced metal-poor stars(Stancliffe & Glebbeek 2008), where heavier elementsaccrete on top of a lightweight atmosphere. In our sce-nario, radiative levitation tends to push heavy elementsupwards, while gravitational settling and thermohalinemixing counteract levitation.To explore the effect of thermohaline mixing on ourWDs, we ran a simulation with thermohaline diffusionenabled. As discussed in Traxler et al. (2011); Zem-skova et al. (2014), thermohaline mixing occurs in re-gions with 1 < R < K/K µ , with R = ( ∇ ad − ∇ ) / ∇ µ ,the thermal/composition gradients have their usual def-inition, and K µ is the composition diffusivity. In ther-mohaline unstable regions, we set the thermohaline dif-fusion coefficient D therm (cid:39) C ( K − R K µ ) / ( R −
1) ac-cording to Equation 4 of Vauclair & Th´eado (2012) (it- self derived from Denissenkov 2010) with their recom-mended coefficient of C = 120. Enabling thermoha-line mixing dramatically changes the chemical compo-sition evolution. The composition remains nearly con-stant despite radiative levitation, maintaining the sameuniform abundances that we start the simulation with.This occurs because any increase in the abundances ofheavy elements due to levitation increases the molecu-lar weight gradient, which enhances thermohaline mix-ing and mixes the heavy elements back down. This alsoexplains why calcium, the least abundant element, cansteadily increase in abundance while the others cannot–because calcium contributes negligibly to the molecularweight.Our thermohaline models run into numerical problemssome time after they pass peak luminosity, leading tothe diffusion solver failing unless the timescales are verysmall. Due to these problems, we leave the full explo-ration of thermohaline mixing under these conditions tofuture work, but the implication of these findings is thatthe surface abundance enhancement of iron group ele-ments will be strongly reduced by thermohaline mixingrelative to the predictions of Section 4.3.We can also analytically estimate the equilibriumcomposition gradient by equating a radiative levitationtimescale t rad = ( v rad d ln µ/dr ) − to a thermohalinemixing timescale t therm = ( D therm d ln µ/dr ) − . Thelength scale (cid:96) on which we expect the composition tovary is then (cid:96) ∼ (cid:18) CKHv rad ( ∇ − ∇ ad ) (cid:19) / . (19)In our models, equation (19) predicts (cid:96) (cid:28) H , so weexpect radiative levitation to produce very weak com-position gradients when competing with thermohalinemixing, in accordance with the results of our MESAmodels. However, we note that rotation and magneticfields, which are not included, may limit the effects ofthermohaline diffusion. We leave a realistic assessmentof these effects to future work. CONCLUSIONWe have modeled type Iax supernova postgenitor starswith MESA with a range of initial conditions, account-ing for uncertainties in their masses and post-explosionstructure. Not surprisingly, we obtained a wide range ofbehaviors. Most of our models followed a canonical be-havior, starting as hot WDs with abnormally high radiithat initially cool and dim. Later, as heat leaks outof the deeper interior, the envelope opacity is reduced,allowing faster radiative diffusion. The stars then be-come much hotter and brighter on timescales of years5to millions of years after the supernova, depending onthe star’s core and envelope mass. At peak brightness,all but the lightest WD models have over-abundances ofiron group elements in their photospheres due to radia-tive levitation. Afterwards, the WDs shrink in radius,cooling and dimming similar to normal WD cooling se-quences. Our highest entropy models became unbound,super-Eddington, or inflate into red giants, indicatingthat some Iax postgenitors could appear as luminouscool stars rather then hot blue stars.Although the prospect for observing these postgeni-tors in the early aftermath of a SN Iax is remote, it isnot unlikely that a known WD inside the Milky Wayis such a postgenitor. In fact, we already have fourcandidates, including LP 40-365, which our lowest massmodels naturally mimic in luminosity and temperatureat a plausible age. Future models for such stars canbe improved with a better implementation of thermo-haline mixing and mass loss, and realistic estimates forthe post-explosion structure. As these models improve,we encourage further deeper observational searches forpeculiar WDs and subdwarf remnant stars of variousflavors of thermonuclear supernovae. ACKNOWLEDGMENTSWe thank Ken Shen, Evan Bauer, and Lars Bildstenfor useful discussions, and we thank the referee for avery constructive and thorough report. This work wassupported by the Heising-Simons Foundation throughGrant
Software: