arXiv is now an independent nonprofit! Learn more
License: CC BY 4.0
arXiv:2607.17488v1 [astro-ph.SR] 20 Jul 2026

Surrogate models for type II supernovae: Probing low-energy explosions and interaction-free regimes

Zhengyang Zhang zhangzhengyang@ynao.ac.cn International Centre of Supernovae (ICESUN), Yunnan Key Laboratory of Supernova Research, Yunnan Observatories, Chinese Academy of Sciences (CAS), Kunming 650216, People’s Republic of China University of the Chinese Academy of Sciences, 19A Yuquan Road, Shijingshan District, Beijing 100049, People’s Republic of China    Shuai Zha International Centre of Supernovae (ICESUN), Yunnan Key Laboratory of Supernova Research, Yunnan Observatories, Chinese Academy of Sciences (CAS), Kunming 650216, People’s Republic of China    Nikhil Sarin Kavli Institute for Cosmology, University of Cambridge, Madingley Road, CB3 0HA, United Kingdom Institute of Astronomy, University of Cambridge, Madingley Road, CB3 0HA, United Kingdom    Takashi J. Moriya National Astronomical Observatory of Japan, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan Graduate Institute for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan School of Physics and Astronomy, Monash University, Clayton, Victoria 3800, Australia    Chengyuan Wu International Centre of Supernovae (ICESUN), Yunnan Key Laboratory of Supernova Research, Yunnan Observatories, Chinese Academy of Sciences (CAS), Kunming 650216, People’s Republic of China    Bo Wang wangbo@ynao.ac.cn International Centre of Supernovae (ICESUN), Yunnan Key Laboratory of Supernova Research, Yunnan Observatories, Chinese Academy of Sciences (CAS), Kunming 650216, People’s Republic of China
(25 June 2026)
Abstract

To address the computational bottleneck of analyzing type II supernova samples from surveys like the Legacy Survey of Space and Time, we present two stella-based neural network surrogates: the interaction model for low-energy explosions with potential circumstellar material (CSM) interaction and the photospheric model for standard interaction-free SNe IIP. Each surrogate follows a two-stage design in which an autoencoder first compresses stella spectral energy distributions (SEDs) into a latent representation and a separate neural-network emulator then maps physical parameters to that latent space. Both models incorporate latent mixup regularization to improve latent-space continuity, while using different network backbones tailored to their respective regimes: ResNet blocks for the interaction model and 2D CNNs for the photospheric model. On test sets, the interaction model and photospheric model achieve normalized SED reconstruction MSEs of ≈9.1×10−5\approx 9.1\times 10^{-5} and 1.0×10−41.0\times 10^{-4}, respectively. For the low-luminosity SN 2005cs, the interaction model favors a low-mass progenitor (MZAMS=10.40−0.05+0.04​M⊙M_{\mathrm{ZAMS}}=10.40^{+0.04}_{-0.05}\,M_{\odot}) with evidence for confined dense CSM, highlighting the likely presence of CSM interaction and offering a physical scenario consistent with direct imaging that helps resolve the historical mass discrepancy. For SN 2012aw, serving as a validation benchmark, the interaction model demonstrates consistency with previous studies by recovering a progenitor mass (MZAMS=11.05−0.06+0.06​M⊙M_{\mathrm{ZAMS}}=11.05^{+0.06}_{-0.06}\,M_{\odot}). For the archetypal SN 1999em, the photospheric model derives a progenitor mass (MZAMS=10.05−0.04+0.07​M⊙M_{\mathrm{ZAMS}}=10.05^{+0.07}_{-0.04}\,M_{\odot}) broadly consistent with direct preexplosion imaging limits without explicit CSM modeling, demonstrating that the photospheric model captures the essential physics of standard type IIP explosions. This model-level implementation can reduce full Bayesian parameter-inference runtime from days to minutes, providing a practical foundation for near-real-time physical characterization in large survey streams.

I Introduction

Time-domain astronomy is entering an era of transient abundance. Surveys such as the Zwicky Transient Facility (ZTF) and the Legacy Survey of Space and Time (LSST) will deliver orders-of-magnitude larger supernova samples over decade timescales, making fast, reproducible, and scalable physics-based modeling and parameter inference a central requirement [1, 2, 3]. Classical event-by-event Bayesian analyses using MCMC or nested sampling can become computationally prohibitive at survey scale, motivating alternative inference strategies, including simulation-based inference with normalizing flows and end-to-end software frameworks for electromagnetic transient inference [4, 5].

Among core-collapse supernovae, hydrogen-rich type II supernovae (SNe II) are among the most common, yet they exhibit substantial diversity in light-curve morphology and spectral evolution [6, 7, 8, 9, 10, 11, 12, 13, 14]. Observational studies have highlighted continuous variations across the classical SN IIP/SN IIL phenomenology in plateau shape, peak luminosity, and velocity evolution, while combined photometric and spectroscopic modeling emphasizes that similar broadband light curves can arise from different physical configurations, exacerbating parameter degeneracies when only limited-band photometry is available [15, 16, 17, 18, 19].

Interaction between SN ejecta and circumstellar material (CSM) strongly shapes early time observables and probes pre-SN mass loss. The underlying radiation-hydrodynamic framework—forward and reverse shocks converting kinetic energy into radiation—has long been established [20, 21, 22]. Although the most obvious signatures appear in narrow-line events (e.g., SNe IIn/Ibn), interaction likely spans a broader SN II population and may be linked to late-stage outbursts, winds, and binary evolution [23, 24, 25, 26, 27, 28]. In particular, confined dense CSM can delay shock breakout and reshape the earliest light-curve rise [29, 30, 31, 32, 33].

Observational and modeling studies now suggest that such confined dense CSM is common among SNe II and can systematically affect the first weeks after explosion [34, 35, 36, 37]. However, broadband light curves alone remain highly degenerate in CSM density structure and radial extent. Adding spectral information can break key degeneracies and better constrain CSM geometry and density profiles, while recent analytic scalings linking peak observables to physical parameters provide additional leverage on the dominant heating mechanism [38, 39].

High-fidelity radiation-hydrodynamics codes such as stella [40, 41, 42] can produce time-evolving SEDs, multiband light curves, and photospheric properties across a broad progenitor and explosion parameter space (e.g., ZAMS mass, explosion energy, 56Ni mass, mass-loss rate, and CSM extent/structure). However, they remain computationally expensive for large-sample inference. Our work builds on the stella-based SN II model grids of Moriya et al. [43], which provide the physical training foundation for fast surrogate modeling in the survey era. This grid has already supported survey-scale inference for ZTF and joint ZTF+ATLAS samples [44, 45, 46]. Related modeling efforts have also been applied to nearby well-observed events such as SN 2023ixf (with complementary radiative-transfer simulations) and to high-redshift SN II samples discovered by JWST/JADES [47, 48, 49, 50, 51, 52, 53].

While precomputed grids dramatically lower the barrier to physics-based modeling, accurate interpolation in continuous parameter space—especially for full time×\timeswavelength SEDs—remains a key bottleneck. Machine-learning surrogate models (emulators) have therefore emerged as a powerful route to accelerate expensive numerical predictions to millisecond-scale evaluations, including recent SN II spectral-emulator work based on tardis [54]. Sarin et al. [55] introduced surrogate models for SN II light curves and photospheric properties trained on a large stella grid and discussed likelihood strategies that explicitly account for surrogate uncertainty in Bayesian fitting. However, extending accurate inference to the low-energy and low-mass regime remains a crucial challenge. Low-luminosity SNe II make up a substantial fraction of the volumetric rate and provide a unique probe of the lower-mass limit of core-collapse progenitors, but they lie near parameter-space boundaries that are often sparsely sampled and numerically unstable [56].

In this work, we focus on the robust generation and inference of SN II light curves and present two independent surrogate models: the interaction model, which is designed to provide stable coverage of low-energy and low-luminosity regimes with possible CSM interaction, and the photospheric model, which provides dedicated support for models without circumstellar interaction. Our central objective is to construct a surrogate capable of accurately generating time-evolving spectra and multiband light curves; these tools have been integrated into redback111https://github.com/nikhil-sarin/redback, an open-source Bayesian inference software package for electromagnetic transients [5]. Both models adopt a two-stage architecture. For the interaction model, we first train an autoencoder to compress SEDs into a 256-dimensional latent space, followed by a parameter-to-latent emulator. For the photospheric model, we employ a 2D convolutional autoencoder to compress high-dimensional SEDs into a low-dimensional latent representation, followed by a similar parameter-to-latent emulation step. To enhance continuity and robustness, both architectures introduce latent-space interpolation-consistency regularization (latent mixup). Furthermore, we incorporate time-weighted losses and gradient constraints to improve reconstruction fidelity, particularly during the critical early time and transition phases.

The paper is organized as follows. Section II details the methodology, including the treatment of the training dataset, the specific architectures of the interaction model and photospheric model surrogates, and the implementation of the latent mixup regularization. Section III presents the validation results, quantifying the spectral reconstruction fidelity and demonstrating the Bayesian parameter inference on three benchmark supernovae (SN 2005cs, SN 2012aw, and SN 1999em). Finally, Sec. IV summarizes our findings, discusses the computational advantages, and addresses current physical limitations.

II Method

II.1 Model overview

Figure 1 summarizes the network architecture used in this work. Both surrogate models adopt a two-stage autoencoder–emulator design: the autoencoder first compresses each time–wavelength spectral energy distribution (SED) into a compact latent representation and learns a nonlinear decoder back to spectral space; meanwhile, an independent emulator maps the physical input parameters to this latent space so that the frozen decoder can generate a complete SED during inference. The interaction model is designed for regimes affected by CSM interaction, including low-mass and low-energy explosions, and uses ResNet-based encoder, decoder, and emulator networks to model the interacting SED manifold. The photospheric model targets the classical interaction-free type IIP regime. It uses a two-dimensional convolutional neural network (CNN) as its backbone. In both cases, the autoencoder input is a two-dimensional synthetic SED image. Before training, the stella spectra are interpolated onto a fixed 100×100100\times 100 time–wavelength grid, so each input sample contains 10,000 normalized values. The two models have the same input dimensionality, but their grids are adjusted according to the relevant physical mechanisms. For the interaction model, the time axis spans 0.1–400 days and the wavelength axis is geometrically spaced from 500 to 49,500 Å. For the photospheric model, the time axis is sampled from 0.1 to 200 days and the preprocessing uses a linearly spaced wavelength grid over the same 500–49,500 Å range. The physical parameter ranges and valid model subsets differ, as summarized in Table 1.

Our implementation differs from the setup of Sarin et al. [55] in three main respects. First, we use the autoencoder latent representation directly instead of adding a PCA compression step; this avoids imposing a linear projection on an intrinsically nonlinear SED manifold and allows the latent coordinates to be optimized for spectral reconstruction without an additional PCA-to-decoder mapping. Second, we regularize the learned manifold with latent-space mixup, motivated by the broader manifold-mixup approach for enforcing smoother interpolations in learned representations [57]; the ablation test for this choice is presented in the Appendix, consistent with previous work showing that interpolation regularization can improve autoencoder latent spaces [58]. Third, we extend the surrogate coverage toward lower-mass and lower-energy explosions and add the photospheric model as a dedicated interaction-free surrogate for standard type IIP light curves.

The interaction model adopts a confined CSM density structure characterized by a wind acceleration law [43]. The CSM density ρCSM​(r)\rho_{\mathrm{CSM}}(r) is defined as

ρCSM​(𝐫)=𝐌˙𝟒​π​𝐫𝟐​𝐯wind​(𝐫),\mathbf{\rho_{\mathrm{CSM}}(r)=\frac{\dot{M}}{4\pi r^{2}v_{\mathrm{wind}}(r)},} (1)

where M˙\dot{M} is the mass-loss rate. The velocity profile vwind​(r)v_{\mathrm{wind}}(r) follows a β\beta-law acceleration

𝐯wind​(𝐫)=𝐯𝟎+(𝐯∞−𝐯𝟎)​(𝟏−𝐑𝟎𝐫)β.\mathbf{v_{\mathrm{wind}}(r)=v_{0}+(v_{\infty}-v_{0})\left(1-\frac{R_{0}}{r}\right)^{\beta}.} (2)

Here, R0R_{0} is the progenitor radius, v0v_{0} is the initial velocity at the surface, and v∞v_{\infty} is the terminal wind velocity (fixed at 10​km​s−110\,\mathrm{km\,s}^{-1}). The parameter β\beta determines the efficiency of wind acceleration; a higher β\beta corresponds to a slower acceleration, resulting in a steeper density gradient in the immediate vicinity of the progenitor.

Table 1: Input physical parameters and their valid ranges for the interaction model and photospheric model surrogates. The last column summarizes how the interaction model differs from Sarin et al. [55, Table 1]. Note that the photospheric model includes derived parameters (MenvM_{\mathrm{env}}, R0R_{0}) dependent on MZAMSM_{\mathrm{ZAMS}} (marked with superscript aa), while the interaction model includes specific CSM interaction parameters. Mixing definitions: “cm” (central mixing) denotes no 56Ni mixing into the hydrogen-rich envelope; “fm” (full mixing) indicates a fully mixed 56Ni distribution throughout the entire ejecta; “hm” (half mixing) represents an intermediate extent. The interaction model inherently adopts the “hm” configuration.
Parameter Symbol Unit Interaction model range Photospheric model range Interaction model vs. Sarin+25
Progenitor & Explosion Parameters
Progenitor Mass MZAMSM_{\mathrm{ZAMS}} M⊙M_{\odot} 9−189-18 10−1810-18 Lower bound extended (10 →\rightarrow 9)
56Ni Mass MNiM_{\mathrm{Ni}} M⊙M_{\odot} 0.001−0.30.001-0.3 0.0001−0.30.0001-0.3 Unchanged
Explosion Energy ESNE_{\mathrm{SN}} 105110^{51} erg 0.1−5.00.1-5.0 0.2−5.00.2-5.0 Lower bound extended (0.5 →\rightarrow 0.1)
Mixing Parameter mixing – Fixed (hm) {cm, fm, hm} Explicitly fixed to hm
CSM configuration (interaction model only)
Mass-loss Rate log10⁡(M˙)\log_{10}(\dot{M}) log⁡(M⊙​yr−1)\log(M_{\odot}\,\mathrm{yr}^{-1}) −5∼−1-5\sim-1 – Unchanged
CSM Radius RCSMR_{\mathrm{CSM}} 101410^{14} cm 1−101-10 – Unchanged
Density Slope β\beta – 0.5−5.00.5-5.0 – Unchanged
Derived Envelope Parameters (photospheric model only)
Envelope Mass MenvM_{\mathrm{env}} M⊙M_{\odot} – 7.175−9.57.175-9.5111These are derived photospheric-model parameters that depend on MZAMSM_{\mathrm{ZAMS}}. N/A
Progenitor Radius R0R_{0} R⊙R_{\odot} – 414−970414-970111These are derived photospheric-model parameters that depend on MZAMSM_{\mathrm{ZAMS}}. N/A

II.2 Interaction model

II.2.1 Data preprocessing and augmentation

The interaction model dataset is based on the stella radiation-hydrodynamics simulation grid of Moriya et al. [43], supplemented by low-energy and low-mass models obtained from Moriya [59]. After filtering invalid or incomplete outputs, the Interaction dataset contains 298,375 original valid SED samples on a uniform 100×100100\times 100 time-wavelength grid. Each sample is defined by six physical parameters, including MZAMSM_{\mathrm{ZAMS}}, MNiM_{\mathrm{Ni}}, ESNE_{\mathrm{SN}}, log10⁡M˙\log_{10}\dot{M}, RCSMR_{\mathrm{CSM}}, and β\beta, together with the corresponding time-dependent stella SED. The spectral data were interpolated onto the uniform 100×100100\times 100 grid and log-normalized to the [0,1][0,1] interval.

To enhance model robustness, we applied dimensionless Gaussian perturbations with σ=0.01\sigma=0.01 in the normalized training space. Specifically, noise was added both to the standardized input-parameter vector XX and to the output array YY (the log10⁡Lν\log_{10}L_{\nu} SED values normalized to the [0,1][0,1] interval), generating five augmented copies for each original sample. Including the original samples and their augmented copies, this gives a total of 1,790,250 samples.

The models were split into training, validation, and test groups with an 80%/10%/10% ratio. Each original model and all of its Gaussian-perturbed augmented copies were assigned exclusively to the same subset. This prevents augmented versions of the same physical model from crossing split boundaries and therefore avoids data leakage between the training, validation, and test sets.

Refer to caption
Figure 1: Schematic illustration of the two-stage surrogate modeling framework. This schematic is AI-generated for visualization purposes only. The two panels display the workflows for the spectral compression (Stage 1; blue background) and physical parameter emulation (Stage 2; green background) phases. The left panel shows the AutoEncoder architecture, which compresses high-dimensional true spectra into a low-dimensional latent vector (𝐳\mathbf{z}) by minimizing the reconstruction loss. The right panel shows the emulator training process, where the network maps physical parameters to a predicted latent vector (𝐳pred\mathbf{z}_{\text{pred}}). Note that the pretrained decoder is frozen (indicated by the padlock icon) during this stage to ensure consistent spectral generation from the predicted latent codes.

II.2.2 Spectral compression via latent mixup

The primary objective of the first stage, as illustrated in the blue left panel of Figure 1, is to compress the high-dimensional spectral input 𝐱∈ℝ10000\mathbf{x}\in\mathbb{R}^{10000} into a low-dimensional latent vector 𝐳∈ℝ256\mathbf{z}\in\mathbb{R}^{256}. Both the encoder and decoder are constructed with 12 residual blocks, utilizing SiLU activation functions to capture nonlinear spectral features.

Standard autoencoders often suffer from discontinuous or fragmented latent spaces [58]. This irregularity makes it difficult for the model to accurately map continuous physical parameters to latent codes during the second stage. To resolve this mapping challenge, we introduce a regularization technique termed “latent mixup consistency” [57]. Its core concept is to ensure that linear interpolation in the latent space directly corresponds to a linear transition in the input spectral space. The total loss function for this stage is defined as

ℒtotal=ℒrecon​(𝐱,𝐱^)+wmixup⋅ℒmixup,\mathcal{L}_{\text{total}}=\mathcal{L}_{\text{recon}}(\mathbf{x},\hat{\mathbf{x}})+w_{\text{mixup}}\cdot\mathcal{L}_{\text{mixup}}, (3)

where ℒrecon\mathcal{L}_{\text{recon}} enforces spectral fidelity to the original stella outputs, while wmixupw_{\text{mixup}} sets how strongly the interpolation regularization is applied. In practice, this balances two goals: accurate spectrum reconstruction and a smooth latent manifold that is easier to map from physical parameters. The Mixup consistency term is defined as

ℒmixup=‖𝒟​(λ​𝐳1+(1−λ)​𝐳2)−(λ​𝐱1+(1−λ)​𝐱2)‖2,\mathcal{L}_{\text{mixup}}=\left\|\mathcal{D}(\lambda\mathbf{z}_{1}+(1-\lambda)\mathbf{z}_{2})-(\lambda\mathbf{x}_{1}+(1-\lambda)\mathbf{x}_{2})\right\|^{2}, (4)

Here, 𝒟\mathcal{D} denotes the decoder; 𝐳1\mathbf{z}_{1} and 𝐳2\mathbf{z}_{2} are the latent codes of two samples; 𝐱1\mathbf{x}_{1} and 𝐱2\mathbf{x}_{2} are the corresponding input spectra; and λ∼Beta​(α,α)\lambda\sim\text{Beta}(\alpha,\alpha) is the random mixing coefficient [60]. The parameter α\alpha controls the shape of the Beta distribution from which λ\lambda is drawn: smaller values place the mixed sample closer to one of the two original samples, whereas larger values concentrate λ\lambda closer to their average. In this work, we set α=0.2\alpha=0.2. Physically, this term encourages intermediate points in latent space to decode into intermediate spectra, rather than producing nonphysical artifacts between neighboring models. This improves interpolation across sparsely sampled regions of parameter space and stabilizes the downstream inference of progenitor and explosion properties.

II.2.3 Physical parameter emulation

The second stage, shown in the green right panel of Figure 1, learns the mapping from the six-dimensional physical parameter space 𝐩∈ℝ6\mathbf{p}\in\mathbb{R}^{6} to the latent code 𝐳\mathbf{z}. For this purpose, we construct an emulator network consisting of 12 ResNet blocks with SiLU activations between layers.

During this phase, the pretrained decoder from Stage 1 remains frozen to preserve the learned manifold (as indicated by the padlock icon in Figure 1). The emulator predicts the latent code 𝐳^\hat{\mathbf{z}}, which is then passed through the frozen decoder to generate the reconstructed spectrum. The training objective minimizes a composite loss function consisting of a Huber loss in latent space and an L1 loss in spectral space.

II.3 Photospheric model

The foundational dataset comprises 1,032 stella models (Moriya [59]), each defined by six physical parameters (MZAMSM_{\mathrm{ZAMS}}, MNiM_{\mathrm{Ni}}, mixing parameter, EexpE_{\mathrm{exp}}, MenvM_{\mathrm{env}}, and R0R_{0}). To achieve robust training despite this relatively limited sample size, we implemented a 10×\times data augmentation strategy by injecting dimensionless Gaussian noise with σ=0.01\sigma=0.01 in the normalized training space, applied both to the standardized input-parameter vector XX and to the [0,1][0,1]-normalized log10⁡Lν\log_{10}L_{\nu} SED output array YY, expanding the total dataset to 10,320 samples. The photospheric model dataset was split into training, validation, and test sets using an 80%/10%/10% ratio, corresponding to 8,256, 1,032, and 1,032 augmented samples, respectively.

Regarding the network backbone, we observed that standard ResNet architectures yielded suboptimal performance on this specific dataset. Under the same setting (ResNet vs. CNN), ResNet shows a 42.15% higher test-set MSE than CNN. Consequently, we adopted a 2D CNN. This architecture retains the established two-stage training strategy (as illustrated in Figure 1) and incorporates latent mixup regularization to ensure the smoothness and continuity of the latent manifold.

III Result

III.1 Model performance

To quantitatively assess the generative accuracy, we evaluated both stages of the surrogate pipeline on held-out test sets: the Stage 1 autoencoder, which measures the intrinsic SED compression–reconstruction accuracy, and the Stage 2 emulator, which measures the full parameter-to-SED surrogate accuracy after mapping physical parameters into the learned latent space. Errors are reported in the normalized training space (YY, the [0,1][0,1]-normalized log10⁡Lν\log_{10}L_{\nu} array) and after conversion back to logarithmic luminosity units, expressed in dex.

For the interaction model, the Stage 1 autoencoder achieves MSEY=8.857866×10−5\mathrm{MSE}_{Y}=8.857866\times 10^{-5}, and the complete Stage 2 emulator gives a very similar MSEY=9.134600×10−5\mathrm{MSE}_{Y}=9.134600\times 10^{-5}, indicating that most of the error budget is already set by the autoencoder reconstruction limit rather than by the parameter-to-latent mapping. For the photospheric model CNN, the Stage 1 and Stage 2 test errors are likewise nearly identical, with MSEY=1.001196×10−4\mathrm{MSE}_{Y}=1.001196\times 10^{-4} and 1.005627×10−41.005627\times 10^{-4}, respectively.

Refer to caption
Figure 2: Spectral reconstruction performance on the held-out test sets. The top row (a–c) shows the photospheric model, while the bottom row (d–f) shows the interaction model. From left to right, the columns show the ground-truth stella SEDs, the emulator-predicted SEDs, and the residual maps. In the true and predicted SED panels, the color scale represents the dimensionless epoch-normalized relative luminosity, Lλ​(t,λ)/Lλ,max​(t)L_{\lambda}(t,\lambda)/L_{\lambda,\max}(t), where Lλ,max​(t)L_{\lambda,\max}(t) is the maximum LλL_{\lambda} at the same epoch. Redder regions therefore correspond to wavelengths close to the normalized SED peak at that epoch. The residual panels show the fractional residual, (Lλ,pred−Lλ,true)/Lλ,true(L_{\lambda,\rm pred}-L_{\lambda,\rm true})/L_{\lambda,\rm true}, computed from the per-epoch normalized relative LλL_{\lambda} maps; low-signal regions with true normalized luminosity below 0.05 are shown in gray because fractional residuals in weak spectral tails are not informative. For visualization, the top-row photospheric model panels use a linear wavelength axis, matching the linear wavelength grid used for its training and preprocessing, whereas the bottom-row interaction model panels use a logarithmic (geometric) wavelength axis, matching the geometric wavelength grid used during training.

Figure 2 illustrates the reconstruction fidelity of our emulators across the full time–wavelength plane using the epoch-normalized quantity Lλ/Lλ,maxL_{\lambda}/L_{\lambda,\max}. This normalization is used only for visualization: it highlights the temporal evolution of the spectral shape. The progressive shift of the brightest normalized regions toward longer wavelengths is expected as type II SN ejecta cool and the SED peak moves redward. The true and predicted two-dimensional SED maps agree closely for both the photospheric model (panels a–c) and the interaction model (panels d–f). The residual panels show the fractional residual (Lλ,pred−Lλ,true)/Lλ,true(L_{\lambda,\rm pred}-L_{\lambda,\rm true})/L_{\lambda,\rm true} in the same normalized space, with weak-tail pixels masked in gray where the true normalized luminosity is below 0.05.

This high fidelity at the spectral level translates directly into the accurate reconstruction of photometric observables, as illustrated in Figure 3.

Refer to caption
Figure 3: Comparison of multiband absolute magnitude light curves among ground-truth stella simulations (solid lines), our emulator predictions (dashed lines), and the Sarin25 model (dotted lines). The colors correspond to the ZTF filters: green for the gg-band, red for the rr-band, and purple for the ii-band. Top row (interaction model): Validation in the low-mass, low-energy, CSM-interacting regime. Bottom row (photospheric model): Validation in the interaction-free regime.

For the top-row interaction model cases, which sample the low-mass, low-energy, CSM-interacting regime, the predictions (dashed lines) closely track the ground-truth stella simulations (solid lines). As demonstrated in panels (a)–(c), the model maintains high fidelity across different parameter combinations, including both the rapid early rise and the transition phase.

For the photospheric model (bottom row), the emulator demonstrates a robust capability to capture interaction-free spectral features despite the limited training data, accurately reproducing the early cooling peak, the flat plateau phase, and the sharp transition to the radioactive decay tail.

III.2 Parameter inference validation

Before presenting the benchmark applications, we describe the inference setup used throughout this section. The benchmark fits are performed in broadband photometric space using the redback/bilby inference interface with the dynesty nested sampler. For a given parameter vector 𝜽\boldsymbol{\theta}, the surrogate first generates a time-dependent stella SED. The SED is then convolved with the corresponding filter response functions to produce synthetic model magnitudes in the observed bands. These synthetic magnitudes are compared with the observed multiband photometry in the likelihood evaluation.

Both surrogate models are implemented through the redback_surrogates222https://github.com/nikhil-sarin/redback_surrogates interface to redback, an open-source Bayesian inference package for electromagnetic transients [5]. The forward models support optional PyTorch GPU evaluation for generating synthetic light curves, but this acceleration acts at the forward-model level rather than as a unified native backend for the full redback sampling stack; consequently, the nested-sampling workflow used here is still primarily CPU executed.

The priors adopted in the benchmark inference are summarized in Table 2. Most physical parameters are assigned uniform priors over the range covered by the corresponding surrogate grid. The nuisance parameter AA is included to absorb residual uncertainty in distance, extinction, calibration, and model normalization.

Table 2: Priors used in the benchmark inference runs. ‘U” denotes a uniform prior, ‘logU” denotes a log-uniform prior, and “Delta” denotes a fixed value.
Parameter interaction model prior photospheric model prior Unit / note
MZAMSM_{\rm ZAMS} U(9,18)(9,18) U(10,18)(10,18) M⊙M_{\odot}
ESNE_{\rm SN} U(0.1,5)(0.1,5) U(0.2,5)(0.2,5) 105110^{51} erg
MNiM_{\rm Ni} U(0.001,0.3)(0.001,0.3) U(0.0001,0.3)(0.0001,0.3) M⊙M_{\odot}
log10⁡M˙\log_{10}\dot{M} U(−5,−1)(-5,-1) Not used M⊙​yr−1M_{\odot}\,{\rm yr}^{-1}
β\beta U(0.5,5)(0.5,5) Not used CSM density-slope parameter
RCSMR_{\rm CSM} U(1,10)(1,10) Not used 101410^{14} cm
mixing Not used Categorical {0,1,2}\{0,1,2\} 0=cm0={\rm cm}, 1=fm1={\rm fm}, 2=hm2={\rm hm}
MenvM_{\rm env} Not used U(7.2,9.5)(7.2,9.5) M⊙M_{\odot}
R0R_{0} Not used U(510,970)(510,970) R⊙R_{\odot}
AA logU(0.01,2)(0.01,2) logU(0.01,2)(0.01,2) Multiplicative flux scale

We adopt a Gaussian likelihood in magnitude space. For each observed photometric data point ii, the effective uncertainty is defined as

σi,eff2=σi,obs2+σadd2,\sigma_{i,\rm eff}^{2}=\sigma_{i,\rm obs}^{2}+\sigma_{\rm add}^{2}, (5)

where σi,obs\sigma_{i,\rm obs} is the reported photometric uncertainty and σadd\sigma_{\rm add} is a fixed additional uncertainty term added in quadrature to account for residual calibration, interpolation, and surrogate-model uncertainties. The log-likelihood is then

ln⁡ℒ​(𝜽)=−12​∑i=1Nphot[(mi,obs−mi,mod​(𝜽))2σi,eff2+ln⁡(2​π​σi,eff2)],\ln\mathcal{L}(\boldsymbol{\theta})=-\frac{1}{2}\sum_{i=1}^{N_{\rm phot}}\left[\frac{\left(m_{i,\rm obs}-m_{i,\rm mod}(\boldsymbol{\theta})\right)^{2}}{\sigma_{i,\rm eff}^{2}}+\ln\left(2\pi\sigma_{i,\rm eff}^{2}\right)\right], (6)

where mi,obsm_{i,\rm obs} is the observed magnitude and mi,mod​(𝜽)m_{i,\rm mod}(\boldsymbol{\theta}) is the corresponding model magnitude generated from the surrogate SED after filter convolution. This likelihood definition is used consistently for all benchmark supernovae analyzed below.

We selected three benchmark supernovae and performed inference in broadband photometric space. Although our surrogates generate time-dependent stella SEDs, these SEDs are multigroup radiation-transfer outputs sampled with a finite number of frequency bins, and are therefore more appropriate for constructing synthetic broadband light curves than for direct line-profile-level comparison with observed spectra. For the interaction model, we analyzed SN 2005cs, which tests the robustness of the model in the presence of potential confined CSM interaction, as suggested by broad spectral “ledge” features. For the photospheric model, we employed the archetypal SN 1999em to characterize the standard regime dominated by pure ejecta evolution. Additionally, SN 2012aw was reanalyzed to facilitate comparison with previous model iterations.

III.2.1 SN 2005cs

Refer to caption
Figure 4: Parameter inference validation for SN 2005cs using the interaction model. Panel (a): Comparison between the observed multiband photometry and the best-fit synthetic light curves (solid lines). The error bars represent a fixed 0.2 mag uncertainty around the light curves corresponding to the median inferred parameters. Panel (b): Corner plot displaying the posterior probability distributions for the inferred physical parameters.

SN 2005cs serves as a critical stress test for our interaction model because of the long-standing tension in its inferred physical parameters. The multiband photometric data used in our fit are adopted from the observations presented by Pastorello et al. [61]. While direct analyses of preexplosion HST images consistently support a low-mass progenitor (6−10​M⊙6-10\,M_{\odot}; Maund et al. 62, Li et al. 63, Eldridge et al. 64), early hydrodynamic modeling required a significantly more massive envelope (∼17.3​M⊙\sim 17.3\,M_{\odot}) to reproduce the plateau properties [65]. Our inference results (Fig. 4) offer a plausible resolution to this discrepancy. The interaction model provides a high-fidelity fit to the multiband light curves, tightly constraining the progenitor mass to MZAMS=10.40−0.05+0.04​M⊙M_{\mathrm{ZAMS}}=10.40^{+0.04}_{-0.05}\,M_{\odot} and the explosion energy to ESN=0.1515−0.0025+0.0018×1051E_{\mathrm{SN}}=0.1515^{+0.0018}_{-0.0025}\times 10^{51} erg. These values agree well with both direct imaging constraints and recent neutrino-driven explosion models [66], favoring the low-mass hypothesis. The inferred 56Ni mass, MNi=0.00639−0.00008+0.00009​M⊙M_{\mathrm{Ni}}=0.00639^{+0.00009}_{-0.00008}\,M_{\odot}, is also consistent with expectations for nickel-poor, low-energy iron core-collapse explosions.

Furthermore, the model characterizes the circumstellar environment as a confined dense shell (RCSM=4.237−0.057+0.055×1014R_{\rm CSM}=4.237^{+0.055}_{-0.057}\times 10^{14} cm) featuring an enhanced mass-loss rate (log10⁡M˙=−2.894−0.012+0.012\log_{10}\dot{M}=-2.894^{+0.012}_{-0.012}) and a steep density gradient (β=2.204−0.042+0.034\beta=2.204^{+0.034}_{-0.042}). These properties provide a self-consistent physical framework for the interaction signatures observed in early spectra. This scenario receives critical observational support from the recent discovery of SN 2024abfl, a physical “twin” to SN 2005cs. Unlike the ambiguous early spectra of SN 2005cs—where the origin of similar features remained debated [67]—SN 2024abfl clearly displayed a broad “ledge” feature near 4600 Å within days of explosion. This feature is interpreted as a blend of shock-accelerated high-ionization lines (e.g., He II and N III) arising from interaction with confined CSM [68, 69]. This indicates that even low-mass red supergiant progenitors (∼9\sim 9–12​M⊙12\,M_{\odot}) can experience enhanced wind mass loss. Consequently, SN 2024abfl provides compelling independent verification that the peculiar early spectral features of SN 2005cs were driven by CSM interaction rather than intrinsic ejecta properties.

III.2.2 SN 2012aw

Refer to caption
Figure 5: Same as Fig. 4, but for SN 2012aw.

As illustrated in Figure 5, we analyze SN 2012aw using the optical photometric dataset presented by Bose et al. [70]. The inferred progenitor ZAMS mass (MZAMS=11.05−0.06+0.06​M⊙M_{\mathrm{ZAMS}}=11.05^{+0.06}_{-0.06}\,M_{\odot}) agrees well with both direct progenitor-imaging constraints (≈12.5​M⊙\approx 12.5\,M_{\odot}; Fraser et al. 71, Kochanek et al. 72, Fraser 73) and the estimate of Sarin et al. [55] (10.61​M⊙10.61\,M_{\odot}). For the 56Ni mass, our model yields a robust estimate: the derived MNi=0.0557−0.0005+0.0006​M⊙M_{\mathrm{Ni}}=0.0557^{+0.0006}_{-0.0005}\,M_{\odot} closely matches the 0.06±0.01​M⊙0.06\pm 0.01\,M_{\odot} measured by Bose et al. [70] and is consistent with the ∼0.06​M⊙\sim 0.06\,M_{\odot} obtained by Dall’Ora et al. [74]. We derive an explosion energy of ESN=0.533−0.005+0.005×1051E_{\mathrm{SN}}=0.533^{+0.005}_{-0.005}\times 10^{51} erg, comparable to the 0.63−0.04+0.05×10510.63_{-0.04}^{+0.05}\times 10^{51} erg reported by Sarin et al. [55]. The high-velocity absorption features identified in early spectra by Bose et al. [70] suggest ejecta-CSM interaction, while radio analysis by Yadav et al. [75] reveals significant inverse Compton cooling, supporting a scenario in which the supernova resided in a dense environment with an intense radiation field. These findings support our inference that SN 2012aw is surrounded by a dense, confined CSM shell.

III.2.3 SN 1999em

Refer to caption
Figure 6: Same as Fig. 5, but for the archetypal type IIP SN 1999em using the interaction-free photospheric model.

SN 1999em is widely regarded as an archetype of normal type IIP supernovae, serving as an ideal benchmark for our interaction-free photospheric model. The photometric dataset used for this benchmark is taken from the optical monitoring presented by Elmhamdi et al. [76]. Our emulator reproduces the observed photometry reasonably without invoking explicit CSM interaction, deriving a progenitor mass of MZAMS=10.05−0.04+0.07​M⊙M_{\mathrm{ZAMS}}=10.05^{+0.07}_{-0.04}M_{\odot}. This estimate is of particular interest regarding the historical mass discrepancy: while previous hydrodynamical studies (e.g., Bersten et al. 77) favored higher masses (∼19​M⊙\sim 19M_{\odot}) to reproduce the light curve, our inference aligns with the direct preexplosion imaging limits (12±1​M⊙12\pm 1M_{\odot}; Smartt et al. 78).

Regarding explosion properties, we derive a synthesized 56Ni mass of MNi=0.033−0.001+0.001​M⊙M_{\mathrm{Ni}}=0.033^{+0.001}_{-0.001}M_{\odot} with negligible mixing. Our derived values are broadly consistent with previous photometric and spectroscopic estimates in the literature [76, 79]. This photospheric model interpretation is also consistent with the literature view of SN 1999em as a normal type IIP event: its photometric and spectroscopic evolution was extensively studied by Elmhamdi et al. [76], and its bolometric light curve and Hα\alpha evolution can be reproduced by conventional hydrodynamic/photospheric modeling [80]. In addition, x-ray and radio observations indicate only weak ejecta–CSM interaction with a low-density red-supergiant wind, corresponding to a mass-loss rate of order 10−610^{-6}–10−6.5​M⊙​yr−110^{-6.5}\,M_{\odot}\,{\rm yr}^{-1} rather than a dense compact CSM component dominating the optical light curve [81, 82]. These results demonstrate that the photospheric model captures the essential physics of standard type IIP explosions. However, we also observe a systematic discrepancy during the first ≲10\lesssim 10 days, where the observed light curve is ≈0.3−0.6\approx 0.3-0.6 mag brighter than our best-fit interaction-free model. This excess may reflect weak CSM interaction or other effects not included in the present photospheric framework, such as cooling-envelope emission, uncertainties in the early color/temperature evolution, line-blanketing effects, or limitations of the simplified model assumptions. Nevertheless, the photospheric model remains a viable alternative for characterizing type IIP events, especially when focusing on the main plateau phase.

IV Discussion and conclusion

While our emulators demonstrate high statistical precision, their physical accuracy is fundamentally bounded by the underlying training data generated by stella. A primary limitation is the assumption of local thermodynamic equilibrium (LTE). This approximation holds well during the optically thick photospheric phase but breaks down as the ejecta expand and transition into the nebular phase. As explicitly noted by Sarin et al. [55], LTE-based predictions become increasingly unreliable at epochs ≳100\gtrsim 100 days postexplosion. Consequently, any inference derived from the late-time radioactive decay tail should be interpreted with caution. Future iterations could integrate non-LTE codes like cmfgen or sumo to address this regime, although the computational cost of generating dense non-LTE grids remains a significant challenge.

The applicability of the interaction model is also bounded by the parameter space of the training grid, specifically the upper mass-loss rate limit of M˙=10−1.0​M⊙​yr−1\dot{M}=10^{-1.0}\,M_{\odot}\,\mathrm{yr}^{-1}. This range is sufficient for events with weak interaction, but it is not suitable for characterizing extreme interaction events such as SN 2013fs [35]. Such objects typically require significantly higher CSM densities (M˙>0.15​M⊙​yr−1\dot{M}>0.15\,M_{\odot}\,\mathrm{yr}^{-1}) to reproduce their strong, persistent narrow emission lines—a regime that remains to be explored in future grid expansions.

In summary, to address the analysis bottlenecks for massive type II supernova datasets, we have presented two independent neural network surrogate models—the interaction model and the photospheric model—designed for distinct physical regimes. The interaction model targets low-energy explosions with possible CSM interaction, whereas the photospheric model provides dedicated support for standard interaction-free events. We tailored specific machine learning architectures to these regimes: the interaction model utilizes a ResNet-based autoencoder for complex interaction scenarios, while the photospheric model employs a 2D CNN optimized for SED reconstruction of standard events. Crucially, both architectures incorporate “latent mixup” regularization to effectively smooth the high-dimensional latent space. We summarize the posterior medians and 68% credible intervals for the three benchmark supernovae in Table 3, which demonstrates that the surrogates recover physically plausible progenitor, explosion, and CSM properties across both interaction and interaction-free regimes.

Table 3: Posterior medians and 68% credible intervals for the physical parameters derived from our surrogate models.
Interaction model
Name MZAMSM_{\rm ZAMS} [M⊙M_{\odot}] MNiM_{\rm Ni} [M⊙M_{\odot}] ESNE_{\rm SN} [foe] log10⁡M˙\log_{10}\dot{M} β\beta RCSMR_{\rm CSM} [101410^{14} cm]
SN 2012aw 11.05−0.06+0.0611.05^{+0.06}_{-0.06} 0.0557−0.0005+0.00060.0557^{+0.0006}_{-0.0005} 0.533−0.005+0.0050.533^{+0.005}_{-0.005} −1.76−0.010+0.010-1.76^{+0.010}_{-0.010} 2.82−0.03+0.022.82^{+0.02}_{-0.03} 1.250−0.004+0.0041.250^{+0.004}_{-0.004}
SN 2005cs 10.40−0.05+0.0410.40^{+0.04}_{-0.05} 0.00639−0.00008+0.000090.00639^{+0.00009}_{-0.00008} 0.1515−0.0025+0.00180.1515^{+0.0018}_{-0.0025} −2.894−0.012+0.012-2.894^{+0.012}_{-0.012} 2.204−0.042+0.0342.204^{+0.034}_{-0.042} 4.237−0.057+0.0554.237^{+0.055}_{-0.057}
photospheric model
MZAMSM_{\rm ZAMS} [M⊙M_{\odot}] MNiM_{\rm Ni} [M⊙M_{\odot}] ESNE_{\rm SN} [foe] MenvM_{\rm env} [M⊙M_{\odot}] R0R_{0} [R⊙R_{\odot}] mixing
SN 1999em 10.05−0.04+0.0710.05^{+0.07}_{-0.04} 0.033−0.001+0.0010.033^{+0.001}_{-0.001} 0.512−0.013+0.0140.512^{+0.014}_{-0.013} 7.66−0.02+0.027.66^{+0.02}_{-0.02} 543−8+7543^{+7}_{-8} 0.00.0

These surrogates are distributed through redback_surrogates. In this implementation, a full Bayesian fit for a single type II supernova can be completed in minutes rather than the days typically required when each likelihood evaluation calls a radiation-hydrodynamics simulation. We are currently applying these surrogates to well-observed type II supernova samples from the CSP and ZTF surveys [83]. This statistical application paves the way for the real-time physical characterization of the tens of thousands of SNe II expected annually from LSST [1], thereby enhancing our capability to derive statistical constraints on massive star explosions.

Acknowledgements.
We thank the anonymous reviewer for the valuable comments and suggestions. We acknowledge the Stellar Physics Group of Yunnan Observatories for providing computational resources. This study was supported by the National Natural Science Foundation of China (Nos. 12288102, 12225304, 12090040/12090043), the CAS Project for Young Scientists in Basic Research (No. YSBR-148), the National Key R&D Program of China (No. 2021YFA1600404), the science research grant from the China Manned Space Project (No. CMS-CSST-2021-A12), the Yunnan Revitalization Talent Support Program (Yunling Scholar Project), the Yunnan Fundamental Research Project (Nos. 202201BC070003, 202501AS070005), the Yunnan Science and Technology Program (Nos. 202605AS350010 and 202601BC070011), and the International Centre of Supernovae (ICESUN), Yunnan Key Laboratory of Supernova Research (Nos. 202302AN360001 and 202505AV340004). S.Z. is supported by the National Natural Science Foundation of China (Nos. 12393811, 12473031), the Yunnan Revitalization Talent Support Program (Young Talent Project), and the Yunnan Fundamental Research Project (No. 202501AS070078). C. Wu is supported by the National Natural Science Foundation of China (No. 12473032), the Yunnan Revitalization Talent Support Program (Young Talent Project), and the Yunnan Fundamental Research Project (Nos. 202501AW070001, 202301AU070039).

Data availability

The software used in this work is publicly available through redback_surrogates at Ref. [84]. Requests for other data or further information should be sent to the authors.

Appendix A Validation of latent space smoothness

For the emulator to accurately map continuous physical parameters to spectra, the underlying latent space (Stage 1) must be smooth and free of discontinuities. A fragmented latent space would lead to prediction artifacts.

Refer to caption
Figure 7: Latent space interpolation test. We compare the reconstruction error of latent midpoints for 1,000 random pairs. The latent mixup model (green) exhibits an order-of-magnitude reduction in MSE compared to the baseline (blue), indicating a significantly smoother and more linear latent manifold.

To verify the smoothness of our learned manifold, we performed a “Midpoint Interpolation Test.” We randomly sampled N=1,000N=1,000 distinct pairs from the test set. For each pair (x1,x2)(x_{1},x_{2}), we compared the decoded latent midpoint against the linear average of the inputs using the MSE:

ℒinterp=‖𝒟​(ℰ​(x1)+ℰ​(x2)2)−x1+x22‖2\mathcal{L}_{\text{interp}}=\left\|\mathcal{D}\left(\frac{\mathcal{E}(x_{1})+\mathcal{E}(x_{2})}{2}\right)-\frac{x_{1}+x_{2}}{2}\right\|^{2} (7)

Figure 7 shows the distribution of this error. The model trained with latent mixup (green) achieves an interpolation error approximately one order of magnitude lower than the baseline (blue), with the median MSE dropping from ∼6.0×10−4\sim 6.0\times 10^{-4} to ∼5.0×10−5\sim 5.0\times 10^{-5}. This confirms that latent mixup effectively enforces a continuous and well-interpolatable latent structure, ensuring robust spectral generation.

References