A class of partly burnt runaway stellar remnants from peculiar thermonuclear supernovae

We report the discovery of three stars that, along with the prototype LP 40-365, form a distinct class of chemically peculiar runaway stars that are the survivors of thermonuclear explosions. Spectroscopy of the four confirmed LP 40-365 stars finds ONe-dominated atmospheres enriched with remarkably similar amounts of nuclear ashes of partial O- and Si-burning. Kinematic evidence is consistent with ejection from a binary supernova progenitor; at least two stars have rest-frame velocities indicating they are unbound to the Galaxy. With masses and radii ranging between 0.20-0.28 Msun and 0.16-0.60 Rsun, respectively, we speculate these inflated white dwarfs are the partly burnt remnants of either peculiar Type Iax or electron-capture supernovae. Adopting supernova rates from the literature, we estimate that ~20 LP 40-365 stars brighter than 19 mag should be detectable within 2 kpc from the Sun at the end of the Gaia mission. We suggest that as they cool, these stars will evolve in their spectroscopic appearance, and eventually become peculiar O-rich white dwarfs. Finally, we stress that the discovery of new LP 40-365 stars will be useful to further constrain their evolution, supplying key boundary conditions to the modelling of explosion mechanisms, supernova rates, and nucleosynthetic yields of peculiar thermonuclear explosions.


INTRODUCTION
LP 40−365 (GD 492;Luyten 1970;Giclas et al. 1970) is a high-velocity star unbound to the Galaxy with a unique composition: it has an O/Ne-dominated atmosphere sprinkled with the ashes of incomplete O-and Si-burning (Vennes et al. 2017). LP 40−365 has been proposed as a partially burnt runaway white dwarf that survived disruption by a thermonuclear supernova, specifically a peculiar group ⋆ E-mail: roberto.raddi@fau.de that is suggested to experience pure deflagrations, leading to subluminous explosions that do not completely disrupt the white dwarf accretors (SN 2002cx-like, or SNe Iax;Li et al. 2003;Phillips et al. 2007;Jordan et al. 2012;Foley et al. 2013;Kromer et al. 2013Kromer et al. , 2015Fink et al. 2014;Jha 2017). In Raddi et al. (2018b), hereafter Paper I, we presented supportive evidence of the SNe Ix survivor interpretation by measuring a super-solar Mn abundance in LP 40−365 that is compatible with theoretical nucleosynthesis yields of near-Chandrasekhar-mass (M Ch ) thermonuclear explosions (Seitenzahl et al. 2013;Seitenzahl & Townsley 2017).
In Raddi et al. (2018a), hereafter Paper II, we used the precise parallax available in the Second Gaia Data Release (Gaia DR2; Gaia Collaboration et al. 2018) to estimate the radius of LP 40−365 to be R ≃ 0.18 R ⊙ , i.e. one order of magnitude larger than canonical white dwarf radii (Tremblay et al. 2017). We interpreted this finding with the star being currently inflated as a consequence of the supernova explosion (Jordan et al. 2012;Kromer et al. 2013;Shen & Schwab 2017). The scenario of a post-explosion expansion is compatible with the low rotational velocity of the star (v sin i ≤ 50 km s −1 ). Studying the kinematics of LP 40−365, we confirmed it is gravitationally unbound from the Milky Way, having a rest frame velocity of v rf ≃ 850 km s −1 that is about 1.5 times larger than the escape velocity at the corresponding Galactocentric radius. Considering that LP 40−365 has benefitted from the Galactic rotation, we estimated it gained an ejection velocity of ≃ 600 km s −1 at the moment of the supernova explosion. In the single-degenerate scenario, considering negligible contribution from asymmetric mass loss, this velocity could be achievable via the ejection from a compact binary (1 hr orbital period), e.g. with a massive He-burning companion (Wang & Han 2009;Wang et al. 2013).
Among many exciting discoveries, Gaia DR2 has led to another possible breakthrough in the field of thermonuclear supernova research with the identification of three new hyper-runaways with v rf 1000 km s −1 (Shen et al. 2018b). These stars form a new, relatively homogeneous spectral class of subluminous, CO-rich objects likely connected to the explosion mechanism of supernovae Type Ia (SN Ia) that is known as dynamically driven double-degenerate doubledetonation scenario (D 6 ; Shen et al. 2018a).
LP 40−365 is interpreted as the formerly accreting white dwarf in a single-degenerate mass-transferring binary, also referred to as the bound remnant. On the other hand, the D 6 stars are interpreted as the former compact donors in close, double-degenerate binaries. Also the D 6 stars are suggested to have reached their rest-frame velocities by conserving their high orbital velocities at the moment of explosion. Like LP 40−365, the D 6 stars are currently thermally bloated; however, their expansion might be caused by the pre-explosion tidal forces and the interaction with the supernova ejecta (Shen & Schwab 2017). Prior to Gaia DR2, a hyper-runaway He-rich subdwarf star, US 708 (Hirsch et al. 2005;Heber 2016), was proposed as the fastest known star that may have been the former subdwarf donor in a thermonuclear supernova event (Geier et al. 2015). Following the identification of the D 6 stars, Shen et al. (2018b) have suggested a connection with US 708, which they proposed as representative for a possible later stage on the way back towards the canonical white dwarf cooling sequence. All together, LP 40−365, the D 6 stars, and US 708 likely represent the first direct evidence of a wider class of binary progenitors for thermonuclear supernovae (Wang & Han 2012;Maoz et al. 2014).
We have performed spectroscopic follow-up of candidate high-velocity stars, which we selected by mining Gaia DR2. Our search has successfully led to the identification of two additional stars that, along with LP 40−365, we propose to form a distinct class of partly burnt, post-supernova runaway stars (hereafter, "LP 40−365 stars"). Furthermore, by mining the Sloan Digital Sky Survey (SDSS) spectro- −3757 are shown as circle, square, and diamond symbols, respectively (correspondence between colours and symbols will be maintained in the following figures). The three D 6 objects, identified by Shen et al. (2018b), are marked by star-shaped symbols. US 708, for which we adopt a spectroscopic distance of 8.5±1.0 kpc (Geier et al. 2015), is plotted as a pentagon symbol. Error bars are smaller or of the order of the symbol sizes. The intrinsic colours of normal stars are interpolated from the Pickles (1998) spectral library (luminosity class V, III, I). Cooling tracks for 0.4, 0.6, 1.0 M ⊙ white dwarfs (from top to bottom) are derived by scaling a grid of synthetic spectra in the T eff = 6000-90,000 K range (Koester 2010) with the Montreal mass-radius relation for pure-hydrogen atmospheres (Wood 1995;Bergeron et al. 2001;Fontaine et al. 2001) and convolving them with the Gaia DR2 filter profiles.
scopic database (Smee et al. 2013;Abolfathi et al. 2018), we have identified two related objects: one is very likely also an LP 40−365 star, while the other requires further follow-up observations. In this manuscript, we present follow-up observations of the Gaia candidates, and the identification of the SDSS candidates. In Section 3, we detail the spectral analysis of the LP 40−365 stars identified with Gaia, and we reanalyse LP 40−365 itself including new Hubble Space Telescope (HST) near-ultraviolet (near-UV) spectroscopy. We present the Galactic orbits of the new LP 40−365 stars in Section 4. In Section 5, we discuss our results in a broader context, focusing on: i) the chemical composition, compared to theoretical simulations of white dwarfs that are predicted to survive peculiar thermonuclear explosions; ii) physical parameters (mass and radius); iii) evolutionary status, drawing a comparison with theoretical predictions; and iv) kinematics, birth places and binary progenitors. Finally, in Section 5.2, we model a population of LP 40−365 stars to estimate the end-of-mission detection limits of Gaia within 2 kpc from the Sun, deriving the distribution of astrometric parameters. In our conclusions, we summarise the remarkable similarities among the observed stars. Table 1. Astrometric and photometric parameters of J1603−6613 and J1825−3757, including broad-band photometry from GALEX (NUV; Morrissey et al. 2007), SkyMapper (uvgriz;Wolf et al. 2018), and Vista Hemisphere Survey (J; McMahon et al. 2013 We searched for high-velocity-star candidates by selecting Gaia objects having relatively large transverse velocities (v t 300 km s −1 ) and falling within the extended colour space below the main sequence in the Gaia colour-magnitude Hertzsprung-Russell (HR) diagram ( Fig. 1) that is occupied both by LP 40−365 and the D 6 stars, i.e. delimited by Gaia colours of B p − R p 0.7 and absolute magnitudes of M G 2.5. Limiting the parallax precision to ̟/σ ̟ > 5, our criteria included LP 40−365 and D6-2, but excluded the two other D 6 stars. The spectroscopic follow-up has led to the identification of several canonical hot subdwarfs and white dwarfs, which we will present elsewhere. Two stars stood out -Gaia DR2 5822236741381879040 and Gaia DR2 6727110900983876096 (hereafter shortened to J1603−6613 and J1825−3757, respectively, based on their equatorial coordinates) -as possessing peculiar spectra, clearly resembling that of LP 40−365 (see Fig. 2). J1603−6613 was previously unknown, while we previously obtained a low-resolution of J1825−3757 (Fig. 2) on the base of which we had classified it as a hot subdwarf (Raddi et al. 2017). Follow-up observations are detailed in the next sections, and the journal of observations is listed in Table 2. The Gaia DR2 astrometry and photometry of the new stars are listed in Table 1. J1603−6613 has an apparent magnitude G = 17.84 that, combined with its parallax ̟ = 0.57 ± 0.10 mas, corresponds to an absolute magnitude M G = 6.62, which is similar to that of LP 40−365. This new star is slightly bluer than LP 40−365, suggesting it may be hotter. The Gaia colour excess factor ( f BP + f R P )/ f G = 1.6 might indicate some problems with the photometry, given that the star is observed in a relatively crowded field. We note an interstellar extinction in the V-band A V = 0.28-0.32 (Schlegel et al. 1998;Schlafly & Finkbeiner 2011), i.e. about four times larger than the that along the LP 40−365 sightline. J1825−3757 is located closer to the Sun, having a parallax of ̟ = 1.08 ± 0.06 mas. It is intrinsically bluer and about 3 mag brighter than LP 40−365 and J1603−6613, implying it is both hotter and larger than the other two stars. The total interstellar extinction is comparable to that of J1603−6613.

Low-resolution spectroscopy
We first observed J1603−6613 on 2018 May 6 with the Goodman Spectrograph (Clemens et al. 2004) mounted on the 4.1 m Southern Astrophysical Research (SOAR) telescope at Cerro Pachón in Chile. We obtained 8 × 300 s consecutive spectra using a 930 line/mm grating with wavelength coverage between 3600−5200Å. We used a 1. ′′ 01 slit and bracketed our exposures with a 60 s Fe arc lamp, since no skylines exist in this spectral region for wavelength calibration. With a dispersion of 0.84Å/pixel our resolution is roughly 2.8Å.
The spectra were optimally extracted (Horne 1986) using standard starlink routines (Currie et al. 2014), after applying bias and quartz-lamp flat corrections using the software pamela. We used the software package molly (Marsh 1989) to wavelength calibrate the spectra, apply a heliocentric correction, and perform a final weighted average of the one-dimensional (1D) spectra, and used the spectrophotometric standard EG 274 for flux calibration. We did not detect any radial velocity shift between the exposures, within the precision of the wavelength calibration, and thus coadded all spectra to obtain a spectrum with a signal-to-noise ratio of S/N ≃ 20 per pixel, shown in Fig.2. Noting the striking similarity with LP 40−365, we used the strongest lines to estimate the radial velocity with respect to the Sun via cross correlation with our synthetic model of LP 40−365, obtaining v rad = −480 ± 10 km s −1 that confirmed J1603−6613 as a high-velocity star.  (Raddi et al. 2017) and the SOAR/Goodman spectrum o J1603−6613 have a resolution of ≃ 8 and 2.8Å, respectively. The WHT/ISIS spectrum of LP 40−365 from Paper I has a 2Å resolution. All the spectra are plotted in the heliocentric frame, so spectral lines are shifted by the radial velocity of each star.

Intermediate-resolution spectroscopy
We followed up J1603−6613 again on 2018 July 9, with the multi-wavelength intermediate-resolution spectrograph X-shooter (Vernet et al. 2011), mounted at the Cassegrain focus of the ESO Very Large Telescope (VLT) UT2 in Cerro Paranál (Chile). We binned the detectors (1×2 in the spatial and dispersion directions, respectively) and we used narrow slits of 1. ′′ 0, 0. ′′ 9, and 0. ′′ 9, in the UVB, VIS, and NIR arms, respectively, taking six exposures in each arm and achieving total exposure times of 7320, 7500, and 7800 s, in order to reduce the impact of cosmic rays and to confirm the absence of radial-velocity variability. The mean seeing was of 0. ′′ 6. The data were reduced using the standard ESO pipeline Reflex (Freudling et al. 2013), and the molecfit software (Smette et al. 2015;Kausch et al. 2015) was used to remove the telluric absorption of Earth's atmosphere. Again, we did not observe any significant radial-velocity variation over short timescales, so we coadded the individual exposures to obtain an average spectrum, confirming the same blue-shifted radial velocity. The coadded spectrum has a S/N ∼ 100 in the blue/visual spectral range, with resolving power R = 5400, 8900, and 5600, in the UVB, VIS, and NIR arms, respectively.

Time-series photometry
We performed time-series photometry of J1603−6613 on the night of 2018 May 10 using the imaging mode of the Goodman spectrograph on SOAR. We used an exposure time of 17 s and had a readout time of 2.0 s, and observed for 2.7 hr through a S8612 broad-bandpass, red-cutoff filter (3300-6200Å). Seeing across the duration of the observation was fairly stable, averaging 1. ′′ -1. ′′ 5.
Images were debiased and flat-fielded using an internal dome lamp. We extracted a light curve using standard circular aperture photometry, with a 4.5-px aperture surrounded by an 8-px annulus. The target light curve was normalised using a bright, nearby comparison star with no close companions, and we subtracted out a second-order polynomial to account for airmass changes. A Lomb-Scargle periodogram showed no significant peaks above 0.2 per cent amplitude for periods between roughly 40-4000 s. This lack of variability suggests that J1603−6613, like LP 40−365, is consistent with not-having an unseen close companion.

Intermediate-resolution spectroscopy
In Raddi et al. (2017) we published a low-resolution spectrum (8Å) of J1825−3757 taken at the New Technology Telescope with the ESO Faint Object Spectrograph and Camera v.2 (NTT/EFOSC2; Buzzoni et al. 1984) at the La Silla observatory (shown in Fig.2). Although we classified this star as a likely hot subdwarf, its spectrum shows numerous unusually strong Mg i-ii lines and other ionised metals. The absence of H and He lines resembles the spectral appearance of LP 40−365 and J1603−6613.
We followed up J1825−3757 on 2016 September 27 with VLT/X-shooter, taking single exposures of 521, 550, and 600 s in the blue, visual, and near-infrared arms, respectively. The instrument setup and data reduction steps were the same as those adopted for J1603−6613. The mean seeing during the observations of J1825−3757 was 0. ′′ 7. The achieved resolving power were 4300, 7400, and 5400, in the UVB, VIS, and NIR, respectively. The extracted spectrum has a S/N > 150 in the UVB/VIS arms and > 100 in the NIR arm.

Time-series photometry
We also observed J1825−3757 through a Bessel-V filter on the night of 2018 July 08 with the imaging mode of the Goodman spectrograph on SOAR. We collected 1140 × 5 s exposures, and each exposure had roughly 2.0 s of dead time from readout. Seeing was steady from 1. ′′ 6-1. ′′ 7. As with J1603−6613 we bias-and flat-field corrected the photometry and extracted the light curve using a fixed 5.5 px circular aperture, and subtracted a second-order polynomial to account for airmass changes throughout the 2.2 hr run. We did not see any coherent variability at any period from 15-3500 s, to a limit of at least 0.1 per cent amplitude, consistent with it currently being an isolated star.

HST spectroscopy
We observed LP 40−365 with HST using the Space Telescope Imaging Spectrograph (STIS; Woodgate et al. 1998) on 2018 May 21 and 25, during Cycle 25, for three and two orbits, respectively. The total exposure time amounted to 15 253 s. We used the G230L grating and the near-ultraviolet (NUV) Multianode Microchannel Array (MAMA) detector with a 52 ′′ × 2 ′′ aperture, which delivers a resolving power R ≃ 1000 across the 1570-3180Å wavelength range. We applied the wavelength and flux calibration via standard iraf calstis tasks (Katsanis & McGrath 1998). We measured an average S/N≃ 35 from the reduced spectrum.

SDSS spectra
We also considered the possibility that similar stars may reside within the 4.9 million spectra of the SDSS database (Abolfathi et al. 2018), which could have been previously misclassified due to their unusual spectral features.
We adapted the template method of Hollands et al. (2017), which essentially fits all SDSS spectra against a library of spectral templates. We created a grid of synthetic LP 40−365 spectra, i.e. with T eff ranging from 5000 to 7000 K in steps of 500 K, and then up to 20 000 K in steps of 1000 K. In all cases the log g was fixed at 5 dex, and the abundances fixed to the average of the already identified stars (as justified by the results presented in Section 3). For each model, additional spectral templates were offset by fixed radial-velocity shifts ranging from −1500 to +1500 km s −1 , in steps of 50 km s −1 , resulting in a total of 1098 templates. Finally, we convolved each template with a Gaussian to match the SDSS spectral resolution of 2.7Å. To filter other easily identifiable spectral types, we complemented the grid of LP 40−365 templates with 1843 white dwarf templates (spectral classes DA, DB, DZ), and 1045 stellar templates of spectral types commonly observed by SDSS (see description in Hollands et al. 2017). For each SDSS spectrum, we interpolated all templates onto the SDSS wavelength scale, flux-scaled them to the data, and then evaluated their χ 2 , adopting the template with the lowest χ 2 as the best match. The scaling factor was the only freeparameter, whose optimal value was determined analytically (see eq. 1 in Hollands et al. 2017). The closest-matching template, the corresponding reduced χ 2 ( χ 2 ν ), and the mean  ν vs. spectral S/N for SDSS spectra that match best to an LP 40−365 template. The red dashed curve indicates a cut, below which the corresponding spectra were visually inspected. Two new candidates are circled in purple, whereas J1240+6710 is circled in orange.
spectral signal-to-noise (S/N) between 4500-6000Å were recorded for each SDSS spectrum.
Out of the 3.1 million SDSS spectra with a mean S/N > 2, 1544 were best matched by a LP 40−365 model. Even so, a closest match did not guarantee a good fit. In Figure 3, we show the χ 2 ν of the best-matching templates vs. the SDSS spectrum mean S/N, with the template-T eff indicated by the point-colours. Towards high S/N, the distribution in χ 2 ν becomes bimodal. Points in the upper-branch correspond to spectra where all templates provided poor matches, but with one of the LP 40−365 templates as the least bad. The lower-branch corresponds to more convincing matches to the templates. Notably, most of the points in the lower branch correspond to high-T eff templates, as their metal-lines become weaker, and thus the number of falsepositives increases. The red-dashed line indicates our cut in the χ 2 ν -S/N plane. The spectra corresponding to the 349 objects below the cut were visually inspected and compared with their best fitting templates.

J0905+2510
With a mean S/N of 5.3, a cursory inspection of the SDSS spectrum of J0905+2510 simply appears as a noisy blackbody. However, degrading the resolution to 8Å (Figure 4 6400 6450 Ne i? Figure 4. Normalised SDSS spectra of J0905+2510 (bottom) and J1637+3631 (top, offset by +0.75). Their best matching "LP 40−365 star" templates are shown in red. For clarity, the spectra and templates have been convolved by Gaussians of FWHM = 8Å to better show the strongest transitions (labelled). For J1637+3631, an inset shows a tentative Ne i detection (at the nominal SDSS resolution). Both spectra are shown in the heliocentric frame. Table 3. Astrometry, photometry, and spectroscopic parameters for the objects identified from SDSS spectra, J0905+2510 and J1637+3631. Reddening values were determined from the Pan-STARRS dust-maps (Chambers et al. 2016;Green et al. 2018). Photometric T eff were calculated from fitting our "LP 40−365 star" models to the GALEX NUV and SDSS photometry including the quoted reddening. Results from our template matching are given at the bottom. close agreement between the strong lines in the data and best-matching template indicates that this star also belongs to the class of "LP 40−365 stars." We also note some discrepancy near 5230Å, corresponding to a feature, currently not accounted for by the models, that is also observed in LP 40−365 as well as in J1603−6613 and LP 40−365 which Vennes et al. (2017) attributed to a resonance in the photoionisation cross-section of Mg. J0905+2510 has only a two-parameter astrometric solution (position) within the Gaia DR2 data. With G = 19.62, it is probable that a parallax will be available in a future Gaia data release. Assuming J0905+2510 has an absolute G magnitude similar to LP 40−365 and J1603−6613 (6.5 mag), the true distance is likely to be about 4 kpc from the Sun. Given the proper motion of (µ α cos δ, µ δ ) = (3.1 ± 3.4, 24.7 ± 3.0) mas yr −1 (Tian et al. 2017), the implied transverse velocity is about 500 km s −1 , which is broadly consistent with the 300 km s −1 radial velocity.
To more accurately measure the T eff , we fitted our LP 40−365 model grid to GALEX NUV and SDSS photometry (interpolating amongst models for estimating intermediate T eff to our grid-steps). We accounted for interstellar reddening, E(B − V) = 0.04 ± 0.02 mag, (as estimated via the Pan-STARRS 3D dust-maps at a heliocentric distance of at 4 kp; Green et al. 2018), obtaining 9840± 170 K. A full abundance analysis and comparison with the other members of this class necessitates higher-quality follow-up spectra.

J1637+3631
This star was first reported as a white dwarf with a C-rich, He-dominated atmosphere (DBQ spectral type) in the SDSS DR7 white dwarf catalogue of Kleinman et al. (2013). However, we neither see evidence for He nor C features. Beyond the general Rayleigh-Jeans slope, the main spectral features matched by the template (Figure 4, top) are the 4129Å Si ii and 4481Å Mg ii lines (rest wavelengths), although redshifted by about 300 km s −1 . Transitions from O i and Ca ii appear stronger than in the template, indicating somewhat different abundances than assumed for our grid of models. The 6402Å Ne i line is also tentatively detected in the unsmoothed SDSS spectrum (Figure 4, inset).
For J1637+3631, Gaia DR2 provides a full fiveparameter astrometric solution, although rather poorly constrained. With a parallax of 0.39 ± 0.60 mas, we place a 3 σ lower limit of about 500 pc on the distance. The implied absolute magnitude is M G = 11.8. Applying a smooth exponentially decreasing prior that accounts for the stellar density distribution in the Milky Way, Bailer- Jones et al. (2018) have estimated a distance of d = 1550 +1115 −630 pc, which implies an intrinsic magnitude of M G = 9.3 ± 1.2. Both these distance estimates place J1637+3631 within the white dwarf sequence of Figure 1.
Of course J1637+3631 could be much further away and still remain consistent with the Gaia parallax. If we assume J1637+3631 is a LP 40−365 star, and conservatively, like LP 40−365 itself, has M G = 6.5, the implied distance would be about 6 kpc. Given the 67 mas yr −1 proper-motion (Table 3), the transverse velocity would be ∼2000 km s −1several times larger than theoretical predictions for SNe Iax survivors (Shen et al. 2018b).
Hence, a more likely explanation is that J1637+3631 could be the second representative of the DOx white dwarf class, after J1240+6710 (Kepler et al. 2016). This hypothesis can be tested with medium-resolution follow-up spectroscopy. One difference we note between these two stars is the intensity of the Ca ii H+K lines: they are plainly visible in the spectrum of J1637+3631, but only an upper-limit of Ca was determined for J1240+6710. Our photometric T eff , determined from fitting GALEX and SDSS photometry (13 500 ± 600 K, Table 3), is cooler than the 21 590 ± 620 K measured by Kepler et al. (2016) for J1240+6710, which could explain the differences in Ca ii strength.

SPECTRAL ANALYSIS
3.1 J1603−6613 Fig. 2 shows that J1603−6613 has a spectrum remarkably similar to LP 40−365, thus providing a natural starting point for the spectral analysis. Owing to the high S/N and resolution of the X-shooter spectrum, we searched for lines from other elements in addition to those previously identified in LP 40−365. We identified new lines from Sc ii, Co i-ii, Zn i, Cu i, and Sr ii. C i lines were also found to be present at optical wavelengths. As with LP 40−365, no evidence for H or He were seen, but with more stringent limits from the much higher S/N ratio and slightly ≃ 1000 higher T eff . We precisely measured the blue-shift of J1603−6613 by fitting Lorentzian profiles to some of the strong, isolated O i and Mg i lines in the VIS arm of the X-shooter spectrum, obtaining v rad = −485±4 km s −1 , which confirms the estimate from the identification spectrum.
As previously noted in Paper I for LP 40−365, a hotter T eff clashes with the non-detection of He i absorption at 5877Å. Given the similarity between the two stars, we reconsidered the feasibility of a ONe-dominated atmosphere, as found by Vennes et al. (2017) for LP 40−365. Our model spectra follow the methods of Koester (2010), but we used a completely rewritten, modernised code, taking into account more recent atomic data 1 .
We initially approached the problem by fitting the Xshooter spectra alone. Given that O, Ne, and Mg, are the dominant elements and main contributors to the continuum opacity, their abundances are necessarily related to each other and coupled with T eff . Hence, due to the inherent inaccuracy of flux calibration and the ill-constrained effect of interstellar reddening, the resulting atmospheric fit was strongly degenerate, in a way that adjustments in T eff could be compensated for by changes in abundances, allowing for reasonable agreement with the data over a wide T eff -range. We therefore tried a different approach.
Instead of allowing T eff as a free parameter, we performed multiple least squares fits to the normalised spectra,   Table 1). Bottom panels: Residuals between observed and synthetic spectra, and observed and synthetic photometry. The best-fit model was obtained by minimising the residuals between observed and synthetic photometry, hence the slope observed for the observed vs. synthetic spectrum comparison. but at fixed T eff in steps of 100 K in the range 9500 K to 11 800 K. The optimal log g and abundances varied across the grid steps. We then used this model grid to fit SkyMapper and GALEX photometry (Table 1) with T eff , the interstellar reddening and stellar radius as free parameters, and the Gaia parallax and its error used as a Gaussian prior. While Gaia photometry were also available, their quoted uncertainties do not include systematic errors, which are expected to be significantly larger than this. Thus we chose to exclude these. The best-fit was achieved by calculating synthetic photometry in each bandpass and minimising the χ 2 ν . This resulted in a T eff of 10590 ± 370 K, a radius of 0.16 ± 0.03 R ⊙ . While the full line-of-sight reddening is about E(B − V) = 0.10, our fit indicated a 99 percent upper limit of 0.09 with a mean value of 0.03 in the posterior distribution, indicating the full line-of-sight extinction to be too high.
To obtain estimates for log g, abundances, and their errors, we fitted the T eff dependence of the element, Z, with a 2nd-order polynomial. The standard deviation of the scatter around the polynomial was taken as σ Z 0 , and the gradient determined at T eff = 10590 K. Finally, the uncertainty on Z was calculated as where σ T eff = 370 K. We show this graphically for Fe as an example (Fig. 6). We compare the observed spectrum and our best-fit model in Fig. 5. The atmospheric parameters and the element abundances are listed in Table 4.

J1825−3757
We confirmed the radial velocity of J1825−3757, obtaining v rad = −47 ± 3km s −1 from Lorentzian fits to isolated O i, Mg ii, and Ca ii lines of the X-shooter spectrum.
Once again we searched the spectrum to identify all elements present in the stellar atmosphere. As inferred from the bluer B p − R p colour, J1825−3757 is hotter than both LP 40−365 and J1603−6613, and so the ionisation balance for each element is expected to be different. Apart from Cu (the higher T eff does not permit detection of Cu i, and no strong optical lines of Cu ii are expected), all elements found at LP 40−365 and J1603−6613 were re-identified. We note that several other elements populate higher ionisation states than found for the other two stars: Al iii was identified from doublets around 3610Å and 5710Å, Si iii from a doublet around 4560Å, O ii at 4350Å and 4650Å, and S ii from a multitude of lines.
For the spectral analysis, we took a similar approach to that described for J1603−6613. We built a model grid between T eff 11 800-14 000 K, in 100 K steps, with the surface gravity and abundances as free parameters. For the photometric fit we made use of the SkyMapper and 2MASS data (Table 1). We excluded GALEX photometry from our fit, because the NUV spectral region is expected to be more heavily line blanketed than found for the optical, and we cannot rule out unaccounted opacity sources, e.g. from elements presently unidentified. Using the Gaia parallax (Table 1) as a prior, we found T eff = 12830±450 K, E(B−V) = 0.080±0.017, and R = 0.60 ± 0.03 R ⊙ .
The values and uncertainties on log g and abundances were calculated in the same manner as described for J1603−6613, and are given in Table 4. Note that, because of the higher T eff , the upper-limits for hydrogen and helium are far more constraining than for LP 40−365 and J1603−6613. This strongly affirms the ONe-dominated nature of these stars. We compare the observed spectrum and our best-fit model in Fig. 7.
A major difference in the results of J1825−3757 vs. LP 40−365 and J1603−6613 is the surface gravity, which is more than one order of magnitude lower. This is not too surprising given the much larger radius we found for the photometric fit, implying the star is more inflated compared to the other two. This is discussed further in Section 5.1.3.
As a final comment, we examine the properties of the updated atomic data that have a more tangible effect in modelling the spectrum of J1825−3757. While our line list already included many thousands of transitions from the NIST and VALD3 databases, we noticed that the profiles of some of the weakest O i and Mg ii lines were poorly reproduced by our models. In particular, three moderately strong Table 4. Composition our best-fitting synthetic spectra of LP 40−365 (Section 3.1), J1603−6613 (Section 3.2), and J1825−3757 (Section 3.3). Abundances are expressed as the number fraction of the total, log (Z/Tot), where the contribution from is not included.

Parameters
LP 40−365 J1603−6613 J1825−3757 9800 ± 350 10 590 ± 370 12 830 ± 450 log g [cgs] 5.5 ± 0.3 5.34 ± 0.20 4.21 ± 0.18 10 −8.14 ± 0.29 −8.40 ± 0.25 lines at 4014, 4094, and 4110Å (rest-frame air wavelengths) were absent in our model. From their intensities and large line widths we concluded these are most likely higher members of the Mg ii series transitioning from the closely spaced 4d/f (93 311/93 800 cm −1 , respectively), as confirmed from wavelengths tabulated by Moore (1959). We found the three missing lines (and some higher members of the series that are less obviously present) in Robert Kurucz's line archive, which includes not only wavelengths and oscillator strengths, but also Stark, quadratic Stark, and van der Waals broadening constants. This update in our models shows a close agreement with the observed spectrum. However, we note that for the highest series members, the strengths are systematically under-predicted, likely owing to the proximity of the upper-levels to the ionisation limit, adversely affecting their estimated occupation numbers. The updated Mg ii and O i lines are generally weaker in the spectra of J1603−6613 and LP 40−365, because of their lower T eff . Hence, their inclusion in the models leads to a minor improvement. Although most of the strong lines are well determined in the literature, our concluding remark is that, because of their extreme photospheric composition, stars in the LP 40−365class necessitate complete and accurate atomic data even for lines usually considered astrophysically irrelevant.

LP 40−365
We re-analysed the prototype star of this class, combining our new HST/STIS spectrum (Fig. 8)   The red-shifted spectrum has been corrected into the laboratory rest-frame. We note the C i detection and the strongest Mg and Si lines, which cause most of the opacity in the spectrum.
and Asiago 1.82-m Copernico/AFOSC spectra from Paper I. This combined, flux-calibrated spectrum has an average spectral resolution of ≃ 2-4Å in the optical and NUV ranges, respectively. The 5350-5700Å region, which is only covered by the Copernico/AFOSC data, has a 13Å resolution. We scaled the optical spectrum to the Pan-STARRS magnitude (g = 15.635; Flewelling et al. 2016). Although the HST/STIS spectrum seamless matches the optical fluxes, we note it is 4 σ brighter than the GALEX NUV magnitude. The optical spectrum, instead, appears slightly bluer than the published Pan-STARRS and Gaia DR2 photometry, corresponding to a small colour term in the range of B − V = −0.02, which is comparable to the photometric uncertainties. In joining the spectra together, we transposed the optical region onto a vacuum wavelength scale to match the HST/STIS data (using the relation given in Morton 2000). We corrected the combined spectrum for radial velocity shift of LP 40−365 (v rad = +500 km s −1 ), which is confirmed by the redshift of the C i line at 1930Å and other strong Mg i-ii lines. Finally, we corrected the combined spectrum for the total interstellar reddening of E(B − V) = 0.026 (Schlegel et al. 1998), corresponding to A g ≃ 0.09 and A NUV ≃ 0.22 magnitudes of extinction in the Pan-STARRS g and GALEX NUV bands, respectively, as estimated via the Fitzpatrick (1999) standard R V = 3.1 extinction law. The combined dereddened spectrum is displayed in Fig. 9.
From an initial appraisal of the combined spectrum, the NUV data indicate a hotter temperature with respect to the optical spectrum, for which we inferred T eff = 8900 K using a He-dominated atmosphere in Paper I. In the previous sections, we have shown the high-quality X-shooter spectra of J1603−6613 and J1825−3757 are well matched by synthetic spectra with ONe-dominated atmospheres. Especially from the analysis of J1825−3757, we have excluded the presence of He with a very stringent limit. Hence, these results favour an ONe-dominated atmosphere also for LP 40−365, as found by Vennes et al. (2017). As we did for J1603−6613 and J1825−3757, we built a grid of synthetic models in the T eff = 8900-10 200 K range, with O, Ne, and Mg as dominating species. In view of the complex nature of the problem, we approached the analysis of LP 40−365 by identifying a set of key spectral features to compare against, in order to evaluate the goodness of fit. These are: i) the level of NUV flux; ii) the spectral lines of Ne, O, and Mg; iii) the continuum slope in the optical range; iv) the Mg i jump near 3800Å; and v) the ionisation equilibrium of Mg i/Mg ii. Because we did not find a solution that satisfies all the above criteria, we have relied more on the flux-calibrated NUV spectrum, obtaining a best-fit synthetic spectrum with T eff = 9800 K and log g = 5.5 that is shown in Fig. 9. This model approximately reproduces the NUV flux between 2200-3000Å, the O i, Ne i, and Mg i lines, and the Mg i jump. The main defect of the synthetic spectrum is that it does not accurately reproduce the ionisation equilibrium of Mg i/Mg ii, over-predicting the intensity of Mg ii lines as shown in the top-centre panel of Fig. 9. In addition, we also note that the residuals between the observed and synthetic spectra, displayed in the lower panel of Fig. 9, suggest that the optical continuum of the synthetic spectrum is slightly too steep, whereas the flux below 2200Å is ∼ 30 per cent too low. As noted in Paper I, models with increasing lower T eff have stronger Ne i lines that overpredict the observed ones.
Based on the strong coupling between O, Ne, and Mg, we estimated a 1σ T eff uncertainty of 350 K, beyond which the absorption caused by these elements becomes too strong or too weak. The surface gravity cannot be constrained more precisely than log g = 5.5±0.3 dex, which is compatible with a ≃ 0.16 R ⊙ radius, as derived from the Gaia parallax and Pan-STARRS g magnitude, and by assuming a plausible mass between 0.15-0.60 M ⊙ (a more detailed discussion, with reference to Paper I, is given in Section 5.1.2).
Although the spectral analysis of LP 40−365 is still plagued by systematic uncertainties, the HST/STIS observations have enabled a better overview of its composition, clarifying the incompatibility between a He-dominated atmosphere proposed in Paper I and an ONe-dominated atmosphere, as proposed by Vennes et al. (2017). The HST/STIS data favour the latter, and our new best-fit model has a T eff that is now more compatible with that estimated by Vennes et al. (2017). However, we note that we have measured the dominant elements in the ratio of O:Ne:Mg ∼ 34:54:6, in contrast with their O:Ne ∼ 50:50 estimate.
The near-UV spectrum has allowed us to detect the C i 1930Å line, enabling us to constrain the abundance of this element. Thus, we have determined abundances for 14 elements (C, O, Ne, Na, Mg, Al, Si, S, Ca, Ti, Cr, Mn, Fe, Ni). The abundances of detected trace elements were measured by iteratively altering them, one by one, and keeping the abundances of the other elements fixed. The mean values and errors are established from the means and standard deviations derived by fitting the strongest lines of each element. We have confirmed a remarkable similarity between element abundances for the three studied LP 40−365 stars (more details are given in Section 5.1.1). We could not estimate plausible systematic uncertainties of trace elements due to the previously mentioned strong coupling of O, Ne, and Mg with T eff , which would cause them to be much larger than statistical errors and difficult to be properly quantified.
We have also estimated upper limits on H, He, P, V, Sc, Co, Cu, Zn, and Sr. We note that Sc, Co, Zn, and Sr are detected in J1603−6613 and J1825−3757, and Cu is detected  Table 4. To conclude, as discussed in the previous sections and in Paper I, it is likely that important physics may be still missing from our models. In addition to the already mentioned atomic line data, we note the lack of absorption coefficients for the bound-free and free-free contribution of Mg − , Al − , and Si − ions could also explain unidentified spectral features.

KINEMATIC ANALYSIS
Following the procedure described in Paper II for LP 40−365, we used Gaia DR2 astrometry and v rad to trace, in a probabilistic fashion, the motion of J1603−6613 and J1825−3757 within the Milky Way potential. Using the python package for galactic dynamics, galpy (Bovy 2015), we modelled the Milky Way's potential with the standard module MWPotential2014, which includes a power-law density profile with exponential cut-off for the bulge, a Miyamoto-Nagai disc, and a dark matter halo (see details in Bovy & Rix 2013;Bovy 2015).
We imposed a Galactic disc rotation of V c = 239 ± 9 km s −1 , with the Sun located at R 0 = 8.27 ± 0.29 kpc (Schönrich 2012), accounting for a one-to-one correlation between R 0 and V c , and a peculiar motion as determined by Schönrich et al. (2010). We sampled 10 000 boundary conditions with a Monte Carlo method, by assuming Gaussian distributions for the Galactic parameters, the radial velocities, and the astrometric parameters of the two stars. For the latter we took into account the full covariance matrix as delivered in Gaia DR2. The orbital parameters of J1825−3757 and J1603−6613 are summarised in Table 5, where we also list those of LP 40−365 for comparison.
Given the lack of or the low precision of Gaia parallaxes for the two SDSS stars, J0905+2510 and J1637+3631, we did not compute their detailed orbits. Nevertheless, we note that the radial and transverse velocities of J0905+2510, which we suggest to be another LP 40−365 star (see Section 2.5 for details), deliver a rest-frame velocity of v rf ∼ 730 km s −1 , making it likely unbound from the Milky Way.

The unbound trajectory of J1603−6613
Although the parallax of J1603−6613 is less precise than that of LP 40−365, we can confirm it as unbound from the Milky Way with a rest frame velocity of v rf = 805 +49 −32 km s −1 , exceeding the escape velocity of v esc = 550 km s −1 estimated with MWPotential2014 at the corresponding Galactocentric radius of R G = 7 kpc (cf with ≈ 520-530 km s −1 in the Solar neighbourhood; Piffl et al. 2014;Williams et al. 2017). We show the Galactic trajectories of the three LP 40−365 stars with precise Gaia parallaxes in Fig. 10. J1603−6613 is escaping from the Milky Way in the southern Galactic hemisphere, while LP 40−365 is departing from the disc at positive latitudes. J1603−6613 has not yet reached its closest approach to the Sun; its spectrum is blue-shifted. Both stars move along the direction of the Galactic rotation, diverging from the tangential vectors by ∼ 20-26 deg on average. They also have grazing trajectories with respect to the plane (≃ 5 and 10 deg for LP 40−365 and J1603−6613, respectively), which imply a significant contribution of the Galactic rotation to the rest-frame velocity.
We note that the two unbound stars may have crossed the Galactic plane at mutually close locations, as displayed in Fig. 10, although with different flight times from the plane, τ fl (Z = 0). We computed the mutual separations between the past trajectories of LP 40−365 and J1603−6613 as a function of time, finding that the minimum separation of 1 kpc might have occurred 1 Myr ago. The orbit integration confirms that their similar spectral appearance is a characteristic of the formation channel and its not related to a specific event that produced both stars at the same time.

J1825−3757 is on a Milky Way bound orbit
Combining the Gaia parallax and proper motions with the radial velocity of J1825−3757, we estimate a rest frame velocity of v rf = 408 +39 −34 km s −1 , which is ≃ 4σ slower than v esc ≈ 560 km s −1 estimated at R G = 7.4 kpc. This implies that J1825−3757 is currently gravitationally bound to the Milky Way. For this star, we followed the orbit evolution up to 10 Gyr from now, in order to estimate the average eccentricity (e), the pericentre and apocentre, the azimuthal period, and the escape probability, which we estimate as the fraction of orbits reaching R G = 100 kpc. These and other relevant kinematic parameters are listed in Table 5. We note that J1825−3757 is relatively close to the pericentre of the predicted orbits, and it is moving towards low Galactic coordinates, i.e. negative Z. It goes around the Milky Way once every ≃ 730 Myr, moving on an eccentric (e = 0.69), halo-like orbit against Galactic rotation. If J1825−3757 had previously crossed the Galactic plane, it would have occurred ≃ 1.5 Myr ago. Its relatively large but negative vertical component of the angular moment, J Z = −2800 kpc km s −1 , indicates its orbit is rather flat (Z max ≃ 7 kpc), but it extends to large Galactic radii (apocentre at R G = 39 +31 −14 kpc). A small fraction of our simulated orbits (< 10 per cent) imply a possible escape of J1825−3757 from the Milky Way in the next 10 Gyr.

DISCUSSION
The evidence presented so far demonstrates that LP 40−365, J1603−6613, J1825−3757, and likely J0905+2510, have unique -but mutually similar -spectral characteristics and peculiar kinematics, which make them clearly distinct from other classes of stars. Our new observations strongly support the interpretation that LP 40−365 and the new stars are the partly burnt white dwarf accretors that survived disruption from a thermonuclear supernova in a singledegenerate scenario (possibly a SNe Iax, as initially advanced by Vennes et al. 2017, and further reiterated in our Paper I; Paper II). In the following sections, we will discuss the properties that identify these stars as members of what appears to be a class of "LP 40−365 stars." We will focus on the three stars with precise Gaia parallaxes for which we performed a detailed spectral analysis, but our conclusions can be extended to J0905+2510, which just awaits Gaia to measure its parallax and proper-motions as well as a higherquality spectrum for a detailed abundance analysis.  Table 4, are displayed relatively to Fe and normalised to the Solar composition (Asplund et al. 2009, dashed line about 0). Black circles, blue squares, and red diamonds represent LP 40−365, J1603−6613, and J1825−3757, respectively. When error bars are not visible, their size is comparable to or smaller than the symbol. Upper limits are shown as downward pointing arrows.

Atmospheric composition
The comparison among the atmospheric composition of the three LP 40−365 stars is displayed in Fig. 11, where we plot the number abundance of each detected element, Z, as [Z/Fe] = log (Z/Fe) − log (Z/Fe) ⊙ , where the Solar reference abundances are derived from Asplund et al. (2009). Two elements, H and He, for which we measure upper limits (Table 4), are excluded from this figure. The remarkable similarity in the average abundance pattern of LP 40−365 stars is confirmed by a typical dispersion in abundance of order 0.25 dex, which is similar to the size of uncertainties, for most of the detected elements. We especially note the characteristic α enhancement and moderately super-Solar abundances of the iron-peak elements Cr, Mn, and Ni. One element that is detected in both J1603−6613 and J1825−3757, Sc, shows the largest mutual scatter (≈ 2 dex). Co and Zn, also detected in J1603−6613 and J1825−3757, have the largest super-Solar abundances among the iron-group elements. Cu was marginally detected in J1603−6613, also with a super-Solar abundance.
Our revised analysis of LP 40−365 with respect to Paper I now indicates a close-to-Solar Mn abundance, [Mn/Fe] = 0.07. Although the uncertainties on the abundances of LP 40−365 are relatively larger than those estimated for the other two stars (cf Section 3), it appears as the least Mn-rich of the three, with J1603−6613 and J1825−3757 having [Mn/Fe] = 0.43 and 0.63 dex, respectively. Given the super-Solar trend, we confirm the evidence in favour of Mn-enhancement. Three other elements, Cr-, Co-and Ni-to-Fe confirm a similar Super-solar trend, which strengthen the suggested connection near-M Ch thermonuclear explosions (Paper I). In fact, we note that nucleosynthesis of these elements is known to be enhanced by electron capture that is thought to occur when the conditions for nuclear statistical equilibrium are reached at the high central densities of near-M Ch explosions (ρ > 2 × 10 9 g cm −3 ), which are expected to occur through mass-accretion from a nondegenerate companion, and not by exploding a sub-M Ch of 0.8 M ⊙ white dwarf through a merger (Iwamoto et al. 1999;Seitenzahl et al. 2013).
Super-Solar detection of Cu and Zn also agree with the observations of normal stars, for which the abundance pattern of these elements is proposed to need a significant contribution from thermonuclear supernovae that previously enriched the interstellar medium with their nucleosyntheis yields (Matteucci et al. 1993). The close-to-Solar detection of the p-nucleus, Sr, in J1603−6613 and J1825−3757 is also interesting, as this heavier element is expected to undergo significant production in the external layers of white dwarfs during thermonuclear explosions (Travaglio et al. 2011). Finally, we also note that ratios of iron-peak elements (e.g. Mn-to-Fe, Ni-to-Fe, and Mn-to-Cr mass fractions) deserve further investigation, as they could be age-estimators for the progenitors of LP 40−365 stars (e.g. via element-ratio vs metallicity relations like those employed for the characterisation of supernova nebular remnants; Badenes et al. 2008;Yamaguchi et al. 2015).

Comparison with theoretical yields
In Paper I, we suggested that LP 40−365 was unlikely to be the donor star in a binary supernova because its α element and iron-peak rich atmosphere is incompatible with the pollution from nucleosynthesis yields of both thermonuclear and core-collapse supernovae. We observed a better correspondence with the results of three-dimensional hydrodynamic simulations of pure deflagrations by Kromer et al. (2013) and Fink et al. (2014), although with some differences. These pure deflagration models consider a CO white dwarf that explodes near M Ch after accreting mass from a non-degenerate companion. The deflagration leads to a partial burning of the accreting CO white dwarf, which becomes what is defined as a "bound remnant." The Fink et al. (2014) simulations take into account a range of internal densities and explosion intensities, of which we selected those forming bound remnants of ≈ 0.1-0.55 M ⊙ , based on the observation that the LP 40−365 stars should be low-mass objects (see Section 5.1.3).
Although the models that form low-mass bound rem-  Fig. 11, is compared to the bulk-composition of bound remnants resulting from threedimensional hydrodynamic simulations of pure deflagrations of CO (Fink et al. 2014) and ONe white dwarfs ), plotted as orange and blue curves, respectively. The circles, squares, and triangle indicate whether three, two, and one star, respectively, have been used to compute the average abundances for a given element. The error bars represent the standard deviation of the variance for each element. Upper limits are not taken into account.

C O NeNaMgAl Si S Ca Sc Ti CrMnCo Ni CuZn Sr
nants are disfavoured by Fink et al. (2014) by comparisons with observed SN Iax lightcurves and spectra, because of their high energy output, they seem to deliver the best match to the average composition of the LP 40−365 stars (Fig. 12). The best match is obtained by a 0.1 M ⊙ bound remnant, which corresponds to a progenitor with the highest central density (5 × 10 9 g cm −3 ). However, we note that the observed abundances of the main atmospheric constituents (i.e. O, Ne, Mg) exceed theoretical predictions by 0.5-1.5 dex. In contrast, C is at least 1.5 dex less abundant than the theoretical curves. The number abundances of trace elements, excluding Sc, Co, Cu, and Zn, are roughly within 1 σ of the locus defined by bound-remnant models. In addition to this comparison, we also examined the recent results obtained by Jones et al. (2016Jones et al. ( , 2018 for mass-accreting ONe white dwarfs. Whether accretioninduced collapse can be avoided is still subject to discussion (Brooks et al. 2017b;Schwab & Akira Rocha 2018). However, the three-dimensional simulations performed by Jones et al. (2018) have shown that thermonuclear explosions in electron-capture supernovae (tECSN) could take place under certain conditions, leading to the formation of partly burnt bound remnants instead of neutron stars. In their work, Jones et al. (2018) already noted an overwhelming discrepancy with the abundance of α elements measured for LP 40−365, in contrast to a good agreement with the iron-peak elements (Cr, Mn, Ni). Their result is confirmed by the inclusion of the new LP 40−365 stars in Fig. 12. Of the new detected elements, Sc, Co , Cu, and Sr, are discrepant with respect to the tECSN model. Zn, however, displays a good match.
To further investigate the comparison with theoretical bound remnants, we plotted histogram bars representing the mass fraction of detected elements (Fig.13). The remarkable similarity among LP 40−365 stars is once again evident, with Ne as the most abundant element by mass fraction (59-65 per cent), followed by O (29-31 per cent), and Mg (3-9 per cent). The remaining elements are just 1-2 per cent of the total mass. The results from theoretical simulations by Fink et al. (2014), instead, contain between 25-40 per cent of C and about 45 per cent of O. These fractions are slightly dependent on the explosion intensity and quite sensitive to central density of the accreting white dwarf. On the other hand, the bound-remnant simulated by Jones et al. (2018) contains 45 per cent of heavy metals, most of which are Fe and Ni.
Whether LP 40−365 stars formed as the result of SNe Iax or tECSN, it would be necessary that either Cburning was more intense in a CO white dwarf or, vice-versa, less C, O, and Ne were burned in a ONe white dwarf. In both cases, we stress that theoretical results for the bulk composition of bound remnants might not reflect the atmospheric abundances, given that internal mixing due to gravitational settling, diffusion, and convection should definitely have a strong impact on their internal structure, atmospheric composition, and spectral evolution. Furthermore, we note that the remarkably similar atmospheric composition of the three LP 40−365 stars contrasts with the large difference in radii between J1825−3757 and the other two objects. Two simple interpretations are that either the three stars are fully mixed or their atmospheres may vary in thickness. For example, heavy elements are theoretically expected to be mostly confined in the stellar interiors. Hence, the Jones et al. (2018) bound remnants could end up representing better the atmospheric composition of LP 40−365 stars (this case would also compatible with the predicted existence of Fe-core white dwarfs; Isern et al. 1991). On the other hand, the Fink et al. Finally, we speculate that a compromise between low C and high O, Ne, and Mg, abundances might be achieved if the progenitors of LP 40−365 stars were hybrid CONe white dwarfs (as first proposed by Vennes et al. 2017), the existence of which is still debated (Denissenkov et al. 2013(Denissenkov et al. , 2015Brooks et al. 2017a). One-dimensional computations performed by Bravo et al. (2016) showed that thermonuclear explosions of such stars could form low-mass bound remnants that might represent LP 40−365 stars. However, we also note that three-dimensional hydrodynamic simulations of pure deflagrations of hybrid CONe white dwarfs suggest that bound remnants would be typically more massive Ne O Mg C others Figure 13. The atmospheres of LP 40−365 stars compared to the predicted bulk composition of three representative boundremnants resulting from the pure deflagrations of CO white dwarfs (Fink et al. 2014) and an ONe white dwarf ). The 0.19 M ⊙ and the 0.09 M ⊙ remnants result from the same initial CO white dwarf model, but the former experiences a weaker deflagration with respect to the latter. The 0.09 M ⊙ and 0.10 M ⊙ bound-remnants, instead, experience deflagrations of equal intensity, but the former has an initial central density that is ≃ 2 times smaller than the latter. The ONe white dwarf leading to the 0.44 M ⊙ remnant has a similar structure as the CO white dwarfs, with Ne in place of C, and external stratified shells of C, He, and H.

Masses and radii
Given the new spectral analysis of LP 40−365 presented here, we have revised our previous estimates of its physical parameters. We estimated the radius as R = (1/̟) f /F = 0.16 ± 0.01 R ⊙ , where f and F are the extinction corrected flux and the integrated Eddington flux, respectively, both considered in the Pan-STARRS g-band. The luminosity is L = 0.20 ± 0.04 L ⊙ , as determined via the Stefan-Boltzmann law. The radius is 11 per cent smaller than what estimated in Paper II, compensating for the 10 per cent increase in T eff and a more realistic spectral modelling. Considering the log g = 5.5 ± 0.3, we inferred a mass of M = gR 2 /G = 0.28 +0.28 −0.14 M ⊙ (24 per cent smaller than that of Paper II), where G is the gravitational constant.
The physical parameters of the three stars are listed in Table 6. Physical parameters of the LP 40−365 stars. The 1σ errors and the 5-95 per cent ranges were estimated via Monte Carlo sampling. The apparent fluxes are listed without errors, because they carry very small errors from precise broad-band photometry and we did not include the effect of T eff uncertainties in the error budget. The errors on radius, luminosity, and mass, account for those on parallax, T eff , and log g.   Boyajian et al. 2013), and white dwarfs in the Solar neighbourhood (white circles; Giammichele et al. 2012). We note that the 1 σ errors on the estimated physical parameters are smaller or of the order of the symbol sizes (on the adopted log-log scale). For reference, we draw the canonical-mass white dwarf cooling sequence (beige-coloured strip Fontaine et al. 2001), the 100 Myr old main-sequence (red curve; Choi et al. 2016), and the evolutionary tracks for 0.47, 0.48, and 0.50 M ⊙ hot subdwarfs (blue curves; Dorman et al. 1993) The luminosities for given stellar radii are shown as a function of T eff (dotted curves). Table 6, which includes both the 1 σ uncertainties and the 5-95 per cent ranges.

Parameters
LP 40−365, J1603−6613, and J1825−3757 are plotted in the theoretical HR diagram (Fig. 14), where we include other classes of stars for comparison. While their T eff span just ∼2500 K, their luminosities cover two orders of magnitude (0.1-10 L ⊙ ). This large portion of the HR diagram is not occupied by any long-lived phase of stellar evolution, but it is occasionally scattered with (pre-) low-mass He-core white dwarfs, which are the products of binary interaction (Istrate et al. 2016), metal-poor A/F-type stars (Pelisoli et al. 2019), and short-lived evolutionary phases with stochastic behaviour, like cataclysmic variables or millisecond pulsar companions. J1603−6613 sits very close to LP 40−365, as expected for two spectroscopically and photometrically similar stars. Although J1825−3757 appears closer to the zero-age extreme horizontal branch and the main sequence (note the log-log scale), it is ∼ 15 000 K cooler but larger than equal-luminosity hot subdwarfs. Hence, it is at least an order of magnitude less luminous than equal-T eff main sequence stars.
Given their unusual, "sub-luminous" location in the HR diagram and their masses, which are below the typical mass of canonical white dwarfs (e.g. 0.621 M ⊙ ; Tremblay et al. 2016), we speculated that the interiors of LP 40−365 stars may not be fully degenerate. We measured their average densities as ρ ≃ 340, 240, and 4.5 g cm −3 for LP 40−365, J1603−6613, and J1825−3757, respectively. Hence, we note that LP 40−365 and J1603−6613 have average densities that are comparable to those of brown dwarfs (or hot subdwarfs), whereas the mean density of J1825−3757 is similar to that of the Sun.
Finally, as discussed in Paper I, we remark that the radii we measure for LP 40−365 stars are in agreement with the observed low-rotation rates. We note that, even if their progenitors were tidally locked in a short period binary (see Section 5.1.5; Fuller & Lai 2012), we expect that mass loss and radius expansion caused by the supernova explosion could have both contributed to slow down a potentially fastrotating star.

Evolutionary timescales
Given their rarity, the current location of the LP 40−365 stars in the HR diagram of Fig. 14 is likely transitory and it is probably characterised by a relatively fast evolution, which could be shorter than the Galactic crossing time for the unbound stars (∼140 Myr; Paper II). Shen & Schwab (2017) modelled the first ∼1000 yr of evolution for post-SNe Iax bound remnants. In this early phase, the stellar photosphere is initially inflated to > 1000 R ⊙ , the energy output is dominated by strong winds powered by the decay of radioactive 56 Ni and 56 Co, and the release of energy takes place over ≃10-100 yr timescales. More recent work by Zhang et al. (2019) covers the longterm evolution of post-SN Iax bound remnants (named postgenitors by the authors), specifically discussing LP 40−365 as one such star. These authors proposed that post-SN Iax bound remnants would initially cool down over relatively short time scales, and then heat up and brighten again over an at least 10 Myr-long timescale until they reach the white dwarf cooling sequence again. This mechanism would be regulated by convection and diffusion of iron-peak elements, which initially would settle in the atmosphere of the bound remnants for 10,000 yr during the early phase and would cause a large opacity. The long-lived brightening phase would be caused by a sharp decrease of ironpeak element abundances when the bound remnant becomes cool enough, characterised by a turn-off at the faintest and coolest point that depends on the mass and envelope size of the bound remnant.
Despite The presence of iron-peak elements could indicate that the decline of stellar luminosity envisaged by Zhang et al. (2019) is characterised by much longer timescales. Hence, the relative ages of the three LP 40−365 stars could be inverted, i.e. J1825−3757 and would be the youngest of the three objects. For quantifying this, we estimated the Kelvin-Helmholtz timescales, τ KH , for the three LP 40−365 stars that, in absence of internal sources of energy production, is a good estimator for the evolutionary timescales. Adopting the formulation τ KH = 3GM 2 /7RL, which holds for a wide range of objects with internal structure reproducible by a polytrope with index n = 3/2, and using the values from Table 6, we obtain τ KH ∼ 32, 12, 0.1 Myr for LP 40−365, J1603−6613, and J1825−3757, respectively. The τ KH of LP 40−365 and J1603−6613 are in the same range of the evolutionary age estimated by Zhang et al. (2019). The τ KH of J1825−3757 suggests a much younger age that is, however, larger than the very short-lived phase of the Zhang et al. (2019) models (days-weeks). We note that such a short τ KH is still compatible with the non-association of J1825−3757 with known diffuse nebulosity caused by young supernova.

Kinematics and birth places
The rest frame velocities of LP 40−365 and J1603−6613 (v rf ≃ 800-850 km s −1 ) are compatible with the two stars being unbound from the Milky Way. This condition implies that, in absence of a powerful birth kick, both stars took advantage of a large boost from the Milky Way rotation so that v 2 rf = V 2 c + v 2 ej +2V c v ej cos α cos β = 800-850 km s −1 , where α and β are the angles drawn by the ejection velocity vector, v ej , with respect to the vertical direction and the Galactic rotation, respectively. Considering the almost-parallel alignment of the trajectories of LP 40−365 and J1603−6613 in the X − Y plane (5-10 deg with respect to direction of the Galactic rotation in the plane, and 20-26 deg with respect to the vertical direction), we found the minimum possible ejection velocities in the range of v ej ≈ 550-600 km s −1 .
In order to achieve the smallest v ej possible, obtainable with the maximum benefit from V c , the supernovae should have exploded at relatively close distances from the Galactic plane and not to far away from the centre (we note that V c is already ∼10 per cent slower above |Z | = 1 kpc, and it stays constant up to R G ≃ 25 kpc; Williams et al. 2013;Huang et al. 2016). This condition is compatible with the prediction for runaway stars ejected by a disc-like distribution of supernova progenitors (Kenyon et al. 2014).
Adopting the Zhang et al. (2019) age estimate of 23 Myr as the flight time for both LP 40−365 and J1603−6613, we could track their simulated trajectories back to the fourth and third Galactic quadrants, respectively. We note that the R G ∼ 25 kpc estimated by Zhang et al. (2019) does not account for the past trajectory of LP 40−365, and it is rather equivalent to the distance covered in 23 Myr at 850 km s −1 . Hence, we find that LP 40−365 and J1603−6613 would have been ejected from possible birth sites at R G ≈ 15 kpc and Z ≈ −1.2 kpc, and R G ≈ 17 kpc and Z ≈ +3.4 kpc, respectively. These coordinates are roughly compatible with birth sites within the thick disc. Furthermore, we note the Galactocentric radii of the identified birth sites are close to the empirical determination for the edge of the Galactic disc, i.e where most of the star formation seems to take place (within a R G = 13-16 kpc cut-off, as observed from photometric surveys of the Milky Way Sale et al. 2010;Minniti et al. 2011), and thus, they are in agreement with SN Iax observations, which are mostly seen in association with regions of recent or ongoing star formation (Foley et al. 2013;Jha 2017;Lyman et al. 2018). If we wanted to impose birth sites that are nearer to the plane, |Z | < 1, to find a better match with the minimum v ej = 600 km s −1 , we would obtain R G = 12, 8 kpc and τ fl = 19, 7.5 Myr for LP 40−365 and J1603−6613, respectively.
We previously noted that J1825−3757 is bound to the Milky Way. The fact that it moves in the opposite direction with respect to the Galactic rotation and it follows a fairly eccentric orbit strongly supports the supernovaejection mechanism. Observing that the direction of motion of this star has an equally small inclination with respect to the Galactic plane (a few degrees) as the other two, and it is a ∼165 deg with respect to the Galactic rotation, we estimated a minimum ejection velocity of v ej ≃ 600 km s −1 , similar to the other two LP 40−365 stars. Finally, we note that if J1825−3757 is as young as τ KH ∼ 0.1 Myr, it would have formed within the thin disc.

Binary progenitors
Given that we find similar v ej for the three LP 40−365 stars, it is reasonable to think they formed from the disruption of similar binary progenitors. In the ideal case of instantaneous mass loss and no interactions between accretor and donor stars (Hills 1983), v ej ∼ 550-600 km s −1 is a lower limit on the orbital velocity at the moment of explosion. This condition could be met if the LP 40−365 stars were white dwarfs with ≃ 0.8 M ⊙ He-burning donors in a 1 hr orbit (Paper II). This result is compatible with a single-degenerate scenario proposed for SN Iax (including progenitors with CO and hybrid CONe cores; Wang & Han 2009;Wang et al. 2013Wang et al. , 2014, and may also work for ONe white dwarfs. We note that kicks due to asymmetric explosions and interactions with the both the donor star and supernova ejecta are expected to take place (Marietta et al. 2000;Jordan et al. 2012;Pan et al. 2012;Liu et al. 2013), adding complexity to this interpretation.
The rest frame velocities of the unbound LP 40−365 stars contrasts with the much larger v rf 1000 km s −1 of the other known hyper-runaway stars, i.e. the D 6 stars (proposed as the former donor white dwarfs in double-degenerate supernova progenitors; Shen et al. 2018b) and the hot subdwarf star US 708 (proposed to be the former low-mass He subdwarf donor in a single-degenerate sub-M Ch supernova; Geier et al. 2015). Hence, the kinematics of LP 40−365 stars confirms their origin from a distinct class of binary progenitors, also exploding as thermonuclear supernovae, but possessing smaller orbital velocities than the D 6 and US 708 progenitors.

Prospects for end-of-mission Gaia detections
Taking into account their models for the long-term evolution of post-SN Iax bound remnants, Zhang et al. (2019) have estimated that four such stars would be detectable within 2 kpc from the Sun. Shen et al. (2018b) estimated the number of D 6 stars in Gaia as 14-40 within 1 kpc. Their results are comparable, considering that empirical SN Iax rates are typically estimated to occur as 1-30 per cent of the SN Ia rates (Foley et al. 2013;Li et al. 2011a;Graur et al. 2017;Jha 2017). Although Zhang et al. (2019) considered LP 40−365 stars as good candidates for the post-SN Iax bound remnants, they assumed these runaways to have v rf ∼ 1000 km s −1 . Such large velocities have only been measured for D 6 stars and the hot subdwarf US 708. Given the smaller v rf of LP 40−365 stars, likely due to their different formation mechanism, we re-estimated their space density.
Assuming the progenitors of LP 40−365 stars have a space density reflecting that of the Galactic thin disc, we distributed them according to an exponentially decaying profile (scale-height and scale-length of 0.3 and 3 kpc, respectively). We scaled the formation rate of LP 40−365 stars to 10 per cent of the Galactic SN Ia rate of 10 −13 yr −1 M −1 ⊙ (Li et al. 2011b), as adopted by Shen et al. (2018b). This rate is compatible with the average of other estimates from the literature, which include a mix of progenitors from both the single-and doubledegenerate channels (i.e. SNe Ia with a variety of delay times; Scannapieco & Bildsten 2005;Mannucci et al. 2005Mannucci et al. , 2006Aubourg et al. 2008;Ruiter et al. 2009;Maoz et al. 2011). Hence, considering a Galactic disc mass of ≃ 5 × 10 10 M ⊙ (Bland-Hawthorn & Gerhard 2016), we created 50,000 objects with birth-times that were uniformly distributed over a time interval of 100 Myr, t 0 = U(0, 100). We assigned Gaussian-distributed initial velocities, v ej = N (600, 25), which we modelled on those we have inferred for the observed LP 40−365. Finally, we used galpy to estimate the kinematic parameters after a flight time of τ fl = (100 − t 0 ) Myr.
Running ten simulations to randomise the initial conditions, we counted 78 ± 8 objects within 2 kpc. Such an estimate included all objects with a τ fl ≤ 100 Myr. However, we know that a large fraction of them might not be detectable for several reasons, one of which is the long-term evolution of their absolute magnitudes. If we considered that LP 40−365 stars evolve as proposed by Zhang et al. (2019), i.e. they are 10-1000 times fainter than LP 40−365 in their first ∼ 10 Myr of life, we would expect that only those older than 20 Myr (i.e. 63 ± 6) are detectable by Gaia with G < 18 (assuming the absolute magnitude of LP 40−365, M G = 6.5, as a lower limit). Hence, we could expect to detect just 1-3 intrinsically fainter, younger LP 40−365 stars in a smaller volume (e.g. d < 1 kpc) down to G = 19.
The estimate of 63 stars might be optimistic, given that so far we have detected just three LP 40−365 stars within 2 kpc, suggesting that the timescales proposed by Zhang et al. (2019) are too short. Following this interpretation, all the observed LP 40−365 stars are still relatively young, having ages of the same order of magnitude or shorter than the τ KH ∼ 32 Myr measured for LP 40−365. In this case, we would expect 22 ± 4 LP 40−365 stars. Although this estimate is still too large in comparison with the observed numbers, we note a better agreement with the predictions by Shen et al. (2018b), if we scale our estimate to the number of D 6 stars within a smaller volume with a 1 kpc radius.
Given their relatively low v ej , most of the LP 40−365 stars could stay bound to the Milky Way, as we find that just ≃20 per cent of the simulated objects have v rf > v esc . The fraction of stars achieving v rf > 1.3 v esc is just 4 per cent of the total. Typical proper motions of LP 40−365 stars would also be relatively small, with an average value of µ = 49 +33 −20 mas, indicating that LP 40−365 and J1825−3757 could be the rarest outliers in the µ distribution.
While we stress that these simulations are affected by several unknowns such as the confirmation of a cooling law and the variety of supernova rates, we note that they are sufficiently straightforward to motivate further spectroscopic follow-up of LP 40−365 star candidates. Finally, we also note that the binary-supernova mechanism implies the existence of ejected donor stars, which may have similar kinematics, but different compositions.

SUMMARY AND CONCLUSIONS
Through detailed spectroscopic, photometric, and kinematic analysis, we have gained further evidence in favour of the idea initially advanced by Vennes et al. (2017) and then reinforced in Paper I and Paper II that the high-velocity star LP 40−365 could be a partly burnt white dwarf accretor that survived to a peculiar class of thermonuclear supernovae occurring in binary systems. Given the striking spectroscopic and kinematic similarity we found among LP 40−365 and three new objects (J1603−6613, J1825−3757, J0905+2510), we propose this star as the namesake of its own spectral class that is defined both by spectroscopic, physical, and kinematic properties: (i) Ne-dominated atmospheres, with O and Mg as the second-and third-most-abundant elements, respectively.
(iii) Low mass, below the canonical white dwarf mass.
Because LP 40−365 stars have small radii, which are comparable to those of hot subdwarf stars, and possess atmospheres rich in metals, we also propose a spectral type "sdZ," which continues the tradition of hot subdwarf and white dwarf classification.
The atmospheric composition of LP 40−365 stars is a strong indication of partial C-, O-, and Si-burning, likely connected to thermonuclear explosions that did not entirely disrupt their progenitors. Of significant note is the nearidentical spectroscopic properties of both LP 40−365 and J1603−6613, which strongly suggests a similar formation mechanism. We found some similarities with the simulations of pure deflagrations in CO and ONe white dwarfs, which could be representative of SNe Iax (Fink et al. 2014) and thermonuclear electron-capture supernovae , respectively. We noted that the models do not reproduce the main atmospheric components (O, Ne, Mg), perhaps because the bulk composition of bound remnants is not represented by the photosphere of these stars. Promising simulations by Zhang et al. (2019) describe, for the first time, the evolution of post-SN Iax stars. However, the temporal evolution does not well match those of the LP 40−365 stars. Either different models are required, or the progenitors of LP 40−365 stars may be neither CO nor ONe white dwarfs. Suitable progenitors, still to be tested, could be hybrid CONe white dwarfs.
The physical properties of LP 40−365 stars are precisely constrained by our spectral analysis and Gaia parallaxes, indicating a relatively small range of masses (median masses between 0.20-0.28 M ⊙ ), which indicates they currently are (not fully degenerate) low-mass objects. The radii of these stars were likely inflated due to the large energy accumulated during and after the supernova explosion, growing to 0.16-0.60 R ⊙ , roughly an order of magnitude larger than typical white dwarf radii. Such a large range is compatible with the stars being caught at different evolutionary stages. While the Zhang et al. (2019) models predicted the LP 40−365 stars could be older than ∼ 23 Myr, we speculated that LP 40−365 stars could be still relatively unevolved and thus younger than the Kelvin-Helmholtz timescale of LP 40−365 (∼ 32 Myr), due to the large abundance of heavy metals in their atmospheres.
The kinematics of LP 40−365 stars are also distinctive, indicative of their formation mechanism. We found that they can be either bound or unbound to the Milky Way, depending on whether the stars were ejected in the direction of Galactic rotation (as with LP 40−365 and J1603−6613) or in the opposite direction (like J1825−3757). In order to obtain the observed kinematics, confined to small Z coordinates, the progenitors of LP 40−365 stars likely resided in the Galactic disc, further constraining their ages to no more than a few 10 Myr. The ejection velocities required to reproduce the rest-frame velocity of the three LP 40−365 stars are in the range of 550-600 km s −1 , consistent with ejection from a relatively compact binary systems containing a He-burning donor star of roughly 0.8 M ⊙ orbiting with a 1 hr period. This result is compatible with one class of thermonuclear supernova progenitors, which are invoked to explain short-delay-time supernovae and subluminous classes like SNe Iax (Wang & Han 2009;Wang et al. 2009aWang et al. ,b, 2013. Furthermore, we noted that the kinematic evidence characterising LP 40−365 stars clearly distinguishes them from other groups of hyper-runaways, i.e. the D 6 stars (Shen et al. 2018b) and the hot-subdwarf US 708 (Geier et al. 2015), which require even more compact progenitors to reach their rest frame velocities of 1000 km s −1 .
We estimate that of the order 22 stars could be detectable by the end of the Gaia mission, imposing the boundary conditions implied by the ejection velocities and flight times determined from our analysis. Although the results of our simulations were susceptible to large variations, e.g. the uncertainty of supernova rates and evolutionary timescales of LP 40−365 stars, we stress they are compatible with the observed numbers and scale as 10 per cent of of predicted D 6 stars (Shen et al. 2018b).
The key message from the simulation of the population of LP 40−365 stars within 2 kpc of the Sun is that most of them are likely gravitationally bound to the Milky Way and will possess relatively small proper motions. The confirmed LP 40−365 stars may not represent the broader population yet to be discovered, suggesting that evolutionary "missing links" as well as their ejected donor stars are yet to be found. It is possible that more evolved LP 40−365 stars would show different atmospheric compositions if their heavy elements sink deeper into the atmosphere as the stars evolve with changing temperatures and radii. We anticipate interesting results from more advanced theoretical modelling, expanding the results presented by Fink et al. (2014), Kromer et al. (2015), and Jones et al. (2018) for the composition of bound remnants and Zhang et al. (2019) for the late cooling of SN Iax bound remnants.
We remark that the LP 40−365 stars that remain bound to the Milky Way will revert back to the white dwarf cooling sequence, still displaying lower log g than canonical white dwarfs due to their smaller masses. These white dwarfs should have stratified atmospheres, containing mostly O, Ne, and Mg. Although not fully investigated, peculiar Orich white dwarfs like J1240+6710 (Kepler et al. 2016) -as well as one of our SDSS candidates, J1637+3631 -could represent or be linked to one of the final evolutionary stages of our newly discovered class of LP 40−365 stars.
Programme (FP/2007(FP/ -2013 / ERC Grant Agreement no. 320964 (WDTracer). OFT and BTG were supported by the UK STFC grant ST/P000495. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Based on observations made with the NASA/ESA HST, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program 15431. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 093.D-0431, 097.D−1029,