# Appendix B: Advanced data reduction

## B.1. Data reduction pipeline

We first run the raw science data through the standard MUSE pipeline (Weilbacher et al. 2020), using the development version<sup>25</sup>. Individual exposures are processed by the *scibasic* recipe with corresponding daily calibrations (flat fields, bias, arc lamps, twilight exposures) and geometry table (one per observing run) to produce the pixel tables (hereafter referred to as *pixtable*).

For the MXDF, we applied two changes related to the dark subtraction and the overscan parameters. To improve the bias subtraction, the overscan parameter was set to *vpoly:15,1.00001,1.00001* for all the recipes using that parameter (Weilbacher et al. 2020 Sect. 3.1). A sequence of long darks from June to August 2018 was used to produce a master dark, and then to derive a smooth model for each CCD as described in Weilbacher et al. (2020) Sect. 3.2.

The pipeline recipe *scipost* is then used to perform astrometric and flux calibrations on each pixable. For the MXDF, we computed a median standard response for each GTO run, using the standard exposures observed during that run and excluding those taken under non-photometric conditions.

The *scipost* recipe is run a first time in a "fast" mode (no sky subtraction, no Raman contamination correction) to produce images that are used to compute the centering offsets. Those offsets are calculated with the PSF fitting algorithm described in Bacon et al. (2017), Sect. 5.1, relative to the HST ACS images from the XDF data release (Illingworth et al. 2013).

The *scipost* recipe is then run a second time with sky subtraction, Raman correction and self-calibration but without applying the astrometry (*-astrometry=false*) to produce pixtables that we can use for the superflat (see Sect. B.4). The same recipe is run again to produce datacubes from these pixtables, using the astrometry and the offsets computed previously. In this case *scipost* will detect that the other steps have been done in the previous run so it will only need to produce the datacubes. The goal of this process is to minimize the computation time since we have to process hundreds of exposures.

For the different steps that require it (self-calibration, sky-subtraction and post-processing sky residual subtraction ZAP software) we use a source mask that we computed on the combined datacube from a previous version of the reduction.

## B.2. New version of the self-calibration

We modified the recipe *scipost* to include the self-calibration algorithm. This algorithm was previously developed in the Python package MPDAF (Piqueras et al. 2019). Since this algorithm uses the average sky level as a reference to calculate the correction, we had to interrupt the *scipost* recipe before the sky subtraction, save the pixtables, reload them into Python and run the self-calibration, save again and resume processing with the MUSE pipeline. The overhead was huge, especially because the pixtables are very large files (8.6 GB) and saving and loading hundreds of them to remote storage takes time. Since the algorithm was already implemented as a C extension, moving it into the MUSE pipeline was an obvious choice to optimize processing, and it was also a great opportunity to make it available to the community.

Before porting the code to the MUSE pipeline, we developed a new version of the algorithm. When comparing the results

<figure><img src="image_42.png" /><figcaption>**Fig. B.1.** Illustration of the new self-calibration. Left panel: background subtracted white-light image for one exposure after the *scipost* recipe. Center panel: the same processed with the old version of the self-calibration. Right panel: with the new version. Flux units are $10^{-20} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ Å}^{-1}$.</figcaption></figure>

of the first version used for DR1 with the flat field correction method *CubeFix* developed by Cantalupo (in preparation; see Cantalupo et al. 2019 for a description), it was clear that we could do better. The new algorithm is inspired by both *Cube-Fix* and the old version, and is described in more detail in Weilbacher et al. (2020). While *CubeFix* operates on datacubes, self-calibration operates directly on pixtables and could therefore be integrated into the *scipost* recipe.

## B.3. New version of ZAP (2.1)

The ZAP software (Soto et al. 2016) is a PCA-based tool developed to remove sky residuals in MUSE empty fields left by the imperfect sky subtraction pipeline process.

With the new version of self-calibration giving better results, it became clear that some of the artifacts we were observing in the reduced datacubes were due to ZAP. It introduced spatial variations for some pixels, which were visible in the white light images, and was also responsible for wiggles in the red part of the spectra. We corrected these problems with a new version of ZAP (version 2.1), where we made several changes: First, for the spatial variations we replaced the custom implementation of the PCA with the implementation from scikit-learn (Pedregosa et al. 2011), which is also much faster. Second, for the wiggles in the red part, we used the median continuum filter instead of the weighted one (which is a median filter weighted by the median sky level) and we increased the width of the filter window to 300 pixels. The default values were changed as well in version 2.0 for the window size and 2.1 for the filter type.

Third, the new version uses only one sky segment by default, which means that the cube is no longer split along the wavelength axis. Originally ZAP used 11 segments, whose purpose was to have coherent groups of sky emission lines, with a smaller number of eigenvalues per segment. It also allowed the computation to be parallelized. But the segments were also responsible for continuum oscillations, and made the choice of the number of eigenvalues per segment very difficult and very sensitive. With only one segment the performance of the sky subtraction is much better, thanks to the higher correlation between sky lines on the whole wavelength range. It also easier to control the power of the sky subtraction since the number of eigenvalues used for the sky reconstruction must be chosen or determined automatically for each segment.

Fourth, another issue was that sometimes ZAP suppressed signal for very bright lines, such as bright [O II] emitters. With the new version this is less likely to happen, thanks to the use of only one sky segment that provides a better estimation of the sky signal.

<footnote n="25"><sup>25</sup> Since v2.8.3 all development features used here are part of the public version.</footnote>
<figure><img src="image_43.png" /><figcaption>**Fig. B.2.** Example of improved performance of the new ZAP version. The top, middle and bottom rows display respectively the non-zapped, old zapped, and new zapped spectra of an [O II] emitter in the MXDF. From left to right, the columns display: the full spectral range, the red part of the spectra, the [O II]λλ3726,3729 emission lines and the difference with the no ZAP spectrum. Flux units are 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹.</figcaption></figure>

As shown in Fig. B.2 this new version of ZAP has significantly improved the accuracy of sky subtraction and reduced the risk of capturing signal from bright emission lines. However, it revealed another issue: the flux at some spaxels are now sometimes rising or falling in the far red (above 8800Å; see the bottom left panel of Fig. B.3). This appears to be an instrumental effect, likely related to 2nd-order effects in the twilight-skyflat exposures used in the first step of the data reduction, and the old version of ZAP was removing this signal with its two last segment, instead of the sky emission in this wavelength range. Fortunately, this problem is now corrected by the superflat process described in the next subsection (see the bottom right panel of Fig. B.3).

Another change with respect to DR1 is that previously we were running ZAP on the individual exposures to get the best sky subtraction. With the superflat correction done on the individual datacubes (see the next section) this appeared not anymore necessary and we could run ZAP only once on the combined cube, which saves some computation time.

## B.4. Superflat

In DR1, we masked the inter-stack holes in each individual exposure before combining them (Bacon et al. 2017, Sect.3.1.3). However, in addition to producing lower S/N region in the final datacube, the masking was never perfect and left some artifacts. To improve the process, we build a *superflat* for each datacube, by combining many exposures where we cancel the rotation and dithering. This is done by adding a *MUSE超级 flat_pos* environment variable in the MUSE pipeline, which allows the RA, DEC, and DROT keywords to be overridden by the values of the exposure for which the superflat is built. This way the "instrumental grid" is the same for all the datacubes and the sources move because of the rotation and dithering of the field. As we have sparse fields, combining those exposures with sigma-clipping removes the astrophysical signal and produces the superflat, where only the instrumental residuals remain.

For each exposure a superflat is built with 36 exposures from the same observing run, possibly supplemented with additional

<figure><img src="image_44.png" /><figcaption>**Fig. B.3.** Example of superflat applied to a single exposure. The top and bottom rows display respectively the white light and red part (9000-9350Å) images. The left and right columns show the correction brought by the superflat (central column). Flux units are 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹.</figcaption></figure>

exposures from the previous or next run when there were not enough exposures in a single run. The superflat is subtracted from the science exposure<fnref n="26" />. To minimize the computation time, we build the datacubes tailored to each exposure from the pixtables that have been computed previously and where the self-calibration, sky-subtraction and Raman correction have already been done. Therefore, it is sufficient to produce the datacubes for each exposure after applying the astrometry with the modified RA, DEC and DROT. Even with this shortcut, applying the superflat correction to hundreds of exposures is computationally intensive.

As shown in the top panels of Fig. B.3 the superflat gives excellent results on the white light image. Most of the residuals and holes in the inter-stacks are corrected, which means that we no longer need to mask these areas. We also note that red effect in the 9000–9350Å band (see the previous section) is also corrected very well (Fig. B.3 bottom panels).

## B.5. *musered*

To deal with the complexity of this data reduction, given the number of exposures and the number of reduction steps, we developed a Python package called *musered* on top of the MUSE pipeline. This package uses `python-cpl`<fnref n="27" /> to run the pipeline recipes, and it uses a database to gather information about the raw files and to keep track of everything that has been done, for each file, each recipe, and each version of the reduction. *musered* takes care of file associations when running a recipe, to find the correct calibration files, and provides a convenient command line interface to run the various recipes, either from the pipeline or from custom Python recipes. It also makes it possible to inspect what was done for each exposure, when it was processed, and with which parameters and which calibration files.

The MXDF data set consists of 69 different nights and 373 exposures. It was then very time consuming to carefully inspect each exposure and each calibration frame. We therefore relied on several measurements and plots to identify problematic expo-

<footnote n="26">²⁶ After a few experiments, we found that subtracting the superflat from the science exposure produced better results than conventional division as it is typically done in imaging.</footnote>

<footnote n="27">²⁷ https://github.com/olebole/python-cpl</footnote>
sures, which sometimes led to the exclusion of a few problematic
calibration sequences<fnref n="28" />. In this case, the sequence is flagged and
another sequence from the previous or next night is used. The
entire process costs 397 min of computing time per exposure on
an 80-core workstation, of which 90% was due to the superflat
computation. It took 103 days to completely reduce the MXDF
data set.

B.6. *Post-processing*

Before combining the datacubes of all individual exposures, we
identify those that should be discarded due to poor quality. The
spatial PSF of each exposure is estimated with the muse-psfr
tool, as described in Sect. 4.2.1. Twenty-two exposures with
a PSF at 7000 Å greater than 0.78 were rejected. The remain-
ing 338 datacubes were combined using a 5σ sigma clipping
scheme. The residuals left by sky subtraction were then removed
using the new version of ZAP, as described in Sect. B.3. An addi-
tional check was performed on the resulting datacube to confirm

that no offset was left in the background level at each wavelength
plane.

As shown in Weilbacher et al. 2020, Sect. 4.6, the pipeline
propagated variance of the datacube does not take into account
the additional correlated noise due to the interpolation process.
In DR1, we replaced the propagated variance of each indi-
vidual datacube with a noise estimate derived from the dat-
acube after masking the light sources, which is then rescaled
to account for the impact of the correlated noise (Bacon et al.
2017, Sect. 3.5.1). The limitation is that this noise estimate is
not valid for bright sources that are not sky-dominated. Another
inaccuracy is that this noise estimate is by construction constant
across the field of view and thus does not account for possible
flat field variation. Thus, to address these limitations, we imple-
mented a different scheme for MXDF: we computed the median
ratio between the propagated variance and the variance estimated
on the datacube itself after adjusting for correlated noise. This
ratio has a slight evolution with the wavelength from 2.45 in the
blue to 2.55 in the red. The propagated variance is then simply
multiplied by this factor.

<footnote n="28">²⁸ 1 BIAS, 1 FLAT, 1 WAVE and 1 TWILIGHT exposures were discarded.</footnote>
# Appendix C: ODHIN

For more details of this whole procedure please refer to chapter 3 of Bacher (2017).

## C.1. Model

### C.1.1. Notations

Let **M** be a matrix composed of the elements $m_{i,j}$ and **v** be a vector. Let diag($s_j$) be a diagonal matrix with diagonal elements $s_j$. $||M||_2$ is the Frobenius norm of matrix **M**. Tr(**M**) is the trace of matrix **M**. $M^+$ is the pseudo-inverse of matrix **M**. Hadamard (element-wise) product between a matrix **A** and a matrix **B** is noted $A \circ B$.

Let **Y** be the data matrix from MUSE of size $n \times \lambda$ where $n = pq$ is the number of spaxels (p,q are the two spatial axis) and $\lambda$ is the number of spectral bands. Let **Z** be the data matrix from HST of size $N \times \Lambda$ with $N$ the number of image pixels, $\Lambda$ the number of filters ($\approx 4$). $\Lambda \ll \lambda$ and $n \ll N$.

Let **X** be the unknown field with a high spatial and spectral resolution of size $N \times \lambda$.

### C.1.2. Model

Data observed by MUSE is assumed to be the addition of spatially degraded **X** data field and measurement noise **N**. Spatial degradation is the multiplication of a convolution operation (by the FSF) and a subsampling operation. The spatial PSF (FSF) varies with the wavelength. We then have the following model:

$$ \text{For } 1 \le l \le \lambda, \quad \mathbf{Y}_l = \mathbf{B}_l \times \mathbf{X}_l + \mathbf{N} \quad (C.1) $$

with $\mathbf{B}_l$ the spatial degradation matrix (subsampling and convolution), which can be written as $\mathbf{B}_l = \mathbf{S} \times \mathbf{C}_l$ where $\mathbf{C}_l$ is the convolution matrix by FSF at wavelength $l$ and $\mathbf{S}$ is the subsampling matrix.

The HST observation, composed of images from different filters, is considered the result of a spectral degradation of **X**:

$$ \mathbf{Z} = \mathbf{X} \times \mathbf{A} \quad (C.2) $$

with **A** the matrix of the HST filter responses.

Let **H** $\in \mathbb{R}^{N \times k}$ be the segmentation matrix at HST resolution. This segmentation is obtained using (Rafelski et al. 2015). Elements of **H** are defined as

$$ h_{ij} = \begin{cases} 1 & \text{if the source } j \text{ is present on pixel } i, \\ 0 & \text{otherwise.} \end{cases} $$

It should be noted that we assume that MUSE and HST data have already been well aligned spatially. We also note that we do not take the LSF into account as it is negligible here.

### C.1.3. Assumptions

We take the following assumptions and approximation:

A 1 (Spatial separability). All sources are spatially separated on HST images.

A 2 (Spatial invariability). The spatial shape of a source does not vary with the wavelength

A 3 (Spectral invariability). Sources are modeled with a unique spectrum.

A 4 (FSF invariability). The transfer function HST-MUSE is known, spatially invariant, and can be approximated as constant over a given spectral band.

Assumption A1 is valid for almost all sources in the considered fields. When it is not valid, the proposed method cannot do further de-blending. Assumptions A2 and A3 are also approximations that allow the sources to be separated. Note however that some sources, like close-by galaxies that are spatially resolved on MUSE, have velocity fields that will deform the observed spectra, breaking the A3 hypothesis. Other sources have emissions at certain wavelengths that differ greatly from the continuum, or even have no continuum, making the HST information irrelevant in these cases. Assumption A4 is in practice quite accurate. The FSF of HST is now well known and MUSE FSF has also been well studied (Villeneuve et al. 2011; Carfantan 2014). The latter is modeled by a 2D circular Moffat 2D function. The Moffat scale parameter varies slowly with the wavelength such that FSF can be considered constant on tens or even hundreds of MUSE spectral bands.

Assumption A2 allows the model to be placed within the framework of a linear mixing model:

$$ \mathbf{X} \approx \mathbf{U} \times \mathbf{D}, \quad \mathbf{U} \in \mathbb{R}^{N \times k}, \quad \mathbf{D} \in \mathbb{R}^{k \times \lambda} $$

where $k$ is the number of sources in the data ($k \ll N$), $\mathbf{D} \in \mathbb{R}^{k \times \lambda}$ is the sources spectra matrix, and $\mathbf{U} \in \mathbb{R}^{n \times k}$ the intensity matrix. We note that here we are talking of intensities, and not abundances like in other hyperspectral contexts (remote sensing), so we cannot use the classical sum-to-one constraint on the rows of this matrix.

All of these assumptions allow us to deduce an intensity matrix **U** using HST images. From the MUSE data and this matrix **U**, we therefore try to directly estimate the matrix of the **D** spectra, without explicitly reconstructing **X**.

## C.2. Method

The strategy is as follows:

### Algorithm 1

1: **for** each (vectorized) image $\mathbf{Z}_i$ of the HST **do**

2: &nbsp;&nbsp;&nbsp;&nbsp;Get intensity matrix $\mathbf{U}_i$ from $\mathbf{Z}_i$.

3: &nbsp;&nbsp;&nbsp;&nbsp;**for** each spectral block $l$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Degrade $\mathbf{U}_i$ to MUSE spatial resolution to get $\tilde{\mathbf{U}}_{i,l}$

5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Reverse the system $\mathbf{Y}_b l = \tilde{\mathbf{U}}_{i,l} \mathbf{D}_{i,l}$ to find $\tilde{\mathbf{D}}_{i,l}$

6: &nbsp;&nbsp;&nbsp;&nbsp;Combine results: $\tilde{\mathbf{D}}_l = \sum_i a_{i,l} \tilde{\mathbf{D}}_{i,l}$, $1 \le l \le \lambda$ where $a_{i,l}$ is the spectral response of the filter $i$ at the wavelength $l$.

Steps of lines (4) and (5) are costly to compute on all 3600 wavelength MUSE. As FSF varies slowly along wavelengths, we work on blocks (200 or 300 spectral bands) where we estimate a mean FSF over the block. This considerably reduces the computational cost since it is sufficient to pseudo-invert a single matrix $\tilde{\mathbf{U}}_{i,l}$ for each block $l$. We note that we nevertheless obtain a full spectrum estimation, only the intensity matrix is assumed constant on the spectral block (A2). Along this process of linear operations we propagates the variance given with the MUSE data. Thus, we get at the end a variance for each estimated spectrum.

### C.2.1. Intensity matrix estimation

First step (line 2) is done using the segmentation map of HST sources. From this map we build a $N \times k$ matrix **H**, each column $j$ of **H** corresponds to the presence/absence binary mask of the
source j in the vectorized HST field. We then apply element-wise product between H and HST image Z_i to get the intensity matrix U_i:

$$U_i = H \circ Z_i. \qquad (C.3)$$

In step (4), we first perform a multiplication by the convolution transfer matrix between HST and MUSE. The convolution kernel from HST to MUSE K_HM is obtained from the convolution kernels of HST K_H and MUSE K_M, assumed to be known. Then the convolution result is resampled to the resolution of MUSE, using a linear interpolation.

We then have the intensity matrix at MUSE resolution $\tilde{U}_{i,l} = B_l U_i$.

Then for each object defined by the segmentation map, we get its own intensity map (even if the objects are not separable on the MUSE data).

### C.2.2. Spectra estimation

Step (c) could be done by simply solving a linear mixing model in the least squares sense. But when the intensity matrix is poorly conditioned, that is to say, when objects are spatially very close, solving the problem in the least squares sense may over-fit the noise and create spectral artifacts (see Sect. C.3).

In order to cope with this phenomenon, a classical approach is to add constraints to regularize the problem. There are a large number of possible regularizations in the literature, among which we can notably mention the penalization regularizations of the ridge type (Hoerl & Kennard 1970) or LASSO (Tibshirani 1996) and the regularizations by informational criterion like the Bayesian information criterion (BIC: Schwarz 1978). To answer to MUSE data specificities, we chose to exploit both ridge and sparse regularizations, as spectra are composed of a smooth spectral continuum and a small set of sharp emission/absorption lines.

We thus separately process lines and continuum. For this, the procedure is as follows: (i) obtaining the line cube $Y^r$ by subtracting the continuum cube estimated by a robust filtering (e.g., a median filter); (ii) reconstruction of the $\tilde{D}^r$ lines associated with the objects from the line cube obtained previously: we will see that a parsimonious regularization approach with model selection has been chosen; (iii) subtraction of the estimated contribution of the lines to obtain the remaining data $Y^c = Y - \tilde{U}\tilde{D}^r$; (iv) estimation of the spectral continuum $\tilde{D}^c$ using a ridge regularization on the data $Y^c$; and (v) combination of the two estimates $\tilde{D} = \tilde{D}^r + \tilde{D}^c$.

### Spectral line estimations

Sparse regularization of spectral lines avoids over-fitting the noise and favors line attribution to a minimal number of spectra, which is the most physically probable. Spectral lines are assumed to extend over several spectral slices. We can thus exploit this assumption to ensure a good lines reconstruction and avoid spectral discontinuities (where part of a line is associated with one object and the other part to another). Thus, we seek to do the model selection jointly over all the spectral support of the spectral line. To do so, it is thus necessary to estimate this spectral support. Strategy is thus as follows: (i) fast detection of all potential spectral supports of lines and (ii) for each segment, selection of objects with nonzero spectrum, then estimation of spectra using only selected objects.

The lines support detection is based on local extrema detection and is detailed in (Bacher 2017). It allows all the spectral supports of the potential lines in the studied area to be obtained quickly. We can then estimate the spectra of objects with regularization by selecting models on each of these supports.

The sparse regularization is done using BIC selection. BIC is written as $BIC = K \log(n) - 2 \log(\hat{L})$ where $\hat{L}$ is the model likelihood maximum, K is the number of free parameters, n is the number of samples (here pixels in the area). In our case we assume that the noise is Gaussian. On a spectral band, for a model M with k objects with nonzero spectrum, we have

$$BIC(M) = K \log(n) + \log(\hat{\sigma}_M^2), \qquad (C.4)$$

where $\hat{\sigma}_M^2$ is the empirical variance of the residuals obtained using the M model chosen. We have a nonzero spectrum free parameter plus one for the noise variance, and hence $K = k + 1$. To minimize the criterion BIC we see that there is a compromise to be found between a model that best explains the data ($\hat{\sigma}^2$ minimal) and a simplest possible model (k minimal). The selection of the best model is then reduced to a combinatorial problem where all the possible combinations of nonzero spectra are tested and the one that minimizes the BIC is retained. For computational reasons, instead of exploring all combinations we use a greedy variant for the BIC.

For each line l, we select the model $M_l$ over the spectral support of l using

$$M_l = \underset{M}{\operatorname{argmin}} \operatorname{BIC}(M), \qquad (C.5)$$

where $\operatorname{BIC}(M)$ comes from Eq. (C.4).

If we then note $\tilde{U}_M$ the intensity matrix composed only of the sources selected in M, we can compute the associated spectra $\tilde{D}_l^r$ for line l with

$$\tilde{D}_l^r = \underset{D}{\operatorname{argmin}} \|Y_l^r - \tilde{U}_{M_l} D\|_2^2. \qquad (C.6)$$

### Continuum regularization

After subtracting the estimated contribution of the lines, we now estimate the spectral continua of the objects using ridge regularization. For a given block of sheets l and an HST image i we note the column vector $y = \underbrace{Y_l}_{n \times 1}, \tilde{U} = \tilde{U}_{i,l}$ and $d = \underbrace{D_{i,l}}_{k \times 1}$. Ridge regularized least squares admits an analytical solution of the form

$$\tilde{d} = (\tilde{U}^T \tilde{U} + \alpha I_k)^{-1} \tilde{U}^T y, \qquad (C.7)$$

with $I_k$ the identity matrix of $\mathbb{R}^{k \times k}$ and where $\alpha \ge 0$ is the regularization parameter to be fixed.

The optimal value (in the Mean Square Error sense) of regularization parameter $\alpha$ is estimated using generalized cross validation (GCV; Golub et al. 1979). In order to increase the robustness of the cross-validation, only one regularization parameter is searched per block of a few tens of spectral sheets. This amounts to considering that the signal-to-noise ratio of the data (deprived of the lines) evolves weakly according to the wavelength. This hypothesis is valid if we consider that the sky's lines have been sufficiently well subtracted.
## Flux correction

A natural drawback of ridge regularization is that the estimated solution is biased toward zero, reflected in our case by a loss of flux on the estimated spectra, which is of course not desired. To overcome this drawback, we take advantage of the redundancy of penalized estimators along the spectrum: by looking for a multiplicative factor per spectral block of sufficiently large size (on the order of a hundred sheets), we can make sure that the average flow (on the block) is preserved, while preserving the benefit of regularization. This correction factor is obtained in the following way: one looks for a diagonal matrix $F = \text{diag}((f_j)_{1 \le j \le k})$ such as one has for continuum $Y^c = \tilde{U}F\tilde{D}^c$. The matrix $F$ being diagonal this amounts to defining

$$f_j = \frac{\left\langle \left( Y^c (\tilde{D}^c)^+ \right)_j, \tilde{U}_j \right\rangle}{\| \tilde{U}_j \|_2^2}, \quad \text{for } 1 \le j \le k. \qquad (C.8)$$

We note that if this least-squares correction was done wavelength by wavelength it would in fact cancel the ridge regularization.

The regularized approach of spectra estimation for an HST image $i$ within a FSF-constant spectral block $j$ is summarized in the algorithm 2. We note that the whole method stays unsupervised as regularizations are without external parameters (for lines estimation) or with auto-estimated parameter (for continuum estimation).

### Algorithm 2 Regularized procedure of spectra estimation

1: Input: data $Y \leftarrow Y_j$, intensity matrix $\tilde{U} \leftarrow \tilde{U}_{i,j}$

2: Estimation of $Y^c$ ▶ Continuum estimation (robust filtering)

3: Lines estimation $Y^r = Y - Y^c$ ▶ Continuum subtraction

4: Line detections

5: **for** line $l$ **do**

6: Computation of $M_l$ using eq. (C.5) ▶

7: Computation of $\tilde{D}_l^r$ using eq. (C.6)

8: Construction of $\tilde{D}^r$ ▶ Spectral concatenation

9: Computation of $Y^c = Y - Y^r$

10: **for** bloc spectral $b$ **do**

11: Computation of $\alpha_{m+1s}$ using GCV

12: Computation of $\tilde{D}^r$ ▶ Continuum estimation using ridge

13: Construction of $\tilde{D}^c$ ▶ Spectral concatenation

14: Computation of $F$ using (C.8) ▶ Flux correction factors

15: Computation of $\tilde{D} = \tilde{D}^r + F\tilde{D}^c$ ▶ Lines + Continuum

16: Output: $\tilde{D}$ ▶ Sources spectra (for image HST $i$ and FSF-constant spectral band $j$)

Finally, in order to take into account the residuals of the atmosphere spectrum subtraction, we add throughout the unmixing procedure, the estimation of a spectrum spatially constant in the considered area, corresponding to the sky background spectrum. This spectrum can be seen as the intercept of the studied regression problem.

### C.3. Validation on simulated data

#### C.3.1. Simulation settings

In order to evaluate the validity of this method, we built a series of data with sources that are increasingly close to each other spatially. To do this, we built a datacube (160×160) in flux units with

<figure><img src="image_45.png" /><figcaption>Fig. C.1. Simulation example. First row: "HST" segmentation map. Second row: "MUSE" white image. Distance of sources is characterized by the unit-less conditioning number $c$, increasing here from 1.2 (left column) to 3.7 (right column).</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>Fig. C.2. Ground truth spectra of the two simulated sources.</figcaption></figure>

high resolution (called HST data by abuse of language), which was then degraded following the MUSE model (same subsampling and core of convolution) to obtain data at low spatial resolution (which will be called MUSE data by abuse of language). An additive Gaussian noise is then applied, slightly correlated spatially by the application of a 3 by 3 spatial core.

We can see in Fig. C.1 the simulated cube white image at MUSE and HST resolution. Simulated spectra are composed of an emission line for each object and a spectral continuum (see fig. C.2).

De-blending difficulty is measured using the conditioning number $c$ of intensity matrix $\tilde{U}$ (at MUSE spatial resolution):

$$c = \frac{s_{\max}(\tilde{U})}{s_{\min}(\tilde{U})}, \qquad (C.9)$$

where $s_{\max}(\tilde{U})$ and $s_{\min}(\tilde{U})$ are the minimum and maximum singular values of the intensity matrix, and $c$ is a bound of estimation error of $\tilde{D}$ relative to a perturbation on data $Y$. The closer $c$ is to 1, the easier the inversion problem is. On the contrary, high values of $c$ will imply instability and over-fitting of the noise.

We then measured the fidelity to the ground truth (the mean of the intercorrelations between each estimated spectrum, $\tilde{d}_j$, and its ground truth, $d_j$),

$$f = \frac{1}{k} \sum_{j=1}^{k} \frac{\langle \tilde{d}_j, d_j \rangle}{\| \tilde{d}_j \| \cdot \| d_j \|},$$
<figure><img src="image_47.png" /><figcaption>Fig. C.3. Ground truth (green dashed-dot) versus estimated spectra (proposed method in orange, least squares in dashed blue) for a difficult case (condition number equals 6)</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>Fig. C.4. Zoom over a part of the non-regularized least squares estimation of both sources spectra, for a difficult case (condition number equals 6). Some anticorrelated artifacts are visible.</figcaption></figure>

and the inter-correlation,

$$ic = \frac{\langle \hat{d}_0, \hat{d}_1 \rangle}{\|\hat{d}_0\| \cdot \|\hat{d}_1\|},$$

between the two estimated spectra.

### C.3.2. Results

We compared the proposed regularization method to a simpler least squares approach that has a closed-form solution given by

$$\hat{D}_{i,l} = (\tilde{U}_{i,l}^T \tilde{U}_{i,l})^{-1} \tilde{U}_{i,l}^T Y_l. \quad (C.10)$$

Figure C.3 shows one of the two estimated spectra, around its spectral emission line, and when the sources are very close spatially ($c \approx 6$). We can see that both the proposed method and a non-regularized least squares approach estimate accurately the emission line but the proposed method is better at avoiding increased variance around the peak.

The results shown in Fig. C.5 describe the evolution of performance as the conditioning degrades. Performances of the proposed method are clearly more robust to ill-conditioning than a simple least squares approach. In particular, variance of the estimated spectra residuals increases and the estimated spectra become strongly negatively correlated. This anticorrelation is explained by the positivity of the $\tilde{U}$ intensity matrix. This phenomenon is also illustrated by Figs. C.3 and C.4, where the unmixed spectra when the sources are very close spatially ($c \approx 6$) are shown. The robustness to anticorrelated patterns in the proposed method allows it to keep a good fidelity to the sources spectra (fig C.5c).

<figure><img src="image_49.png" /><figcaption>Fig. C.5. Performance evolution function of conditioning. The dashed blue line indicates no regularization and the orange line with regularization. Results are averaged over ten Monte Carlo runs.</figcaption></figure>
## Appendix D: Emission and absorption lines

The lines input table used in `pyPlatefit` is given in Tab. D.1. We note that the default line table can also be retrieved directly from the `pyPlatefit` python interface<fnref n="29" />.

The lines are grouped into 3 different families (column Set): the Balmer series (bal), non-balmer emission lines (for) and ISM absorption lines (ism). The vacuum rest wavelength in Å is given in the column Wave. The multiplet (Mult) column is used to group lines together when performing the fit. Mult is an integer value following the adopted notation in the literature (e.g., [OII]3727 for the [O II]λλ3726,3729 doublet). The Main boolean column can be used to filter main and secondary lines. The two boolean columns, Emi and Abs, indicate if the line is in emission and absorption, respectively. The Res column is used to perform an independent fit for some major resonant lines. We note that some emission lines like OVI and NV, although resonant, are not fitted separately because they were too faint in our galaxy sample.

<table><caption>**Table D.1.** `pyPlatefit` lines definition.</caption><thead><tr><th>Name</th><th>Set</th><th>Wave</th><th>Mult</th><th>Main</th><th>Emi</th><th>Abs</th><th>Res</th></tr></thead><tbody><tr><td>OVI1032</td><td>for</td><td>1031.91</td><td>1033</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>OVI1038</td><td>for</td><td>1037.61</td><td>1033</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>LYALPHA</td><td>bal</td><td>1215.67</td><td>0</td><td>T</td><td>T</td><td>F</td><td>T</td></tr><tr><td>NV1238</td><td>for</td><td>1238.82</td><td>1240</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>NV1242</td><td>for</td><td>1242.80</td><td>1240</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>SiII1260</td><td>ism</td><td>1260.42</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>OII3702</td><td>ism</td><td>1302.17</td><td>1303</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>SiII1304</td><td>ism</td><td>1304.37</td><td>1303</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>CII1334</td><td>ism</td><td>1334.53</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>SiIV1394</td><td>for</td><td>1393.76</td><td>1403</td><td>F</td><td>T</td><td>T</td><td>F</td></tr><tr><td>OIV1397</td><td>for</td><td>1397.23</td><td>1403</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>OIV1400</td><td>for</td><td>1399.78</td><td>1403</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>SiIV1403</td><td>for</td><td>1402.77</td><td>1403</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>NIV1486</td><td>for</td><td>1486.50</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>SiII1527</td><td>ism</td><td>1526.71</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>CIV1548</td><td>for</td><td>1548.20</td><td>1549</td><td>T</td><td>T</td><td>T</td><td>T</td></tr><tr><td>CIV1550</td><td>for</td><td>1550.77</td><td>1549</td><td>T</td><td>T</td><td>T</td><td>T</td></tr><tr><td>FeII1608</td><td>ism</td><td>1608.45</td><td>1610</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>FeII1611</td><td>ism</td><td>1611.20</td><td>1610</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>HeII1640</td><td>for</td><td>1640.42</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>OIII1660</td><td>for</td><td>1660.81</td><td>1663</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>OIII1666</td><td>for</td><td>1666.15</td><td>1663</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>AlII1671</td><td>ism</td><td>1670.79</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>NIII1750</td><td>for</td><td>1749.67</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>AlIII1854</td><td>ism</td><td>1854.10</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>AlIII1862</td><td>ism</td><td>1862.17</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>SiIII1883</td><td>for</td><td>1882.71</td><td>1886</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>SiIII1892</td><td>for</td><td>1892.03</td><td>1886</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>CIII1907</td><td>for</td><td>1906.68</td><td>1909</td><td>T</td><td>T</td><td>F</td><td>F</td></tr><tr><td>CIII1909</td><td>for</td><td>1908.73</td><td>1909</td><td>T</td><td>T</td><td>F</td><td>F</td></tr><tr><td>CII2324</td><td>for</td><td>2325.40</td><td>2326</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>CII2325</td><td>for</td><td>2326.11</td><td>2326</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>CII2326</td><td>for</td><td>2327.64</td><td>2326</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>CII2328</td><td>for</td><td>2328.84</td><td>2326</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>FeII2344</td><td>ism</td><td>2344.21</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>FeII2374</td><td>ism</td><td>2374.46</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr></tbody></table>

<table><caption>**Table D.1. Continued.**</caption><thead><tr><th>Name</th><th>Set</th><th>Wave</th><th>Mult</th><th>Main</th><th>Emi</th><th>Abs</th><th>Res</th></tr></thead><tbody><tr><td>FeII2382</td><td>ism</td><td>2382.76</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>NeIV2422</td><td>for</td><td>2421.83</td><td>2424</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>NeIV2424</td><td>for</td><td>2424.42</td><td>2424</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>OII2470</td><td>for</td><td>2471.02</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>FeII2586</td><td>ism</td><td>2586.65</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>FeII2600</td><td>ism</td><td>2600.17</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>MgII2796</td><td>for</td><td>2796.35</td><td>2799</td><td>F</td><td>T</td><td>T</td><td>T</td></tr><tr><td>MgII2803</td><td>for</td><td>2803.53</td><td>2799</td><td>F</td><td>T</td><td>T</td><td>T</td></tr><tr><td>MgI2853</td><td>ism</td><td>2852.97</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>NeV3426</td><td>for</td><td>3426.85</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>OII3726</td><td>for</td><td>3727.09</td><td>3727</td><td>T</td><td>T</td><td>F</td><td>F</td></tr><tr><td>OII3729</td><td>for</td><td>3729.88</td><td>3727</td><td>T</td><td>T</td><td>F</td><td>F</td></tr><tr><td>H11</td><td>bal</td><td>3771.70</td><td>0</td><td>F</td><td>T</td><td>T</td><td>F</td></tr><tr><td>H10</td><td>bal</td><td>3798.98</td><td>0</td><td>F</td><td>T</td><td>T</td><td>F</td></tr><tr><td>H9</td><td>bal</td><td>3836.47</td><td>0</td><td>F</td><td>T</td><td>T</td><td>F</td></tr><tr><td>NeIII3869</td><td>for</td><td>3870.16</td><td>0</td><td>T</td><td>T</td><td>F</td><td>F</td></tr><tr><td>HeI3889</td><td>for</td><td>3889.73</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>H8</td><td>bal</td><td>3890.15</td><td>0</td><td>F</td><td>T</td><td>T</td><td>F</td></tr><tr><td>CaK</td><td>ism</td><td>3933.66</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>CaH</td><td>ism</td><td>3968.45</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>NeIII3967</td><td>for</td><td>3968.91</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>HEPSILON</td><td>bal</td><td>3971.20</td><td>0</td><td>F</td><td>T</td><td>T</td><td>F</td></tr><tr><td>HDELTA</td><td>bal</td><td>4102.89</td><td>0</td><td>T</td><td>T</td><td>T</td><td>F</td></tr><tr><td>CaG</td><td>ism</td><td>4304.57</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>HGAMMA</td><td>bal</td><td>4341.68</td><td>0</td><td>T</td><td>T</td><td>T</td><td>F</td></tr><tr><td>OIII4363</td><td>for</td><td>4364.44</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>HBETA</td><td>bal</td><td>4862.68</td><td>0</td><td>T</td><td>T</td><td>T</td><td>F</td></tr><tr><td>OIII4959</td><td>for</td><td>4960.30</td><td>0</td><td>T</td><td>T</td><td>F</td><td>F</td></tr><tr><td>OIII5007</td><td>for</td><td>5008.24</td><td>0</td><td>T</td><td>T</td><td>F</td><td>F</td></tr><tr><td>MgB</td><td>ism</td><td>5175.44</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>HeI5876</td><td>for</td><td>5877.25</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>NaD</td><td>ism</td><td>5891.94</td><td>0</td><td>F</td><td>F</td><td>T</td><td>F</td></tr><tr><td>OI6300</td><td>for</td><td>6302.05</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>NII6548</td><td>for</td><td>6549.85</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr><tr><td>HALPHA</td><td>bal</td><td>6564.61</td><td>0</td><td>T</td><td>T</td><td>T</td><td>F</td></tr><tr><td>NII6584</td><td>for</td><td>6585.28</td><td>0</td><td>T</td><td>T</td><td>F</td><td>F</td></tr><tr><td>SII6717</td><td>for</td><td>6718.29</td><td>0</td><td>T</td><td>T</td><td>F</td><td>F</td></tr><tr><td>SII6731</td><td>for</td><td>6732.67</td><td>0</td><td>T</td><td>T</td><td>F</td><td>F</td></tr><tr><td>ARIII7135</td><td>for</td><td>7137.80</td><td>0</td><td>F</td><td>T</td><td>F</td><td>F</td></tr></tbody></table>

<footnote n="29">²⁹ `pyPlatefit.readthedocs.io/en/latest/tutorial.html#Selecting-and-updating-emission/absorption-lines</footnote>