HOLISMOKES -- V. Microlensing of type II supernovae and time-delay inference through spectroscopic phase retrieval

We investigate strongly gravitationally lensed type II supernovae (LSNe II) for time-delay cosmography incorporating microlensing effects, which expands on previous microlensing studies of type Ia supernovae (SNe Ia). We use the radiative-transfer code ${\rm \small TARDIS}$ to recreate five spectra of the prototypical SN 1999em at different times within the plateau phase of the light curve. The microlensing-induced deformations of the spectra and light curves are calculated by placing the SN into magnification maps generated with the code ${\rm \small GERLUMPH}$. We study the impact of microlensing on the color curves and find that there is no strong influence on them during the investigated time interval of the plateau phase. The color curves are only weakly affected by microlensing due to the almost achromatic behavior of the intensity profiles. However, the lack of non-linear structure in the color curves makes time-delay measurements difficult given the possible presence of differential dust extinction. Therefore, we further investigate SN phase inference through spectral absorption lines under the influence of microlensing and Gaussian noise. As the spectral features shift to longer wavelengths with progressing time after explosion, the measured wavelength of a specific absorption line provides information on the epoch of the SN. The comparison between retrieved epochs of two observed lensing images then gives the time delay of the images. We find that the phase retrieval method using spectral features yields accurate delays with uncertainties $\small {\lesssim}$2 days, making it a promising approach.


Introduction
Over the last few years, several new and independent techniques for measuring the Hubble constant, H 0 , have been presented to address the problem of the >4σ discrepancy (Verde et al. 2019) between the value measured from observations of the cosmic microwave background (Planck Collaboration et al. 2020) and distance ladder measurements from the supernova H 0 for the equa-tion of state (SH0ES) program (Riess et al. 2016(Riess et al. , 2018; see also Freedman et al. 2019Freedman et al. , 2020. One of these techniques is to conduct time-delay cosmography on strongly lensed variable objects (Refsdal 1964) to infer H 0 through time-delay measurements, lens mass modeling, and reconstructing the lineof-sight mass distributions. The calculated time-delay distance is inversely proportional to H 0 .
The method can be applied to different variable sources. One possibility is to use quasars, which has already been done in the H0LiCOW program Birrer et al. 2019;Sluse et al. 2019;Chen et al. 2019;Wong et al. 2020;Rusu et al. 2020) in cooperation with the COSmological MOnitoring of GRAvItational Lenses (COSMOGRAIL; Eigenbrod et al. 2005;Courbin et al. 2017;Bonvin et al. 2018), the Strong lensing at High Angular Resolution Program (SHARP; Chen et al. 2019), and the STRong-lensing Insights into the Dark Energy Survey (STRIDES) collaboration . The newly established Time-Delay COSMOgraphy (TDCOSMO) organization is further testing systematic effects and expanding the lensed quasar sample for future studies (Millon et al. 2020b,a;Gilman et al. 2020;Birrer et al. 2020). A second way to conduct time-delay measurements is through strongly lensed supernovae (LSNe), which was first envisioned by Refsdal (1964) when he introduced the concept of time-delay cosmography.
Until today, only two LSNe with multiple resolved images have been observed. The first, called "SN Refsdal," is a type II supernova (SN II), discovered by Kelly et al. (2016a,b), which was lensed by the galaxy cluster MACS J1149.5+222.3. The second LSN, iPTF16geu, was discovered by Goobar et al. (2017) and is a type Ia supernova (SN Ia) at redshift 0.409 lensed by a galaxy at a redshift of 0.216. Several lensing studies on SN Refsdal have been performed (e.g., Rodney et al. 2016;Kelly et al. 2016b;Treu et al. 2016;Pierel & Rodney 2019;Grillo et al. 2016Grillo et al. , 2020Baklanov et al. 2021); iPTF16geu has been stud-ied extensively (Goobar et al. 2017;More et al. 2017;Yahalomi et al. 2017;Dhawan et al. 2018;Johansson et al. 2020) and shows strong signs of microlensing. Further studies on SNe Ia and their usability to measure time delays have been conducted by Goldstein et al. (2018), Foxley-Marrable et al. (2018), Huber et al. (2019), . The main focus of those papers has been the characterization of the perturbations introduced by microlensing.
Microlensing of SNe results from the additional lensing of strongly lensed SNe induced by the compact masses that the lens galaxy is composed of (e.g., stars). The influence of microlensing is independent for each observed image of the SN. This causes the spectrum to be magnified to an unknown amount and makes time-delay measurements difficult. Goldstein et al. (2018) and Huber et al. (2020) show that microlensing of LSNe Ia is mostly achromatic. In three of five color curves, Huber et al. (2020) find characteristic minima or maxima within the achromatic phase of SNe Ia, making time-delay measurements with LSNe Ia feasible.
With this work, which is part of Highly Optimised Lensing Investigations of Supernovae, Microlensing Objects, and Kinematics of Ellipticals and Spirals (HOLISMOKES; Suyu et al. 2020), we take a closer look at SNe II for time-delay measurements that incorporate microlensing effects. Characteristic features in the light curves and spectra of SNe II could be used to conduct time-delay measurements. With the upcoming Rubin Observatory Legacy Survey of Space and Time (LSST; LSST Science Collaboration et al. 2009) we expect to detect over 1000 LSNe II (Wojtak et al. 2019;Goldstein et al. 2018;Goldstein & Nugent 2017). This number is larger than the number of LSNe Ia expected to be observed, which is estimated to be around 500 to 900 (Oguri & Marshall 2010;Wojtak et al. 2019). Therefore, it would be a missed opportunity to not include SNe II for timedelay cosmography.
We use radiative-transfer models to investigate the potential of LSNe II for cosmographic measurements. In recent years computational capabilities have grown, allowing us to conduct more and more complex calculations for modeling SNe of different types. The study of LSNe Ia by Huber et al. (2019) used the Applied Radiative Transfer In Supernovae (artis) code, which was developed by Kromer & Sim (2009).
For this study of LSNe II, we use the Monte Carlo radiative transfer code tardis (Kerzendorf & Sim 2014), as modified by Vogl et al. (2019). The spatial specific-intensity profiles for each wavelength are microlensed by magnification maps generated with the inverse ray-tracing code gerlumph (Vernardos et al. 2015), and deformed spectra are constructed. We then evaluate the microlensed spectra, light curves, and color curves on their usability for time-delay measurements. In doing so, we investigate a new approach to inferring phases on the basis of the spectral evolution of the SN. Throughout the paper, we assume a radiation-free flat Lambda cold dark matter (ΛCDM) cosmology with H 0 = 72 km s −1 Mpc −1 , Ω m = 0.26, and Ω Λ = 0.74, as assumed by Oguri & Marshall (2010).
The structure of the paper is as follows. First, we discuss the tardis simulations and present the synthetic SN II spectra in Sect. 2. In Sect. 3 we explain our application of microlensing to those spectra and show the resulting light curves and color curves. Further, we show the impact of microlensing on the absorption features. In Sect. 4 we investigate the phase inference of SNe from the absorption wavelengths. Finally, we discuss our results in Sect. 5 and conclude in Sect. 6. Fig. 1: Schematic light curve of a typical SN II-P. The red hatched area is the region of interest for the models in this paper and is part of the plateau phase marked by the gray area.

Type II supernova models
In this section we describe the creation of the SN II models we use for the microlensing investigation in Sect. 3. We performed the calculations with the extended version of the radiative-transfer code tardis from Vogl et al. (2019) and applied it to the prototypical type II-plateau supernova (SN II-P) 1999em. The reconstructed spectra are part of the plateau phase. The characteristic temporal evolution of the light curve of a SN II-P is shown in Fig. 1 with the region of interest marked with a red crisscross pattern (Barbon et al. 1979;Anderson et al. 2014;Valenti et al. 2015;Gall et al. 2015;Utrobin et al. 2017;Hicken et al. 2017).

tardis simulations
tardis is a Monte Carlo based radiative transfer code based on the indivisible energy packet scheme of Lucy (1999Lucy ( , 2002Lucy ( , 2003, originally written for spectral modeling of SNe Ia (Kerzendorf & Sim 2014). Vogl et al. (2019) modified the code to perform spectral synthesis for SNe II and specifically applied it to spectra of SN 1999em .
The original tardis code considers Thomson scattering and bound-bound line interactions as an approximation for SNe Ia (Barna et al. 2017;Boyle et al. 2017;Magee et al. 2016Magee et al. , 2017. To account for the hydrogen-rich composition of SNe II, the modified code additionally accounts for bound-free, freefree, and collisional processes (Lucy 2002(Lucy , 2003. The modified tardis code accounts for nonlocal thermodynamic equilibrium (NLTE) effects in the ionization and excitation balance of hydrogen and, in selected cases, helium. Additionally, the thermal structure is calculated from the equilibrium of the heating and cooling of the shells.
To synthesize the spectra, the properties of the propagating energy packets are tracked and finally binned. The resulting spectra are affected by Monte Carlo noise. To improve the signal-to-noise ratio (S/N), multiple virtual Monte Carlo packets (v-packets) are emitted for each emitted or interacting packet. These additional packets propagate along randomly selected rays and get attenuated by the integrated optical depth. When the v-packet exits the computational region through the outer boundary, the emission radius r, the propagation direction cos(θ), the frequency ν, and the energy are saved. These are used to calculate the specific intensity and the impact parameter p to construct the specific intensity profiles in Sect. 2.3.
We assumed power-law density profiles and a homogeneous composition to calculate the tardis models for this study. The abundances of H, He, C, N, and O were set to their CNO cycle equilibrium values (taken from Prantzos et al. 1986). All other elements follow the solar chemical composition (Asplund et al. 2009). Due to computational expense, the code produces only parameters at a specific epoch and no temporal development. In order to provide a time-dependent evolution of spectra for the analysis of microlensing effects, we computed several spectra at different epochs and used these snapshots as the temporal evolution.

Observed and model spectra of SN 1999em
For the analysis of the effect of microlensing, we used the tardis models for SN 1999em. We selected SN 1999em as it is a prototypical SN II-P (e.g., Dessart & Hillier 2006). Additionally, it is one of the most well-observed and well-studied SNe II. Its distance has been determined several times, and it has been observed over a long period of time, resulting in many spectra (Leonard et al. 2002;Baron et al. 2004;Hamuy et al. 2001). Given the availability of a huge amount of data, it is a good candidate for SN II modeling (e.g., Utrobin 2007).
The observed spectra of SN 1999em (Leonard et al. 2002;Hamuy et al. 2001) were used as a reference for the simulated spectra as tardis only creates snapshots at specific epochs and does not perform fully time-dependent radiative transfer computations. To achieve a sensible temporal evolution of the snapshots, they were fit to the observed spectra of a specific SN.
The observed spectra are shown in Fig. 2 in comparison to the model spectra selected from extensive model grids, including the one from Vogl et al. (2020). The observed spectra have been de-reddened under the assumption of the Cardelli et al. (1989) extinction law. The color excess is set to E(B − V) = 0.08 and the total-to-selective extinction ratio to R V = 3.1, as in Vogl et al. (2019). The flux of the model spectra was calculated for the SN at a distance of 10 pc and shifted by eye to match the observed spectra. The wavelength resolution of the spectra was chosen to be about 3 Å per bin.
At early times, around day 11 after the explosion, the spectrum is characterized by a broad hydrogen Balmer series and a He i line at 5876 Å. Roughly two weeks after the explosion, Fe ii absorption lines become visible and stronger with time. All five spectra were taken during the plateau phase of SN 1999em (Leonard et al. 2002;Hamuy et al. 2001).

Specific intensity profiles
We calculated the specific intensity profiles from the tardis simulations using a similar approach as Huber et al. (2019). In contrast to these authors, we used v-packets instead of regular Monte Carlo packets. We further neglected flight-time differences between packets that are emitted at different distances to the observer. The parameters of the v-packets used for the calculation are the emission radius r, the energy , the frequency ν, the time elapsed since explosion t, and cos(θ), where θ is the angle between the position vector of the observer and the position vector of the v-packet leaving the computational domain.
In Fig. 3 the spatial distributions of the normalized intensity profiles for the six LSST filters, u, g, r, i, z, and y (LSST Science Collaboration et al. 2009), are shown. The profiles for the different filters are calculated by convolving the transmission profile of each filter with the specific intensity. Based on the plot, we chose 25 bins of the impact parameter p as an appropriate binning resolution for each epoch. The resulting intensity was normalized by the maximum of the specific intensity for each filter. The normalized specific intensity profiles are very achromatic, especially at the first three investigated epochs, which show barely any difference between the profiles for the different filters.

Microlensing of type II supernovae
With the upcoming LSST, thousands of LSNe will be detected, including SNe Ia and SNe II. The possible use of SNe Ia for time-delay measurements, and the influence of microlensing on it, has already been investigated by Huber et al. (2019) and Goldstein et al. (2018) using characteristic structures within color curves to infer a time delay between two lensing images. However, microlensing on SNe II is still an open question from a theoretical viewpoint and is the motivation for this work.
The microlensed flux was calculated according to Huber et al. (2019) using specific intensities, retrieved with the tardis simulation at the source plane, which were microlensed using magnification maps. To investigate the effect of microlensing on time-delay measurements of LSNe II, we inserted 10000 random SN positions into a microlensing map.

Microlensing magnification maps
We used the magnification maps from gerlumph (Vernardos et al. 2015) to calculate the microlensed SN spectra. These maps, which were calculated using the inverse ray-shooting technique (Kayser et al. 1986;Wambsganss et al. 1992;Vernardos & Fluke 2013), contain the magnification µ(x, y) in the source plane as a function of the Cartesian coordinates x and y in units of the Einstein radius R Ein . The characteristic scale of the map, R Ein , is defined as: where D s is the angular-diameter distance from the observer to the source, D ds the angular-diameter distance from the lens to the source, D d the angular-diameter distance from the observer to the lens, and c the speed of light. We further assumed a Salpeter initial mass function with a mean mass of the point mass microlenses of M = 0.35 M as described by Huber et al. (2020). The magnification maps are primarily defined by three parameters: the convergence κ, the shear γ, and the smooth matter fraction s = 1− κ * κ , where κ * is the convergence of the stellar component. The maps used in this work have a resolution of 20000 × 20000 pixels, corresponding to a size of 10 R Ein × 10R Ein .
In this work we consider two SN images of a lens system with different microlensing maps. Galaxy-scale strong lens systems typically have two or four images of the background source. For a two-image lensed system, the first and second images that appear are time-delay minimum (type I) and saddle (type II) images 1 , respectively. For a four-image lensed system, the first two images are of type I and the next two images are of type II. For our current study, we consider two SN images from either a double or a quad system, with the first SN image as type I and the second one as type II. Therefore, we simulated microlensing for these two SN images using two different microlensing magnification maps, which are shown in Fig. 4.
For the first image, we chose κ = 0.36 and γ = 0.35, which are the median values of type I lensing images of the OM10 catalog (Oguri & Marshall 2010). For the second image, we chose κ = 0.7 and γ = 0.7, which are the median values of type II lensing images of the OM10 catalog Suyu et al. 2020). We adopted s = 0.5 for both cases. The observed microlensed flux F λ,o (t) was computed by placing the SN into the magnification map and solving: where D lum is the luminosity distance to the source and I λ,e (t, x, y) is the specific intensity in the source plane. For the calculation we interpolated the specific intensity onto a 2D Cartesian grid of the same spatial resolution as the magnification map and integrated over the SN (Huber et al. 2019). The AB magnitudes were then calculated by:

Type II supernova light curves and color curves
In this subsection we study the effect of microlensing on the light curves and color curves retrieved from the tardis simulations and investigate if they can be used for time-delay measurements. The impact on the color curves is of particular interest as the achromatic behavior of the specific intensity profiles suggests that microlensing has little influence on the color curves. The curves were created using the filters u, g, r, i, z, and y of LSST (LSST Science Collaboration et al. 2009). The transmission of each filter was convolved with the supernova flux to obtain the light curve for that filter. The microlensed flux was calculated according to Eq. (2) with D lum = 10 pc to get the absolute magnitudes when calculating the light curves. For the non-microlensed flux, the magnification was set to µ(x, y) = 1. The luminosity distance was chosen such that further computations could be made with the flux instead of the luminosity density. The microlensed and non-microlensed AB magnitudes for the light curves were then calculated via Eq. (3). As an example for the microlensed light curve, we chose a position in the magnification map of type I lensing where the expanding layers cross a caustic with time. The position is shown in Fig. 5, with the circles representing the radius of the SN containing 99.9 % of the luminosity for the five different epochs. We see that the caustic is crossed from day 13 onward. The caustic crossing can also be seen within the light curves shown in Fig. 6 for the non-microlensed and microlensed case. We normalized the light curves by shifting the magnitude of the first epoch to zero. The light curves show that after day 11 the difference in the absolute magnitude between the nonmicrolensed and microlensed case is increased. This means that the SN is magnified by the caustic crossing after day 11. The intrinsic light curves do not have characteristic structures, such as peaks, as they include only epochs of the plateau phase of SN 1999em. The plateau phase is the part in the light curve of a SN II-P where the magnitude barely changes over a certain period of time. The light curves within the plateau phase of our investigation can be strongly affected by microlensing, which adds large uncertainties to time-delay measurements.
Even though the non-microlensed light curves are almost linear, which leads to nearly linear color curves, we wanted to show the influence of microlensing on these color curves, following Goldstein et al. (2018), Huber et al. (2019), . Due to the achromatic specific intensity profiles, the influence is expected to be low. We calculated the LSN signal at 10000 different random positions in the microlensing map. For each position, we calculated the magnitude difference, m AB,X (t i ) − m AB,Y (t i ), for two filters, X and Y. For each filter combination, we plot the median and the 1σ and 2σ deviations retrieved from the 10000 individual color curves. Four comparisons of the microlensed color curves to the non-microlensed color curves are shown in Fig. 7. Additional color curves are shown in Appendix A. It can been seen that the influence of microlensing on the color curves is small overall. Yet, these color curves cannot be used in time-delay measurements in a straightforward way because they do not show large characteristic structures that would simplify the time-delay measurements, given the possibility of differential dust extinction between the multiple SN images (Elíasdóttir et al. 2006). Therefore, to use color curves, a longer time span, covering more than the first part of the plateau phase that we consider in this work, would be helpful for inferring the time delays. We defer this to future studies.

Type II supernova spectra and absorption features
As the light curves and color curves lack significant features within the plateau phase, we investigate another approach to perform time-delay measurements on SNe II. We used the temporal evolution of the minima of the absorption lines of Hα, Hβ, and Fe ii as they are common features used for velocity measurements that can be detected in at least three epochs. The hydrogen lines can even be detected in all five of the epochs (see Fig. 8).
We applied a two-step procedure to determine the wavelength of an absorption line from noisy spectra. In the first step, we visually estimated the absorption minimum and selected a range of wavelength bins around that minimum to refine the estimate through the function signal.find peaks of the python module scipy (Virtanen et al. 2020). We adopted the output of signal.find peaks as λ init , the initial wavelength of the absorption minimum. In the second step, we fit a Gaussian to the absorption line profile to further refine the wavelength measurement, especially in the presence of noise. As the lines have individual asymmetrical P-Cygni profiles, different fitting ranges were selected for each feature. For Hα we chose 15 bins toward smaller wavelengths and 25 bins toward longer wavelengths, for Hβ 8 bins toward smaller and 18 bins toward longer wavelengths, and for Fe ii 20 bins toward smaller and 15 toward longer wavelengths. The bins are counted from λ init , and each bin has a width of 3Å. The success of the fitting is dependent on the fitting range. If the range is too small, then the fit might fail. Noisy spectra in particularly can sometimes have two small minima within the absorption feature. If the range is too large, then the fit might shift to longer wavelengths as a result of the asymmetric P-Cygni line shape. The optimal wavelength range to use for the Gaussian fit was determined through simulations of noisy spectra. We applied different fit windows onto the noisy spectra and chose the window that best reproduced the known absorption minimum of the model spectrum. 2 The result of the Gaussian fit yields the wavelength of the absorption line, as detailed below.
The Gaussian fit was performed using the optimize.curve fit function of the python module scipy, a least-squares minimization using a Levenberg-Marquardt optimizer with a fitting function: The input x is the wavelength range over which the fit is computed. The parameter A is the square root of the amplitude, M the mean of the Gaussian, σ the standard deviation, and h the offset. If the fitting failed, which is indicated by a determined minimum that is larger than λ init +60 Å (outside the fitting range) or smaller than λ init − 60 Å, the fitting was repeated with a fit range shifted by five bins, which equals 15 Å, to longer wavelengths. In these few cases, it is the effect of microlensing that deforms the spectral feature and shifts it to the red or blue once the SN crosses a (a) (b) (c) (d) Fig. 7: Color curves of the non-microlensed and microlensed spectra with median, 1σ, and 2σ ranges obtained from 10000 random positions in the type I lensing magnification map. In addition, the bottom part of each panel (7(a) to 7(d)) shows the deviation of the microlensed color from the non-microlensed color with 1σ uncertainties. caustic. The shift to redder wavelengths cannot be compensated for by the fitting method itself, making this additional calculation necessary. In Fig. 9 examples of fits to the non-microlensed absorption lines are shown. We see in these plots that some fits, especially for later times, tend to the bluer or redder part of the absorption line due to the asymmetrical line profile and fitting range. As this shift is always in the same direction, this does not affect the following phase retrieval method. The fitting was performed for every position in the microlensing map, for each of the five epochs and for each absorption line in the microlensed spectra without noise. The temporal evolution of the absorption minima of Hα, Hβ, and Fe ii is shown in Fig. 10 for the two magnification maps of Fig. 4. The plots include the median of all determined minima of the positions in the microlensing map and the 1σ and 2σ uncertainties. The temporal evolution of the wavelengths of the absorption minima of the spectra without microlensing is shown in yellow for reference. The evolution is similar for the two images, but microlensing tends to introduce larger uncertainties for the second lensing image. In addition to the noise-free case, we repeated the calculations assuming a S/N per pixel of 10 for the spectra. To simulate a noisy spectrum, we added Gaussian noise to the flux of each wavelength bin of the microlensed spectrum; in other words, we drew a random Gaussian number with standard deviation given by σ G = flux/(S/N) and added it to the flux. We simulated a noisy microlensed spectrum for each position in the microlensing map and applied a Savitzky-Golay filter to make the fitting procedure easier. We obtained the distribution of the absorption minima shown in red in Fig. 10. In all panels of Fig.  10 (particularly the top panels for Hα images), the uncertainties due to microlensing (gray bands) are smaller than the uncertainties due to noise in the spectra (red bands) when S /N = 10.
Furthermore, we investigated the correlation between the absorption minima of the three lines for different S/N in Fig.  11. We show the plots for S /N = 10, 20, 30, and ∞ at day 22. We consider these three additional S/N values besides the noiseless case as feasible for future observations. We also use these S/N and day 22 for the phase retrieval in the following section. In contrast to the noisy cases, the noiseless case has a very tight data distribution, and therefore we only show the 2σ contour. For S /N = ∞, the absorption minima of Hα, Hβ, and Fe ii are strongly correlated within the 1σ and 2σ contours. For S /N = 30, the 1σ contour already shows no correlation, whereas in the 2σ contour slight correlations are still present. As the S/N gets smaller, the 1σ stays uncorrelated and the 2σ contour become increasingly more uncorrelated.

Supernova phase inference from spectra and time-delay measurement
In this section we investigate the possibility of inferring phases from the wavelengths of the absorption minima of the spectral lines. The basic idea is based on Fig. 10, which shows that the wavelengths of the absorption minima evolve over time. When we detect a LSN system based on its first appearing image, which we denote as SN-A, we can obtain spectra of this SN-A image at different epochs to determine the evolution of the absorption features (as shown in Fig. 10). A subsequent spectrum taken at time t B of the next appearing SN image, which we denote as SN-B, can then be used to infer the phase of SN-B based on the measured wavelength (i.e., the time t A at which image SN-A had the same absorption wavelength value as that measured in the SN-B spectrum). The difference between these two times, t B − t A , then gives the time delay between the SN-B and SN-A images. Below, we describe in more detail how we quantitatively inferred the phase and time delays, taking microlensing and noise in the spectrum into account.
Here we consider the scenario where we have obtained spectra of SN-A at five epochs (as in the previous sections), denoted by t i , where i = 1, 2, ..., N epoch and N epoch = 5. The spectra have the same S/N. Through the microlensing simulations described in the last two sections, we obtained the distribution of absorption line wavelengths, accounting for microlensing and noise as shown in the left panels of Fig. 10. In particular, for the microlensing of the SN-A image, we adopted the magnification map with κ = 0.36 and γ = 0.35 since the first appearing SN-A image must be a time-delay minimum image. By sampling Fig. 10: Temporal evolution of the absorption minima of the Hα, Hβ, and Fe ii lines for non-microlensed and microlensed spectra using the two different magnification maps of Fig. 4 (left: Fig. 4(a); right: Fig. 4(b)). For the microlensed spectra without noise (labeled S/N = ∞), the median of the 10000 positions and the 1σ and 2σ ranges are shown in black and gray, respectively. The temporal evolution of the absorption lines with added noise of S /N = 10 is shown in orange and red. 10000 positions on the microlensing map and thus simulating 10000 noisy microlensed spectra for the given S/N at each t i , we obtained the distribution of absorption line wavelengths λ(t i ) following the absorption line fitting method described in Sect.

3.3.
We approximated the distribution of λ(t i ) at each epoch as a Gaussian with mean λ d,i and standard deviation σ d,i . To obtain a continuous distribution of λ d (t) and σ d (t) for any time t, we linearly interpolated between the discrete λ d,i and σ d,i values. This allows us to describe the likelihood of the data d A (i.e., the spectra and microlensing maps of SN-A) for a given time t and absorption wavelength λ as The solid 1σ contours in the 2D histograms contain 68% of the data and the dashed 2σ contours 95%. The 1D histograms contain the 1σ (vertical solid lines) and the 2σ (vertical dashed lines) ranges. For S /N = ∞ we only show the contour containing 95% of the data and the 2σ range. Above the 1D histograms the mean wavelength of the absorption minimum with 1σ uncertainties is indicated. While there are significant correlations in the absorption minima of the different lines in the case of noiseless spectra, the correlations become weak or negligible with noisy spectra.
where C A is the normalization constant.
Next we consider a spectrum taken of the SN-B image at t B . Similar to the procedure described above for the SN-A image, we repeated the microlensing simulations to produce 10000 microlensed spectra of SN-B for a given S/N to account for uncertainties associated with microlensing and noise in the spectrum. The main difference from the procedure of SN-A is that the microlensing map for SN-B is for the time-delay saddle image with κ = 0.70, γ = 0.70. This is appropriate since at least one of the next appearing SN images is a time-delay saddle. Furthermore, considering the saddle image allows us to probe uncertainties due to both types of lensing images (time-delay minimum and time-delay saddle images) given that the microlensing maps are sensitive to the image types. From the distribution of 10000 microlensed noisy realizations of the SN-B single epoch spectrum, we can measure the mean wavelength λ B and standard deviation σ B . We used the wavelength measurement λ B ± σ B to infer the times in SN-A that are consistent with these values.
Quantitatively, we are interested in obtaining the probability distribution of time t of SN-A given the data d A and mea- surement d B (of λ B ± σ B ) -that is, P(t|d A , d B ) -which can be expressed by where λ are the absorption line wavelengths associated with SN-A. Using Bayes' rule, we can rewrite the integrand as: where we have separated dependences (d A and t do not depend on d B ). The expression for P(d A |λ, t) is given by Eq. (5). Through our Gaussian approximation of the measurement d B , we can write Given that our data d A is fixed, P(d A ) is a constant. Assuming a uniform prior on t (i.e., P(t) is constant), we can simplify Eq. (6) to where C AB is a constant. We can compute P(t|d A , d B ) through importance sampling (e.g., Lewis & Bridle 2002). We first generated samples of P(d A |λ, t) by (1) randomly drawing a value of t, assuming they follow a Gaussian distribution, which we denote as t samp , (2) computing λ d (t samp ) and σ d (t samp ), and (3) using these sampled values to draw a random Gaussian deviate, λ samp . We repeated steps (1)-(3) to obtain 31981 samples of (t samp , λ samp ) (or 19980 samples for the specific case of the Fe ii absorption line). The number of samples was the result of setting the amount between the third and fourth epoch to 10000 samples. The number of samples within the other epochs resulted from choosing the same density of samples along the wavelength axis. Each of these samples was then weighted by P(λ samp |d B ) as given by Eq. (8). The weighted distribution of t finally yields P(t|d A , d B ).
As an example, we applied this procedure with t B set equal to the fourth epoch at day 22.0 of the SN-A image modeled with tardis. We chose this epoch as it is the central epoch of the temporal evolution of the Fe ii absorption line and enabled us to evaluate all three absorption lines of interest, Hα, Hβ, and Fe ii. For the noise addition, we again selected the three different S/N values of 10, 20, and 30.
The retrieved phases of all three absorption lines are shown in Fig. 12 Table 1: Weighted average and the 1σ uncertainties of the phase retrievals considering all five epochs of the first image. We list the results for the three absorptions lines (Fe ii, Hα, and Hβ) for the three different S/N values (10, 20, and 30) and for the noiseless case, indicated by S /N = ∞. The last row shows the uncertainty on the phase when using all three lines, assuming the noise and the microlensing in these lines are uncorrelated. recovered very well by our phase retrieval method. In particular, assuming a S /N = 20 or higher, the histogram shows a sharp peak centered around the input phase value. We additionally combined the three separate phase retrievals in Fig. 12(d). This assumes the best case scenario with negligible correlation in the noise and microlensing, which means that we can simply multiply the three distributions with one another. The assumption of uncorrelated data is justified based on the correlation plots in Fig. 11, which show that in the case of noisy spectra the absorption minima are only very weakly correlated. Only the noiseless spectra show correlations, which we can safely disregard since S /N = ∞ is a limiting case, only included in this study to give an upper limit for the accuracy. Additionally, its median and width are very similar to the case of S /N = 30, suggesting that the correlation does not matter much for the present work. The peak around the epoch that should be retrieved is sharper than for the single line retrievals. Another possibility for reducing the uncertainties would be to combine several spectra of the SN-B image, if they are available, but taking possible correlations introduced by microlensing into account.
In Table 1 we list the weighted average and the corresponding 1σ deviations of the phases retrieved for the different absorption lines and their combination. For Fe ii and Hβ the averages deviate more than the average of Hα from the input value of 22.0 days. For Hα the input value is recovered within its 1σ uncertainties. For all considered S/N, we recover the input value when combining the phase retrieval results of the three absorption lines.
Even though a S/N of 10 already yields a good result for the phase retrieval with uncertainties smaller than 2.5 days, the uncertainties of S /N = 20 and 30 are both smaller than 2 days. The difference in 1σ of S /N = 20 and 30 is marginal or even negligible, with S /N = 20 also being very close (within 0.2 days in uncertainty) to the noiseless 1σ deviation, suggesting that a S/N of 20 would be sufficient to apply our phase retrieval method.
Our method also does not require very high resolution spectra. Well-calibrated spectra with resolution R 500 are sufficient for determining the relative phase of SN images to fractions of a day since the line profiles of SNe are broad.
To further explore the effect of microlensing on the phase retrieval, we investigate a scenario with extremely high micromagnification using two additional magnification maps in Appendix B. We find the phase retrieval to be only negligibly affected by microlensing even in this high microlensing scenario since we obtain the same retrieved phases within 0.2 days for the original and the high-magnification microlensing scenarios.
Furthermore, we consider a case where the fourth epoch that we want to retrieve with the second image is not measured within the first image to check how well the method works if the available data are limited. Another reason for this scenario is that with real measurements we might not have data of the second image covering the same epochs as for the first image. For this case we only used the third and fifth epoch and interpolated linearly between them to create samples as explained previously in this section. We justify the linear interpolation by the linear behavior of the absorption minima of the last three epochs as shown in Fig. 10. In particular, Hα follows a linear trend to a high degree of precision. We thus expect Hα to give the best results in this case. The resulting phase retrieval histograms are shown in Fig. 13 with the weighted average and the 1σ uncertainties for the different absorption lines, their combination, and for the different S/N realizations shown in Table 2. As expected, the recovered weighted average of the phase using the Hα absorption line agrees well with the exact value of 22.0 days. Even though the precision is almost unchanged, the accuracy is strongly biased in the case of Fe ii and Hβ. This most likely results from the nonlinear evolution of the last three bins, which is not as linear as for Hα.
In general we detect three sources for bias: microlensing, fitting ranges, and the cadence of the observations. The cadence is an especially important factor for the first image. At the moment we can only estimate the strength of each of these bias contributions from our results. How much each bias affects the measurements will depend on the details of actual data once they are available. Combining the phase inference from multiple spectral features can potentially reduce the bias compared to the measurements for a single line. As illustrated in Fig. 13 and Table 2, when we multiply the phase-retrieval histograms together, the bias cancels out since the underprediction of the phase from the Hβ absorption line counterbalances the overprediction of the phase from Fe ii. Depending on the spectral features used, such a perfect cancelation of the biases might not occur. Therefore, when applying our method to real data, we recommend investigating individual spectral features through simulations to quantify potential biases, as we have done here.
For accurate phase retrievals, we ideally need as many measured epochs of the first image as possible to ensure that we cover the range of phases where we might be able to measure the second image and use it in the phase retrieval. This further ensures that linear interpolations only cover short periods, which minimizes the bias induced by a long gap between epochs. From the example in Fig. 12 we can see that a cadence of around five days for the first image already reduces the bias significantly compared to a cadence of approximately ten days as used in Fig.13. Having as many epochs as possible for the first image also ensures that we do not need to extrapolate (which would lead to another source of bias) to cover an epoch of the second image that could be missing in the measurements of the first image. As the combination of the phase retrieval of several different absorption lines can lead to the canceling of the bias, we recommend using at least three different absorption lines and also considering them individually, which gives a general hint of the amount of bias to expect. 22.5 ± 2.3 days 22.7 ± 1.9 days 22.7 ± 1.8 days 22.8 ± 1.7 days Hα 21.9 ± 2.0 days 21.9 ± 1.5 days 21.9 ± 1.4 days 21.7 ± 1.3 days Hβ 21.6 ± 2.1 days 21.5 ± 1.8 days 21.5 ± 1.8 days 21.4 ± 1.7 days Fe ii + Hα + Hβ 22.0 ± 1.4 days 22.0 ± 1.0 days 22.0 ± 1.0 days 21.9 ± 0.9 days Table 2: Weighted average and the 1σ uncertainties of the phase retrievals considering the third and the fifth epoch of the first image. We list the results for the three absorptions lines (Fe ii, Hα, and Hβ) for the three different S/N values (10, 20, and 30) and for the noiseless case, indicated by S /N = ∞. The last row shows the values when using all three lines, assuming the noise and the microlensing in these lines are uncorrelated.

Discussion
As of today, only one strongly lensed SN II has been observed, namely SN Refsdal, which was first detected by Kelly et al. (2016a,b). With its five detected and resolved images, time-delay cosmography with a SN II is being conducted for the first time. The system has so far primarily been used to crash test and validate lens modeling codes' predictions of the place and time of the reappearance of the last trailing SN image. The short time delays between the first four appearing images have been measured by Rodney et al. (2016) based on their light curves and by Baklanov et al. (2021) using model light curves, and there are ongoing efforts to measure the time delay of the last trailing SN image using light curves (P. Kelly, private communications).
Type Ia supernova iPTF16geu, discovered by Goobar et al. (2017), has four spatially resolved images from which timedelay measurements have been conducted. The delays have been calculated through two independent channels: from photometry of consecutive image appearances by Dhawan et al. (2020) and from spectral features by Johansson et al. (2020). The authors mention that microlensing effects (More et al. 2017;Mörtsell et al. 2020) on the light curves cause inaccuracies in the timedelay measurements. Therefore, they used the spectra for the first time to retrieve time delays. They fit the blended spectra of three of the four SN images and, separate from that group, the spectrum of the remaining SN image with two different sets of templates: the template spectra of Hsiao et al. (2007) and the spectral energy distribution template of SN 2011fe (Amanullah et al. 2015). Significant features of the fits were then used to measure the time delay between the image groups. In the case of iPTF16geu, the spectroscopic time delay is not as precise as the photometric time delay, but the spectra seem to be less sensitive to microlensing than the color curves . In contrast to our new method, which takes advantage of the temporal evolution of the blue shift of absorption lines, the method by Johansson et al. (2020) needs very good sets of spectral templates. These only exist for SNe Ia.
In the context of our new method of retrieving time delays from absorption lines of SNe II, the studies of iPTF16geu indicate that light curves and color curves with strong features are in general more promising than spectral lines for the calculation of time delays. Microlensing has a minor impact on the color curves, especially if the intensity profiles appear to be mostly achromatic. The spectral approach is also highly dependent on the observational data quality. If the measurements are dominated by noise, features and absorption lines may become undetectable, making spectra unusable for time-delay measurements. Nevertheless, in cases where the color curves do not show useful nonlinear structures for delay measurement, which is the case in the early plateau phase of SNe II-P, we demonstrate that the phase of a spectrum can be accurately inferred. The uncertainties of the retrieved phase using a single absorption line lie within ±2 days (provided several spectra of one of the other LSN images in the system have been obtained to determine the evolution of the absorption line). Considering that delays of galaxy-scale LSNe typically range from days to weeks, our new method is very promising for obtaining delays with <10 % uncertainty for systems with delays longer than 20 days.
By comparing our LSNe II-P work to previous LSNe studies (e.g., Goldstein et al. 2018;Foxley-Marrable et al. 2018;Huber et al. 2019Huber et al. , 2020, we find that LSNe Ia and II-P have separate advantages for measuring time delays. The homogeneity of SNe Ia facilitates the comparison with templates, making it easier to retrieve time delays with them. In contrast to SNe Ia, the plateau phase of SNe II-P lacks characteristic maxima or minima within the light and color curves, which makes this phase not ideal for time-delay measurements using color curves. However, this phase showed less "chromatic" evolution within the intensity profiles than SNe Ia (Huber et al. 2019, making LSNe II-P less susceptible to microlensing effects within the early plateau phase. Therefore, our newly developed time-delay measurement method using the absorption lines of the spectra of LSNe II-P is a promising way of conducting cosmography in the LSST era (LSST Science Collaboration et al. 2009) as more LSNe II than LSNe Ia are expected to be detected (Oguri & Marshall 2010;Wojtak et al. 2019;Goldstein et al. 2018;Goldstein & Nugent 2017). Methods of conducting time-delay measurements during the plateau phase of SNe II-P will be most useful as this is the time range where observations will be easiest to perform due to the brightness and duration of that phase.

Conclusion
In this work we investigated strongly lensed SNe II and the possible impact of microlensing on time-delay measurements. We used model spectra of SN 1999em from the tardis code modified by Vogl et al. (2019). These were microlensed using magnification maps generated with the code gerlumph (Vernardos et al. 2015). We then investigated the effect of microlensing on light and color curves. We find that the color curves are mostly unaffected by microlensing since the specific intensity profiles in different bands are very similar during the plateau phase. During the investigated time span, the light curves, and therefore also the color curves, are almost linear, with no strong characteristic features that could be used as easily as for SNe Ia around the same epochs for time delay retrieval. Further investigation of time-delay measurements using SNe II-P color curves is deferred to future studies. We thus developed a new method for calculating the time delay of two SN images from the absorption lines in their spectra. For both images, we simulated spectra that were distorted by microlensing and had noise with different S/N values. We discovered that the phase can be determined from a single spectral feature with a 1σ uncertainty of <2 days and the weighted average reproducing the input phase value within the 1σ uncertainty, assuming a realistic S/N of 20. If several measurements of different absorption lines are combined, the uncertainty can be further minimized and the accuracy improved as possible biases can cancel one another out. By neglecting any correlations between the errors of the different absorption lines, the 1σ uncertainty can be minimized to ± 1.0 days for S /N = 20. Given the complementarity between the spectral and photometric delay inference, we advocate for the combination of both approaches in future measurements of LSNe.

Appendix B: Microlensing maps with high magnification
In addition to the already investigated microlensing case using the magnification maps with κ = 0.36 and γ = 0.35 (image I) and κ = γ = 0.70 (image II) from Fig. 4, in this section we investigate a high magnification scenario. Therefore, we applied our method following the steps explained in Sects. 3.3 and 4 to gerlumph magnification maps with κ = γ = 0.43 (image I) and κ = 0.57 and γ = 0.58 (image II) as shown in Fig. B.1. These two microlensing maps correspond to the high magnification case of the 1σ contours of LSNe in the OM10 catalog for images I and II, respectively ).
(a) Type I microlensing map (b) Type II microlensing map We show the temporal evolution of the absorption minima of Hα, Hβ, and Fe ii of the high magnification case in Fig. B.2. Compared to the case discussed in Sect. 3.3, there is barely any visible difference.
The histograms of the phase retrieval using this second microlensing scenario are shown in Fig. B.3. In this scenario the input epoch is also recovered very well, especially for S /N = 20 and higher S/N cases.
In Table B.1 we list the weighted average and the 1σ deviations of the retrieved phases. Compared to the averages and uncertainties of the originally investigated microlensed spectra in Table 1, we see that this second case has the same or smaller uncertainties and retrieves the phase with the same accuracy. As the second scenario microlensing maps correspond to a high magnification case, we expected higher or at least equal uncertainties. The smaller uncertainties can be explained by fluctuations caused by the spectral fitting procedures. As the absorption lines have the asymmetric P-Cygni profile, the fitting ranges influence the phase retrieval uncertainties by 0.2 days. The chosen fitting ranges discussed in Sect. 3.3 appear to work better for the second scenario than for the first scenario, especially for Hβ, resulting in smaller phase retrieval uncertainties in some cases. In general, both investigated microlensing cases have comparable phase retrieval results.  21.8 ± 2.2 days 21.9 ± 1.5 days 21.9 ± 1.3 days 21.9 ± 1.2 days Hβ 21.5 ± 2.2 days 21.7 ± 1.7 days 21.7 ± 1.6 days 21.7 ± 1.6 days Fe ii + Hα + Hβ 22.0 ± 1.3 days 22.0 ± 1.0 days 22.0 ± 0.9 days 22.0 ± 0.8 days Table B.1: Phase retrieval under high microlensing magnifications. We list the weighted average and the 1σ uncertainties of the phase retrievals using the three absorption lines (Fe ii, Hα, and Hβ) for the three different S/N values (10, 20, and 30) and for the noiseless case, indicated by S /N = ∞. The last row shows the uncertainty on the phase when using all three lines, assuming the noise and the microlensing in these lines are uncorrelated.