© The Authors 2026

1. Introduction

The evolution of a galaxy is imprinted in its stellar structures. By detecting and characterising these structures with increasing accuracy, we can make increasingly detailed inferences about a galaxy’s history. Integral field spectroscopy (IFS; Bacon et al. 1995) allows us to map out stellar kinematics and populations in extragalactic systems. The growing size and quality of IFS galaxy surveys (Sánchez et al. 2012; Bundy et al. 2015; Scott et al. 2018; Sarzi et al. 2018; Emsellem et al. 2022; van de Sande et al. 2024) demand increasingly sophisticated methods for measuring these stellar properties optimally.

Full-spectrum fitting (e.g. Cappellari & Emsellem 2004; Tojeiro et al. 2007; Koleva et al. 2009; Cid Fernandes 2018) is the prevailing method for measuring stellar kinematics and populations. In this approach, the spectrum is modelled directly in pixel space using a wider wavelength range than in earlier methods which focussed on a few equivalent widths or spectral indices (Worthey 1994; Kauffmann et al. 2003). Full-spectrum fitting typically combines single stellar population templates (e.g. Bruzual & Charlot 2003; Vazdekis et al. 2010; Maraston et al. 2020) with line-of-sight velocity distributions (LOSVDs) that are either described via Gauss Hermite expansions (Rix & White 1992; Saha & Williams 1994) or non-parametrically (Rix & White 1992; Saha & Williams 1994). These developments have cemented the position of full-spectrum fitting as a flexible and important tool for characterising stellar systems. However, one drawback of the method is that it treats the 3D datacube as a number of 1D spectra that are modelled independently.

Three-dimensional methods are IFS analyses in which spatial and spectral dimensions are modelled simultaneously. This feature allows us to leverage spatial regularities during reconstructions to provide more accurate measurements. 3D modelling is an established technique for modelling gas components, and several tools are available for modelling gas disks (Józsa et al. 2007; Bouché et al. 2015; Di Teodoro & Fraternali 2015; Varidel et al. 2019) alongside specialised tools for handling more complex structures such as bars (Davis et al. 2013) and gravitationally lensed observations (Rizzo et al. 2018). The relative dearth of 3D methods for stellar reconstructions may be explained by the fact that this is a much more computationally demanding problem, because while gas components are characterised by single emission lines, information on the stellar content is spread throughout the full spectrum.

The projected Nesterov Kaczmarz reconstruction method (Hinterer et al. 2023) was the first algorithm that was developed for a 3D full-spectrum fitting. PNKR leverages the mathematical structure of the forward model to efficiently reconstruct stellar properties from IFS data. Additionally, PNKR models the joint distribution over populations and kinematics, in contrast to most full-spectrum fitting codes, which assume that at a given position in the galaxy, kinematics and populations are independent. This means that PNKR has the freedom to infer local relations between populations and kinematics. Previously, these types of measurements have been made via decomposition techniques for galaxies with spectroscopically resolved substructure, for example, galaxies with large-scale stellar counter-rotation (Coccato et al. 2011; Johnston et al. 2013; Coccato et al. 2013; Pizzella et al. 2018). We test PNKR using mock data experiments, and in particular, we explore its ability to make joint population-kinematic inferences.

We also introduce a second approach for a 3D full-spectrum fitting based on Bayes-LOSVD (Falcón-Barroso & Martig 2021), which is one of a number of recent tools (Mehrgan et al. 2023; Gasymov & Katkov 2024) for non-parametric LOSVD recoveries. The use of non-parametric LOSVDs has enabled the detection of weak kinematic signatures (Mehrgan et al. 2023) and been shown to significantly affect dynamical inferences. For instance, Reiter et al. (2025) applied orbit-superposition modelling to the archetypal counter-rotating galaxy NGC4550 (Rubin et al. 1992) and found that non-parametric kinematics results in significantly colder disks and lower masses than predicted when using Gauss Hermite LOSVDs. Despite their benefits, the additional flexibility of non-parametric methods comes at the cost of greater degeneracy. Typically, this is handled using regularisation in velocity space. Fig. 7 of Falcón-Barroso & Martig (2021) for example illustrates the effect of velocity regularisation on LOSVDs recovered from NGC4550. The bottom two rows of that figure show that velocity regularisation tends to oversmooth LOSVD bimodality, and the authors therefore recommend non-regularised fits over regularised ones. We forego regularisation in velocity space and instead explore the effect of spatial regularisation on LOSVD recoveries.

The main aim of this work is to assess the effect of 3D methods for stellar recoveries. Section 2 introduces the recovery algorithms by first describing the existing PNKR algorithm and then introducing a novel implementation of spatial regularisation within the Bayes-LOSVD framework. Section 3 describes the generation of mock data, and then, Sect. 4 presents the results of reconstructions on mock data and compares the performance of PNKR and Bayes-LOSVD. We discuss the implications and limitations of these results in Sect. 5. Finally, in Sect. 6, we summarise our main conclusions and discuss the current status and future plans for the public release of our recovery codes.

2. Methods

This section introduces methods for reconstructing stellar population and kinematic distributions from IFS datacubes. Concretely, we aim to recover the joint distribution f(x, v, z, t) of stars over 2D position x, line-of-sight velocity v, stellar metallicity z, and stellar age t. The forward model introduced by Hinterer et al. (2023) to relate f to the 3D datacube signal $\overline{y}(\mathit{x},\lambda )$ at position x and wavelength λ is given by

$$\begin{array}{c}\hfill \overline{y}(\mathit{x},\lambda )={\displaystyle \int \int \int \frac{1}{1+v/c}}S(\frac{\lambda }{1+v/c},z,t)f(\mathit{x},v,z,t)\mathrm{d}t\mathrm{d}z\mathrm{d}v,\end{array}$$(1)

where the kernels S(λ; z, t) are single stellar population (SSP) templates. Equation (1) represents an idealised model assuming that effects unrelated to the stellar light (e.g. dust absorption and gas emission) can be otherwise modelled or corrected for.

We investigate two reconstruction algorithms, Bayes-LOSVD and PNKR. While the former is limited to LOSVD recovery, PNKR models the full joint density f(x, v, z, t) over position, velocity, and populations. Both algorithms have 1D versions, where spectra at different locations are fit independently (PNKR-1D and BLOSVD-1D), as well as versions that fit multiple spectra simultaneously by accounting for spatial correlations (PNKR-3D and BLOSVD-2D). The following sections introduce the two algorithms and demonstrate their 1D versions on a toy problem of recovering a Gaussian ground-truth LOSVD. For this toy problem, we used a single template spectrum consisting of a single absorption line instead of SSP templates, and we used mock data with a signal-to-noise ratio (S/N) of 50. Applications to a more realistic scenario are given in Sect. 4.

2.1. PNKR

The projected Nesterov Kaczmarz Reconstruction method (PNKR, Hinterer et al. 2023, henceforth H23) was the first algorithm developed for 3D full-spectrum fitting. A full mathematical description of the algorithm is given in H23; the remainder of this section provides a qualitative overview of the existing algorithm and introduces one minor modification.

Given an observed datacube and set of SSP templates, PNKR inverts Eq. (1) to return an estimate of the full joint density f(x, v, z, t), which is described non-parametrically as a 5D array of weights. PNKR is an iterative regularisation method. This means that the density f is iteratively updated, with later iterations adding increasingly detailed refinements. By terminating the iterations at a suitable time, the solution should be regularised so as to avoid spurious details that are driven by over-fitting to noise.

The PNKR iteration scheme proceeds as follows. The density is initialised as f₀ = 0, and all wavelengths λ in the observed datacube are marked as active. In later iterations, PNKR cycles through all active wavelengths λ_k and applies an update to the density f_i, which is a linear function of the datacube slice y(x, λ_k) and a pre-defined step size. If the reconstructed slice y(x, λ_k) fits the observed datacube within a provided S/N-dependent noise tolerance, the wavelength λ_k is marked as inactive and is then ignored in later iterations. Iterations are continued until one of two stopping criteria is achieved: either the number of active wavelengths reaches zero, or it has plateaued within some tolerance. Figure 1 illustrates the PNKR iterative regularisation scheme on the toy problem.

Fig. 1. Toy demonstration of LOSVD recovery with PNKR. Left: Number of active wavelengths per iteration of the algorithm, with iterations 300 and 1000 highlighted. Right: Corresponding LOSVD recoveries. At iteration 300, the recovery is smooth and matches the truth well, while extending to iteration 1000 introduces noise.

The name projected Nesterov Kaczmarz reconstruction describes the following aspects of the algorithm. After each iteration, the solution is projected to non-negative values to satisfy the physical constraint that f ≥ 0. The algorithm employs Nesterov acceleration (Nesterov 1983) to speed up convergence. Lastly, PNKR is an example of a Kaczmarz method Engl et al. (1996), Kaltenbacher et al. (2008), Kindermann & Leitao (2014). This class of methods iteratively solves large linear problems with a specific system structure. PNKR is a Kaczmarz method by virtue of the fact that we cycle through wavelengths at each iteration, which has the benefit that each update comprises a relatively cheap computation.

PNKR can be run in two modes: PNKR-1D, which is equivalent to independent spectral fits in each spaxel, and PNKR-3D, which accounts for spatial correlations during the reconstructions. PNKR-1D follows the iterative scheme described above, which imposes no explicit coupling between adjacent spaxels in the recovery. PNKR-3D incorporates spatial correlations in the theoretical framework of adjoint Sobolev embedding operators (Hubmer et al. 2023). For our problem, this means that we filter the density f with a spatial smoothing kernel between iterations. The smoothing kernel K(x) is parametrised to allow us to explore different degrees of spatial smoothing. It takes a shape parameter α and a dimensionless scale length β. Roughly, these can be interpreted as follows: when α = 1(2), the kernel is Laplacian (Gaussian), and when β = 1, the kernel scale is comparable to the full spatial domain. In practice, the kernel is defined in Fourier space via

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathit{x}}\{K(\mathit{x})\}(k)={(1+\frac{\beta }{2}{\left|\left|\frac{k}{k_{\max }}\right|\right|}_2^2)}^{-\alpha },\end{array}$$(2)

where k_max = 1/Δx and Δx is the spaxel width. The convolution of f with K is then numerically performed in Fourier space using fast Fourier transforms.

Lastly, we note that we have made one modification to the original H23 algorithm in order to account for a spatially varying S/N. Specifically, we provided the algorithm with the true noise level used to generate our mock data, that is, σ(x, λ) from Eq. (6), which PNKR uses as a weighting factor when computing the overall residual between reconstructed and observed datacube slices. For real IFS data, standard data reduction pipelines (e.g. Weilbacher et al. 2020) provide an estimate of this noise level.

2.2. Bayes-LOSVD

Bayes-LOSVD (Falcón-Barroso & Martig 2021, henceforth FBM21) is a framework for probabilistic recoveries of non-parametric LOSVDs from 1D spectral fits. This section describes the existing 1D algorithm before we introduce our novel extension, which incorporates spatial regularisation.

Bayes-LOSVD imposes two assumptions on Eq. (1) to reduce the problem size. Firstly, in common with most full spectral fitting codes, FBM21 assume that at a given position x, all stellar populations share the same kinematics. This is equivalent to factorising the 5D density f as

$$\begin{array}{c}\hfill f(\mathit{x},v,z,t)=f_1(\mathit{x})f_2(v|\mathit{x})f_3(z,t|\mathit{x}),\end{array}$$(3)

which says that, locally (i.e. at a given x), the LOSVD (f₂) and stellar population distribution (f₃) are independent. Secondly, FBM21 replace SSP templates, of which there are typically several hundred, with a few templates distilled from the SSP library via principal component analysis (PCA). These two modifications reduce the size of the problem so that it is amenable to posterior sampling via Markov chain Monte Carlo (MCMC) algorithms. From a 1D spectral fit at position x, the main outputs of Bayes-LOSVD are posterior summaries (median and credible intervals) of the LOSVD f₂(v|x), which is described non-parametrically as an array of weights in evenly spaced velocity bins.

Bayes-LOSVD provides a number of different options for regularisation of the LOSVD in velocity space. The need for regularisation is illustrated in Fig. 2, which shows un-regularised LOSVD recoveries on the toy problem. While the LOSVD is well recovered when the template spectrum has a one-pixel wide absorption line (left panel), when the line width is two pixels (right panel), noise is amplified, and the credible intervals grow significantly. At resolutions of interest for galactic astronomy, absorption lines in SSP templates are typically resolved into several pixels, hence un-regularised LOSVD recoveries are prone to exhibit the spike artefacts seen in the right panel of Fig. 2. Regularisation in velocity space offers one way to manage this, but it can oversmooth genuine kinematic substructures, as discussed in Sect. 1. We therefore did not use any velocity regularisation in Bayes-LOSVD, but instead explored spatial regularisation as an alternative approach.

Fig. 2. Demonstration of LOSVD recovery with BLOSVD-1D. Each panel shows BLOSVD-1D recoveries via their posterior median (thick blue line), 99% credible interval (shaded region), and ten randomly selected posterior samples (thin blue lines). For the left panel, the template spectrum contains a single absorption line exactly one pixel wide; here, the recovery tightly encloses the truth. For the right panel, the absorption line is two pixels wide, which introduces a degeneracy between neighbouring velocity bins, evidenced by the zig-zagging posterior samples. This degeneracy results in noisy median recoveries and inflated credible intervals.

Our novel implementation of spatial regularisation in Bayes-LOSVD uses a conditional autoregressive (CAR) prior (Besag 1974). CAR priors are used in Bayesian hierarchical models of spatial data (Morris et al. 2019). They encode the prior belief that adjacent sites in the model are normally distributed with respect to one another with some degree of correlation. Mathematically, this is related to the spatial smoothing in PNKR for a Gaussian kernel. For our application, we used a CAR prior to promote spatial smoothness of the LOSVD. Specifically, when we denote the LOSVD value in velocity bin i at spatial bin j as y_ij := f₂(v_i|x_j) and let 〈y_i〉_Nbrs(j) be the mean LOSVD value in all spatial bins neighbouring j, then the CAR prior encodes that

$$\begin{array}{c}\hfill p(y_{\mathit{ij}}|{\langle y_i\rangle }_{\mathrm{Nbrs}(j)},{\rho }_{\mathrm{CAR}},{\sigma }_{\mathrm{CAR}})=\mathcal{N}({\rho }_{\mathrm{CAR}}{\langle y_i\rangle }_{\mathrm{Nbrs}(j)},{\sigma }_{\mathrm{CAR}}^2),\end{array}$$(4)

where the two hyper-parameters are a correlation 0 < ρ_CAR < 1 and a scale σ_CAR > 0. The full CAR prior extends Eq. (4) to the joint distribution for all spatial bins using an efficient sparse matrix implementation to encode adjacency relations¹. We used CAR priors to spatially smooth the LOSVD in each velocity bin i, but used the same hyper-parameters (ρ_CAR, σ_CAR) across all i for simplicity.

To keep problem sizes manageable for these initial experiments, we did not use the CAR prior across the full cube, but instead restricted its use to columns of spaxels. We therefore denote this method BLOSVD-2D, rather than 3D. The choice of restricting to columns was made for convenience, but is not a fundamental limitation of this method. It would be straightforward, for example, to replace columns with clusters of adjacent spaxels.

We newly re-implemented Bayes-LOSVD with modes for 1D fitting (BLOSVD-1D) and 2D fitting (BLOSVD-2D) using the probabilistic programming language NumPyro (Phan et al. 2019). This provided speed-ups of a factor of ≈10 compared to the original Stan (Carpenter et al. 2017) implementation used by FBM21. In keeping with FBM21, we used the no-U-turn sampler (NUTS; Hoffman & Gelman 2014) algorithm for sampling. This is a gradient-based MCMC algorithm that employs a warm-up stage for calibration, followed by a sampling stage during which posterior samples are efficiently drawn. We checked convergence via the standard NUTS diagnostics (Stan Development Team 2024): the potential scale reduction factor $\widehat{R}$, effective sample size N_eff, and number of divergent MCMC transitions N_div. Unless explicitly mentioned, all Bayes-LOSVD results shown henceforth are converged subject to the criteria that all parameters satisfy $|\widehat{R}-1|<0.03$, N_eff > 50 and N_div = 0.

3. Mock data

This section outlines the generation of mock datasets we used to test our reconstruction algorithms. In brief, we used simple analytic models to describe a ground-truth model f(x, v, z, t), evaluate the resulting datacube y(x, λ) via Eq. (1), and we then added noise to achieve a desired S/N. We did this all using the software popkinmocks (Jethwa 2023)².

Figure 3 illustrates the default ground-truth model used throughout this work. The model consists of two counter-rotating disk components: a young metal-rich thin disk (blue contours), and an older metal-poor thick disk (orange contours). The superposition of these two components gives rise to a mean velocity map (right panel) that exhibits an inversion, with negative velocities along the major axis swapping to positive velocities in the upper and lower left corners. The central and right panels of Fig. 3 illustrate the field of view and spatial sampling used for our mock observations: we simulated one side of the galaxy, sampled with a 30 × 30 grid of spaxels.

Fig. 3. Our ground-truth galaxy model in age-metallicity space (left) as a light-weighted image (centre) and mean velocity map (right). The model consist of two components. The dominant component (blue contours) is a young metal-rich thin disk with negative LOS velocity. The counter-rotating weaker component (orange contours) is older, more metal poor, and more extended vertically. Successive contours in the central panel show changes in flux by 15%. The grid of crosses in the right panel indicates the spaxels used for illustration in Figs. 4, 7, and 10.

Given this ground-truth model f(x, v, z, t), we next evaluated the noise-free datacube $\overline{y}(\mathit{x},\lambda )$ according to Eq. (1). To do this, we adopted MILES SSPs (Vazdekis et al. 2010) with BaSTI isochrones for base α-content (Pietrinferni et al. 2004) and a Chabrier IMF (Chabrier 2003), using 12 metallicity bins spanning −2.51 ≤ [M/H]≤0.47 and 38 age bins spanning 0.45–14.25 Gyr, and a wavelength range of 4800–5700 Å. For the velocity sampling, we adopted a range of [−750, 750] km/s, which fully covers the ground-truth model, divided into 41 bins with a width of ≈37 km/s. The popkinmocks software re-samples the SSPs in the equivalent log λ, and evaluates Eq. (1) using Fast Fourier Transforms to perform convolution of SSPs with LOSVDs. This gives the noise-free datacube $\overline{y}(\mathit{x},\lambda )$.

We generated five realisations of observed datacubes at varying S/N levels that we analyse independently throughout this work. To arrive at an observed datacube y(x, λ), we added noise ϵ(x, λ),

$$\begin{array}{c}\hfill y(\mathit{x},\lambda )=\overline{y}(\mathit{x},\lambda )+ϵ(\mathit{x},\lambda ),\end{array}$$(5)

which we sampled from a normal distribution with a mean 0 and a variance σ²(x, λ), that is,

$$\begin{array}{c}\hfill ϵ(\mathit{x},\lambda )\sim \mathcal{N}(0,{\sigma }^2(\mathit{x},\lambda )).\end{array}$$(6)

We took σ²(x, λ) proportional to $\overline{y}(\mathit{x},\lambda )$ to mimic shot noise and selected a constant of proportionality in order to achieve a desired S/N in the brightest spaxel. Throughout this work, we investigate five S/N values logarithmically spaced between S/N = 20 and 200. These S/N values refer to the maximum S/N in the cube. The noise model described above gives rise to spatial variations in S/N of around 30%. We did not apply any line-spread function beyond that inherent to the SSP templates, nor did we apply any point-spread function or other sources of noise.

4. Results

This section presents reconstructions using the mock datasets described in Sect. 3. Sections 4.1 and 4.2 focus purely on LOSVD reconstructions, first comparing the 1D versions of PNKR and BLOSVD (Sect. 4.1), and then, their 2D/3D equivalents (Sect. 4.2). In Sect. 4.3 we explore the ability of PNKR to recover the joint distribution over stellar populations and kinematics. Additional details of the implementation and fit convergence are given in Appendix A, and the spectral fits are shown in Appendix B.

Before presenting results, we quantify the size of the reconstruction problems we describe. We aim to reconstruct an unknown density f(x, v, z, t) with sampling in spatial, velocity, and metallicity dimensions as described in Sect. 3. In order to keep the run-times for PNKR within a reasonable range (≲ 48 hours), we chose to thin the age sampling of the SSP grid for the reconstructions by a factor of 3, which reduced the number of ages from 38 to 13. Nevertheless, in total, we sampled f(x, v, z, t) with 30 ⋅ 30 ⋅ 41 ⋅ 12 ⋅ 13 ≈ 6 × 10⁶ unknown entries. The simplifying assumptions of Bayes-LOSVD reduced this number significantly, but we stress that PNKR simultaneously solves for 6 × 10⁶ unknowns, which constitutes a significantly high dimensional problem.

4.1. LOSVD recovery from 1D fits

Figure 4 shows the recovery of light-weighted LOSVDs (defined in Appendix C.1) at 9 spaxels whose locations are shown by the crosses in Fig. 3. In spaxels 3–9, the two-component galaxy model gives rise to bimodal LOSVDs, while in spaxels 1 and 2, the components overlap in v_LOS, leading to unimodal LOSVDs. At S/N 200 (left grid), PNKR-1D recoveries (red lines) are good everywhere, showing only minor deviations from the truth, for instance at v_LOS = 0 in spaxel 3. The PNKR-1D recoveries are also smooth, suggesting that the iterative regularisation scheme of the algorithm has worked well. In contrast, the BLOSVD-1D recoveries (blue lines, showing posterior medians) at S/N = 200 show large spike artefacts similar to those seen in Fig. 2. At S/N 20 (right grid), both sets of recoveries are worse. Both algorithms produce erroneous flat wings extending to the edge of the velocity range. Additionally, PNKR-1D recoveries are over-smoothed in spaxels 5–9, where, in contrast, BLOSVD-1D does recover the strong dip between the two components. At this low S/N, BLOSVD-1D no longer produces spike artefacts, but the recovered LOSVD shapes are visibly worse than at high S/N.

Fig. 4. Comparison of LOSVDs recovered from 1D spectral fits. True LOSVDs (black) are shown alongside recoveries from PNKR-1D (red) and BLOSVD-1D (blue). The 3 × 3 grids in the sub-panels correspond to the spaxels highlighted in Fig. 3. At S/N 200 (left grid), PNKR recovers the LOSVD well everywhere, while BLOSVD recoveries show spike artefacts in some cases e.g. panels 1, 2 and 9. At S/N 20 (right grid) the recoveries are worse. Both algorithms produce LOSVDs with extended, flat wings which reach the edge of the velocity range. For spaxels 5-9 we see that BLOSVD-1D captures the dip between the two components whereas PNKR-1D oversmooths this feature.

Figure 5 summarises the LOSVD recovery capabilities of the 1D algorithms. It shows the LOSVD recovery error relative to the ground truth averaged over all spaxels and velocity bins (fully defined in Appendix C.2) as a function of S/N. When we compare the two methods, PNKR-1D (red line) achieves the lowest error for S/N > 40, while BLOSVD-1D (blue line) performs best at S/N = 20. These results show a smooth transition between the two S/N regimes depicted in Fig. 4, where oversmoothing affects the PNKR-1D recoveries at low S/N, and spike artefacts affect BLOSVD-1D at high S/N. Visually inspecting LOSVD recoveries at the intervening S/N, we find that oversmoothing appears for PNKR-1D below S/N ≈ 40, while spikes affect BLOSVD-1D above S/N ≈ 100. We reiterate here that we used BLOSVD-1D with no regularisation in velocity space. Including this would alleviate LOSVD spikes at the cost of oversmoothing the LOSVD bimodality, as discussed in Sect. 2.2.

Fig. 5. Error in recovered LOSVDs from 1D spectral fits as a function of S/N. Comparing the two methods, PNKR (red) achieves the smallest error at high S/N, while BLOSVD (blue) is best at low S/N. By modifying PNKR to retain all wavelengths throughout the recovery (pink line), we see an improvement at low S/N, with a recovery error similar to BLOSVD, at the expense of larger errors at high S/N.

We next quantify the prominence of the flat wing artefacts (seen in the right panel of Fig. 4) at different S/N values. We measured the wing prominence as the median value of recovered LOSVDs for |v_LOS|> 600 km/s averaged over all spaxels. We found that the prominence of wing artefacts is very similar for BLOSVD-1D and PNKR-1D. For both algorithms, the wing prominence decreases with S/N in a manner that is well described by a power law, with a wing prominence ∝S/N^−0.85 (shown in Appendix D). This power-law dependence of the wing prominence with S/N effectively suppresses the wing artefacts to a level where they are no longer visible on the plotted scale for S/N = 200 (left panel of Fig. 4).

We next investigated why PNKR-1D tends to oversmooth at low S/N. We varied the algorithm step-size and tolerance settings (see Appendix A.3 for details), but found that this had no effect. We then modified the algorithm itself. In the standard algorithm as described in Sect. 2.1, when a wavelength is well fit, it is deactivated and not used to update the density f in later iterations. The pink line in Fig. 5 shows LOSVD recovery errors for a modified algorithm without deactivation. This modification improves PNKR-1D reconstructions at S/N = 20, bringing errors in line with BLOSVD-1D, but for S/N > 60, the reconstructions become worse. This test suggests that while wavelength deactivation benefits the iterative regularisation of PNKR at high S/N, it is likely connected to oversmoothing at low S/N.

4.2. LOSVD recovery from 2D/3D fits

We next ran reconstructions accounting for spatial correlations between spaxels. We ran 16 PNKR-3D recoveries over a (5 × 5) grid of smoothing parameters (α, β) detailed in Appendix A.2. For BLOSVD-2D, it is possible in principle to infer the values of the smoothing parameters (ρ_CAR, σ_CAR) by assigning them suitable hyper-priors and then sampling their posteriors alongside the LOSVD. In practice, we found that this approach works well for the correlation parameter ρ_CAR, but leads to problems when applied to the scale parameter σ_CAR. Specifically, MCMC chains exhibit many divergent transitions when they explore low values of σ_CAR. For BLOSVD-2D, we therefore ran two reconstructions using two fixed values of σ_CAR = 0.001, 0.03, but retained the hyper-prior over ρ_CAR. All fits converged with σ_CAR = 0.03, but only the S/N = 200 cube converged with σ_CAR = 0.001, achieving a smaller recovery error with this lower value of σ_CAR (further details are given in Appendix A.4). Henceforth, in this section, all results for S/N < 200 are shown for σ_CAR = 0.03, while S/N = 200 results use σ_CAR = 0.001 unless otherwise stated.

Figure 6 shows LOSVD recovery errors for 2D/3D fits (dashed lines) compared to their 1D equivalents (solid lines, repeated from Fig. 5). PNKR-3D produces very little improvement (≈0.01%) compared to PNKR-1D at all S/N. This result is not limited by the range of smoothing parameters we explored (see Appendix A.2). Furthermore, we verified that this result does not change when (i) the algorithm step-size and tolerance settings were varied, (ii) we initialised PNKR-3D fits not from zero, but instead from the outputs of the PNKR-1D fits, and (iii) we used a smoothing filter based on wavelets and not on Fourier convolution. We conclude that PNKR-3D leads to very little improvement over PNKR-1D, but we also note that both versions of PNKR achieve lower reconstructions errors than BLOSVD for S/N > 40. We discuss this result further in Sect. 5.

Fig. 6. Error in recovered LOSVDs from 2D/3D fits as a function of S/N. The solid lines are repeated from Fig. 5, showing the results of 1D spectral fits with PNKR (red) and BLOSVD (blue). The dashed lines show results when fitting multiple spaxels simultaneously and promoting spatial smoothness in the recoveries. PNKR-3D (dashed red), which fits the full 3D datacube, produces results that are almost indistinguishable from the 1D PNKR fits. In contrast, BLOSVD-2D (dashed green), which fits columns of spaxels simultaneously, achieves a significant improvement over the 1D case for all S/N values.

For BLOSVD, the effect of accounting for spatial correlations in LOSVD reconstructions is very different. Fig. 6 compares 1D results (blue) with the 2D results (green). BLOSVD-2D significantly outperforms BLOSVD-1D at all S/Ns, with an improvement of 0.3% at S/N = 20 and 0.1% at S/N = 200. The remainder of this section focuses on a more detailed comparison of BLOSVD-1D and BLOSVD-2D.

Figure 7 compares median LOSVD reconstructions for BLOSVD-1D (blue) and 2D (green) for 9 example spaxels. At S/N = 200 (left grid), the most severe spike artefacts are effectively regularised away when using BLOSVD-2D. This is most evident for spikes occurring at the LOSVD maximum (e.g. spaxels 1, 2, and 5), whereas spikes in the LOSVD of the subdominant component with v_LOS > 0 show little improvement (e.g. spaxels 2, 3, and 4). At S/N = 20 (right grid), BLOSVD-2D recoveries show slight improvements in a few spaxels (e.g 4, 6, and 9), again more so for the v_LOS < 0 component than for the v_LOS > 0 component. The flat wing artefacts in the S/N = 20 recoveries are unchanged between 1D and 2D modelling (see Appendix D for a full characterisation of the wing artefacts as a function of S/N).

Fig. 7. Comparison of LOSVDs recovered from with BLOSVD-1D (blue lines) against BLOSVD-2D (green lines). At S/N 200 (left grid), the spike artefacts seen in the 1D fits are alleviated in the 2D fits (e.g. panels 1, 2, 5, and 9). At S/N 20 (right grid), the 2D fits improve the recovered LOSVD shapes in many cases, particularly for the dominant component with vLOS < 0 (e.g. panels 4, 6, and 9). At high and low S/N, the effect of 2D fitting is weaker for the subdominant component with vLOS > 0.

Figure 8 compares BLOSVD-1D and 2D reconstructions at S/N = 200 via channel maps. These show the spatial distribution of stars in a given velocity bin and provide an intuitive way to visualise the effect of spatial regularisation on LOSVD reconstructions. In most cases, the channel maps from 1D reconstructions (second column) are noisier than their 2D equivalents (third/fourth columns for σ_CAR = 0.03/0.001). For the negative velocity channel (top row), corresponding to the dominant galactic component, BLOSVD-2D provides some smoothing for either choice of σ_CAR, with σ_CAR = 0.001 (top right panel) being the most accurate. For the positive velocity channel (bottom row), the BLOSVD-2D reconstruction with σ_CAR = 0.03 provides little change with respect to the 1D case, while for σ_CAR = 0.001, the reconstructions are smoother, but remain clumpier than the truth. These results suggest that the effectiveness of spatial regularisation for a given velocity channel depends on σ_CAR and on the LOSVD density values in that channel.

Fig. 8. Channel maps showing the spatial distribution of stars in velocity bins labelled on the left. The first column shows true distributions, which are smooth by construction. The remaining columns show reconstructions for S/N = 200 data. The second column shows BLOSVD-1D reconstructions, which show large spaxel-to-spaxel variations. The final two columns show BLOSVD-2D reconstructions for two choices of σCAR. These display vertical stripes arising from spatially regularising over columns of spaxels. The degree of smoothing is greater for the dominant negative velocity component (top row) compared to the positive velocity component (bottom), and also for σCAR = 0.001 (fourth column) vs. σCAR = 0.03 (third). Spatial regularisation has proved most effective in the top right panel, where σCAR = 0.001 is lower than the density present in the channel map i.e. 0.002–0.004, as seen in the colour bar.

Figure 9 compares the uncertainties derived from BLOSVD-1D (blue) versus BLOSVD-2D (green). The example reconstructions at S/N = 20 (left column) and S/N = 200 (right) illustrate that BLOSVD-2D produces significantly smaller uncertainties that enclose the truth more tightly. At S/N = 20 (left column), the uncertainty intervals also shrink noticeably within the dashed vertical lines, but they remain largely unchanged outside these lines, that is, a reduction in uncertainty is more prominent within the distribution body than the wings. The central panel shows the uncertainties averaged over all spaxels for a range of S/N. In general, BLOSVD-2D significantly decreases the uncertainty in the LOSVD body at all S/N (solid lines). The additional dip in the BLOSVD-2D uncertainty for S/N = 200 arises from the difference in σ_CAR, which is 0.001 for S/N = 200 and 0.03 for all other S/N values. For the LOSVD wings (dashed lines), accounting for spatial correlations decreases uncertainty much less than in the body, but increasing the S/N has a much stronger effect.

Fig. 9. Uncertainty of LOSVDs recovered with BLOSVD (blue) and BLOSVD-2D (green). The side panels show an example from the central spaxel in Fig. 7: the true LOSVD (black line) is compared against the median recovery (solid line) and 99% credible intervals (shaded regions) for S/N 20 (left panels) and 200 (right). All four side panels share the same y-axis range. For 1D fits with BLOSVD (top left/right), the uncertainties are so large that almost the entire LOSVD body (i.e. within the vertical dashed lines) is consistent with zero; increasing from S/N 20 to 200 reduces uncertainty in the wings, but much less so in the body. BLOSVD-2D, in contrast, gives credible intervals that enclose the truth more tightly: for S/N 20 (bottom left), the dominant negative-velocity component is recovered with high statistical significance, while for S/N 200 (bottom right), this is also true for the weaker component. The central panel shows the uncertainty (measured by the width of the 99% credible intervals) averaged over all spaxels and velocity bins for a range of S/N. BLOSVD-2D shrinks the uncertainties in the LOSVD body at all S/N (solid lines). In the wings (dashed lines), the difference is weaker, but varying S/N affects the uncertainties much more strongly. A horizontal dotted line at y = 10−3 has been added in all five panels to aid comparison.

We end this section with a more detailed comparison of LOSVDs recovered with BLOSVD-2D and PNKR at high S/N. Fig. 6 already showed that PNKR outperforms BLOSVD-2D at S/N = 200. One reason for this seems to be the failure of BLOSVD-2D to effectively smooth the v_LOS > 0 component. A close comparison of S/N = 200 recoveries from PNKR-1D (red lines, left panel of Fig. 4) and BLOSVD-2D (green lines, left panel of Fig. 7) reveals another discrepancy. Focussing on spaxel 8, we see that BLOSVD-2D slightly over-predicts the amplitude of the v_LOS < 0 component and under-predicts the v_LOS > 0 component, while PNKR-1D predicts the relative amplitude of both components more accurately. Although this difference is small, a similar effect can be seen in most spaxels. To quantify this, we calculated the median bias over all 900 spaxels, finding that BLOSVD-2D over-predicts the v_LOS < 0 peak by 8% and under-predicts the v_LOS > 0 peak by 7%, while the equivalent values for PNKR-1D are 3% and 1%. One possible cause for this difference is that PNKR has the freedom to assign different populations to different kinematic components, while BLOSVD does not. We next ascertained whether PNKR is able to effectively capitalise on this additional freedom to recover dependences between stellar populations and kinematics.

4.3. Recovery of metallicity-velocity relations with PNKR

We examined the ability of PNKR to recover dependences between stellar populations and kinematics encoded in the 5D joint distribution f(x, v, z, t). This section shows results for PNKR only because stellar-population recovery is not possible with BLOSVD. Furthermore, all results shown here use PNKR-1D; as with the LOSVD recoveries in Sect. 4.2, we found no improvement when using PNKR-3D. This section focusses on reconstructions of metallicity-velocity relations; equivalent results on age-velocity relations are shown in Appendix E.

Figure 10 shows the PNKR recoveries of the metallicity-velocity relations. It shows LOSVDs from 9 spaxels coloured by the mean metallicity of stars at the given position and velocity (defined in Appendix C.3) for the ground-truth model (left grid) and the PNKR recovery at S/N = 200 (right grid). Overall, PNKR successfully recovers that the negative velocity component is more metal rich than the positive velocity component. Even for spaxels where the two components overlap (e.g. panel 1) and are therefore spectroscopically unresolved, we were able to recover a metallicity difference between the two components. While we correctly inferred the metallicity ordering of the two components, the absolute values of the recovered metallicity are less extreme than the true values. Additionally, some spaxels (e.g. 5, 8) show an erroneous upturn in metallicity at high positive v_LOS.

Fig. 10. Recovery of metallicity-velocity relations with PNKR. The 3 × 3 grids of sub-panels correspond to the spaxels highlighted in Fig. 3. LOSVDs are coloured by the mean metallicity of stars with a given vLOS, for the truth (left) and PNKR recovery at S/N 200 (right). Overall, PNKR successfully recovers that the negative-velocity component is more metal rich than the positive. This occurs where the two components are well separated in vLOS (e.g. spaxels 7-9), but also where they overlap (spaxel 1). Although the metallicity ordering of the two components is successfully inferred, the absolute values of the recovered metallicities are less extreme than the true values, and some spaxels (e.g. 5 and 8) show an erroneous upturn in recovered metallicity at high positive vLOS.

Figure 11 shows spatially integrated metallicity-velocity relations (defined in Appendix C.4) at different S/N levels. For velocities which contribute negligibly to the density (dotted-lines), metallicities are unconstrained and the recoveries diverge significantly from the truth. Focussing instead on the solid lines, Fig. 11 shows that at S/N = 200 (blue line) PNKR successfully recovers a sharp decline in metallicity, but this becomes washed out for lower S/N. Although the qualitative behaviour of the metallicity-velocity relation is well-recovered at S/N = 200, absolute metallicities are offset from the truths by 0.25/0.2 dex at negative/positive v_LOS. We defer the discussion of these absolute offsets to Sect. 5.2, and we next explore whether PNKR can infer the qualitative behaviour of metallicity-velocity relations in different scenarios.

Fig. 11. Spatially integrated metallicity-velocity relations. The true relation (black line) is compared against recoveries at varying S/N (coloured lines). The dotted lines show velocities with a negligible contribution to the galaxy (selected via a LOSVD threshold of 0.01) where the metallicity is unconstrained. A sharp drop in metallicity is recovered for S/N = 200, marginally so at S/N = 112, but it washes out for lower S/N.

Figure 12 shows PNKR recoveries of the spatially integrated metallicity-velocity relations for nine different ground-truth models. The models vary in the velocity separation (columns) and metallicity difference (rows) between the two components. For all nine models, we generated S/N = 200 mock data using popkinmocks (Jethwa 2023), and we ran PNKR-1D recoveries with settings as described in Appendix A.1. Fig. 12 shows that although absolute metallicities are often discrepant from the true values, in eight of nine cases, PNKR can successfully infer whether the metallicity-velocity relations are decreasing, flat, or rising. The sole exception is shown in the top right panel, where a flat relation is inferred while the truth decreases. This exception occurs where the two components overlap in velocity, presenting the most challenging recovery problem. The overall conclusion, however, is that PNKR generally can infer the qualitative behaviour of metallicity-velocity relations with S/N = 200 data.

Fig. 12. Recoveries of spatially integrated metallicity-velocity relations for nine different ground-truth models. The models vary in the separation of the two components in velocity space illustrated by the LOSVDs above each column (corresponding to spaxel 7 in Fig. 10), and the metallicity difference between the two components (in the top, middle and bottom rows, metallicity is a decreasing, flat or rising function of velocity). All nine recoveries were performed at S/N = 200. The dotted lines show velocities with a negligible contribution, where the metallicity is unconstrained.

5. Discussion

We have tested two algorithms (BLOSVD and PNKR) for non-parametric reconstructions of stellar population/kinematic distributions from IFS data. A key novelty of the two algorithms is that in addition to modes for 1D spectral fitting (BLOSVD-1D and PNKR-1D), they provide a novel functionality to fit multiple spectra simultaneously while accounting for spatial correlations (BLOSVD-2D and PNKR-3D). We begin this section with a summary of the results shown in Sect. 4 and then discuss our key findings related to BLOSVD (Sect. 5.1) and PNKR (Sect. 5.2). We then discuss the limitations of our experimental setup (Sect. 5.3) and finally summarise our practical recommendations (Sect. 5.4).

Kinematic recoveries with PNKR and BLOSVD were, to first order, successful, but both algorithms produced various artefacts and biases in different S/N regimes. The most accurate recoveries were achieved by PNKR at S/N > 40, but for a lower S/N, PNKR significantly over-smoothed the LOSVD bimodality. The 1D and 3D versions of PNKR produced indistinguishable results. For BLOSVD-1D (which uses no regularisation in velocity or spatial dimensions), we found spike artefacts for S/N >  100. The inclusion of spatial regularisation with BLOSVD-2D fixed the most severe spikes, but did not effectively smooth the subdominant v_LOS > 0 component. Additionally, PNKR and BLOSVD both produced flat wing artefact whose prominence decreased with S/N (∝S/N^−0.85). For S/N = 200, the BLOSVD recoveries show a bias in the relative amplitudes of the two kinematic components that is significantly reduced when using PNKR.

Stellar population recoveries with PNKR were biased towards central age/metallicity values in the SSP grid (i.e. less extreme than the truths) for all S/N. Nevertheless, for S/N = 200, we found that PNKR can recover the qualitative behaviour of metallicity-velocity relations, that is, whether the metallicity is a rising/falling/flat function of velocity. Stellar populations are not explicitly modelled in BLOSVD, but it adopts the implicit assumption that stellar populations and kinematics are locally independent, which likely explains the bias of the relative amplitudes of the two kinematic components seen at S/N = 200.

5.1. BLOSVD

Our tests with BLOSVD demonstrated the benefits of accounting for spatial correlations for reconstructing stellar kinematics. Our LOSVD reconstructions with BLOSVD-2D are more accurate (Fig. 6) and more precise (Fig. 9) than their 1D equivalents. The improvements in accuracy arise from a reduction in spike artefacts at high S/N and from subtler improvements to the LOSVD shape at lower S/N (Fig. 7). Regarding precision, we showed that BLOSVD-2D significantly reduced the LOSVD uncertainty at all S/N values we tested, and this occurred preferentially in the LOSVD body compared to the wings, as evidenced by the shrunken credible intervals in Fig. 9. These improvements might enable the detection of weak stellar components that would reach a relatively low statistical significance from 1D spectral fits, and we therefore advocate for the use of some spatial regularisation when modelling IFS data.

While our implementation of a spatial regularisation in BLOSVD-2D has demonstrated significant improvements, there is scope for further enhancements. The spatial regularisation has been much more effective in smoothing the dominant (v_LOS < 0) component than the sub-dominant (v_LOS > 0) component (Fig. 8). We observe that this is related to the σ_CAR hyper-parameter of the CAR prior in the sense that velocity channels with a density > σ_CAR are most effectively smoothed. This observation motivates a modification of our current implementation, which assigns all velocity bins the same σ_CAR, to a more flexible implementation with a velocity-dependent σ_CAR. We caution, however, that σ_CAR strongly affects the stability of the MCMC sampling (see Appendix A.4), so that added flexibility might require an additional tuning of the inference algorithm. Additionally, our current 2D experiment, which was restricted to fitting columns of spaxels, should be extended to fully 3D. This is a straightforward extension of our current implementation.

We showed BLOSVD reconstructions for kinematics, but not for stellar populations. As described in Sect. 2.2, this is because BLOSVD uses PCA to replace several hundred SSP templates with a handful of PCA templates. BLOSVD treats inferred PCA-weights as nuisance parameters, foregoing all information on stellar populations, but we note that this is not a fundamental limitation of this approach. Parzer et al. (2023) showed that PCA-weights can be de-projected to SSP-weights, which can be used to reconstruct physical quantities such as mean ages, metallicities, and star formation histories in a post-processing step. In future work, we will investigate stellar population properties derived in this way.

Despite its successes, we found that BLOSVD achieves less accurate LOSVD reconstructions than PNKR for S/N > 40 in its 1D or 2D varieties (Fig. 5). One reason for this is that BLOSVD-2D did not effectively smooth the v_LOS > 0 component. Additionally, BLOSVD imposes the ansatz that, locally, stellar-populations and kinematics are independent, that is, Eq. (3). As discussed in Sect. 4.2, this results in the over-prediciton by BLOSVD of the amplitude of the v_LOS > 0 component, while it under-predicts the v_LOS < 0 component. This trade-off arises because BLOSVD lacks the freedom to assign different stellar populations to the two components. PNKR imposes no such assumption, and our experiments with PNKR support the claim that this difference contributes to the relative underperformance of BLOSVD at high S/N.

5.2. PNKR

PNKR has the freedom to model the joint distribution over stellar populations and kinematics. The results of Sect. 4.3 show that with S/N = 200 data, PNKR effectively uses this freedom to infer the qualitative behaviour of metallicity-velocity relations. We demonstrated that PNKR can make these measurements locally (Fig. 10) and, by considering spatially integrated measurements, we showed that this requires S/N = 200 data (Fig. 11) and is robust against a variety of different ground-truth models (Fig. 12). These results explain why PNKR recovers the relative amplitudes of the two kinematic components more accurately than BLOSVD at S/N = 200, which contributes to the overall higher accuracy of the kinematic recoveries of PNKR at high S/N.

Our tests with PNKR demonstrated that it is possible to infer dependences between stellar populations and kinematics in a non-parametric way. To our knowledge, the only previous work to demonstrate non-parametric³ recoveries of joint population-kinematic distributions is STECKMAP (Ocvirk et al. 2006). The overall quality of our kinematic recoveries are broadly comparable with that achieved using STECKMAP. Regarding joint population-kinematic recoveries, Figs. 5 and 6 of Ocvirk et al. (2006) illustrate non-parametric recoveries of age-velocity distributions, successfully inferring distinct LOSVDs for components with ages 30 Myr and 3 Gyr. We extended this result by demonstrating that even given two relatively old components, both with ages > 1 Gyr, PNKR simultaneously recovered their kinematics and metallicities. Measurements such as these will be particularly useful for extragalactic archaeology (van de Ven et al. 2025), since by separating the properties of old stellar populations, we can learn about a galaxy’s evolution over timescales of billions of years.

While PNKR is able to recover the qualitative behaviour of metallicity-velocity relations at S/N = 200, our recovered metallicities are biased by ≈0.2 dex (Fig. 11). Validation studies of full-spectrum fitting (Koleva et al. 2008; Ge et al. 2018; Cid Fernandes 2018; Wang et al. 2024) typically recover mean metallicities with higher accuracy. For instance, Woo et al. (2024) reported a metallicity bias of only 0.07 dex when using pPXF (Cappellari & Emsellem 2004; Cappellari 2017). Compared to most validation studies, we modelled a relatively complex case of a galaxy with two counter-rotating components, but in additional tests with only a single component, the bias did not decrease significantly. Instead, this bias might be due to the iterative regularisation scheme of PNKR. Just as this caused oversmoothing of LOSVDs (right panel of Fig. 4), PNKR likely oversmooths the age-metallicity distributions, leading to mean metallicity and ages that are biased toward central values of the SSP grid. This phenomenon might also explain the erroneous upturns in metallicity seen at high positive velocities in panels 5 and 8 of Fig. 10: once again, the recoveries are biased towards the central values of the SSP grid. While we only observed LOSVD oversmoothing at low S/N, this might affect stellar-population distributions even at high S/N because there are fewer age/metallicity bins than velocity bins in the problem discretisation. We identified the PNKR wavelength deactivation scheme as a possible cause of oversmoothing (see the pink line in Fig. 5), hence modifications to this scheme (e.g. via wavelength importance sampling) may remedy the issue. Alternatively, using SSP grids with finer sampling and incorporating prior information on age-metallicity relations may also prove helpful.

We found that accounting for spatial correlations with PNKR-3D did not lead to any improvements compared to PNKR-1D. This conclusion was valid for LOSVD reconstructions (Fig. 6) and metallicity-velocity relations. The reason for this is unclear. H23 proposed two different approaches to incorporate spatial regularisation in PNKR. In the first approach, the unknown density f is built from hat-shaped basis functions that span across adjacent spaxels. This approach led to some qualitative improvements in LOSVD recoveries (see Fig. 5.6 of H23), but was prohibitively slow to run. For this work, we therefore adopted the alternative approach based on spatial smoothing filters described in Sect. 2.1 and named the reduced PNKR method in H23. Theoretical arguments suggest that the two approaches are equivalent (Hubmer et al. 2023), but we found that the second approach provides little benefit over 1D reconstructions.

5.3. Limitations of our experimental setup

We tested the specific scenario of a galaxy with a large-scale stellar counter-rotation. This scenario was chosen as it is complex enough to illustrate key features of our recovery algorithms, but it is far from typical. For example, in a systematic search of MaNGA galaxies, Bevacqua et al. (2022) found that just 64 of 4000 showed evidence of stellar counter-rotation. Nevertheless, with the improved spectral resolution of future IFS instruments (e.g. Thatte et al. 2016; Richard et al. 2024; Bryant et al. 2024), we will be able to identify far more typical varieties of stellar kinematic substructures, such as co-rotating disks, bars, and bulges. The benefits of the recovery algorithms shown in this work should generalise to these other situations, but dedicated tests will be useful to understand the limitations in different contexts.

The mock data we used in this work are idealised. We omitted instrumental effects including the PSF, LSF, and uncertain flux calibrations, as well as physical effects such as dust-attenuation, nebular emission, and sky residuals. While this suffices for the theoretical investigation we presented here, progress on real data will require us to deal with these additional effects. Several existing tools are available to create more realistic mock data (e.g. Harborne et al. 2020; Bottrell & Hani 2022; Cornejo-Cárdenas et al. 2025). BLOSVD and PNKR can both be extended to handle more complex forward models. Since the PNKR algorithm is tailored to the mathematical structure of the model, modifications to the forward model may require alterations of the inference algorithm. BLOSVD, on the other hand, is built in the framework of a probabilistic programming language, where the tasks of model specification and inference are decoupled. This means that additional effects can be added during the modelling process without the need to modify the inference algorithm.

5.4. Practical recommendations

The experiments presented in Sect. 4 allow us to provide a limited set of practical recommendations for full-spectrum fitting. Regarding kinematics, for moderate S/N (≲40) and/or applications focussed on on low-order kinematic moments, established tools such as pPXF (Cappellari & Emsellem 2004; Cappellari 2017) remain well suited. At higher S/N, where the data begin to support more detailed inferences, non-parametric approaches become warranted, but with insufficient regularisation, these are prone to produce spike artefacts that risk being misinterpreted as genuine kinematic structure. Existing tools for non-parametric LOSVD recovery (Ocvirk et al. 2006; Falcón-Barroso & Martig 2021; Mehrgan et al. 2023; Gasymov & Katkov 2024) provide options for regularisation in velocity space. Our results with BLOSVD showed that spatial regularisation is an effective alternative that may prove useful in cases where velocity regularisation risks oversmoothing kinematic substructure.

For stellar population recoveries, existing tools such as pPXF (Cappellari & Emsellem 2004; Cappellari 2017) and STARLIGHT (Cid Fernandes 2018) remain well suited for most current applications, but they are typically used under the assumption that stellar populations and kinematics are locally independent. Our tests with PNKR suggest that for S/N ≳ 200, it might be fruitful to drop this assumption and move towards a joint population-kinematic modelling. Since PNKR recovers absolute population parameters with significant bias, we currently recommend its use primarily as an exploratory tool for identifying qualitative population–kinematic trends rather than for precision population measurements. STECKMAP (Ocvirk et al. 2006) provides some similarity functionality, but we note that it is limited to recovering kinematics with an explicit dependence on age, but not on metallicity.

6. Conclusions

We have compared two algorithms for the reconstruction of stellar kinematics and populations from IFS data, PNKR and BLOSVD. In addition to modes for 1D spectral fitting, the two algorithms have extensions that simultaneously fit several spectra accounting for spatial correlations across spaxels. The PNKR algorithm Hinterer23 was the first tool developed for this task. We introduced a second approach that extends Bayes-LOSVD (Falcón-Barroso & Martig 2021) to incorporate spatial regularisation. We evaluated the algorithms using idealised mock data at varying S/N levels.

Accounting for spatial correlations led to significant improvements in the accuracy and precision of LOSVD recoveries with BLOSVD (Sect. 4.2). These benefits might enable the detection of weak stellar kinematic substructures that would evade detection from 1D spectral fits. We therefore advocate for further development and usage of 3D modelling tools for stellar recoveries. While we saw significant benefits in BLOSVD, 3D modelling with PNKR did not lead to any noticeable improvements over 1D fits. We conclude that the implementation of 3D modelling within PNKR, based on smoothing the solution between iterations of the optimisation algorithm, is less effective than that in BLOSVD, which couples adjacent spaxels together using a hierarchical Bayesian prior to encode spatial regularities. Although our implementation of spatial regularisation in PNKR proved ineffective, it still achieved smaller LOSVD reconstruction errors than BLOSVD for S/N > 40 (Fig. 6). This result can partly be explained by the other key novelty of the PNKR algorithm.

The PNKR algorithm has the freedom to model the full joint density in stellar populations and kinematics, forgoing the common ansatz that these are locally independent. With a sufficiently high S/N, we found that PNKR makes use of this freedom to successfully recover the qualitative behaviour of metallicity-velocity relations (Sect. 4.3). Although the qualitative behaviour was well recovered, absolute metallicities were recovered with a significant bias, which we speculate was caused by the iterative regularisation scheme of PNKR. Nevertheless, this is a promising result that suggests that we might be able to separate the population properties of stellar substructures in a wider range of contexts than has been attempted previously.

This work has demonstrated two promising strategies to better detect and characterise stellar structures in galaxies: joint population-kinematic modelling with PNKR, and 3D modelling with BLOSVD. Future development efforts will focus on incorporating these ideas into a unified tool that will be accompanied by an open-source software release. In particular, extending the BLOSVD formalism to allow for joint population-kinematic modelling offers a promising path toward a flexible user-friendly tool.

Data availability

Existing code for PNKR is provided at https://github.com/prashjet/PNKR_python for reproducibility, while the 3D extension of BLOSVD will be made available in the future software release.

Acknowledgments

This research was funded in part by the Austrian Science Fund (FWF) SFB 10.55776/F68 “Tomography Across the Scales”, project F6805-N36 (Tomography in Astronomy) and F6811-N36 (Advancing Extragalactic Archaeology through Novel Inversion Techniques). For open access purposes, the authors have applied a CC BY public copyright license to any author-accepted manuscript version arising from this submission. The computational results have been achieved using the Austrian Scientific Computing (ASC) infrastructure. The authors thank Dr. Fabian Hinterer for his contributions.

References

Appendix A: Additional details

This appendix contains implementation details and comments on the convergence of the recoveries presented in Sect. 4.

A.1. PNKR-1D

For PNKR-1D recoveries in Sect. 4.1, we use step-size 100 and tolerance of 1.3 to determine wavelength deactivation i.e. if a wavelength is fit to within 30% accuracy it is deactivated. All PNKR-1D fits reach zero active equations before the algorithm terminates. The number of iterations required to reach this point increased with S/N, from 74 iterations for S/N = 20 to 560 for S/N = 200, with computational run times correspondingly increasing from 1 hour at S/N = 20 to 15 hours for S/N = 200. This behaviour is expected from the theory of iterative regularisation methods (Kaltenbacher et al. 2008). Intuitively, this is because for small S/N we cannot iterate many times before the we begin to over-fit to noise and worsen reconstructions, while larger S/N requires more iterations to reach this point. In attempts to to fix LOSVD oversmoothing for S/N = 20, we additionally ran recoveries with smaller step-sizes (50, 10) and a less generous tolerance (1.1) but found these led to no improvement.

A.2. PNKR-3D

For PNKR-3D recoveries, we use the same setting and tolerance listed in Appendix A.1. For the parameters of the smoothing kernel, we use a (5×5) grid over shape parameter α and scale-length β. We take 5 values of α spaced linearly between 0.2 and 1.0. For S/N = 200, we take 5 values of β logarithmically spaced between 10⁻⁵ and 10⁻³. For lower S/N, we inversely scale the β values with S/N, so that lower S/N data can benefit from more spatial regularisation.

Figure A.1 shows the change in LOSVD recovery error for PNKR-3D versus PNKR-1D. For all S/N and for all choices of α, the difference in recovery error shows a U-shape: for small β, the error difference approaches 0, while for large β, PNKR-3D becomes oversmoothed and it’s error exceeds PNKR-1D. These U-shape curves reach a minimum for some intermediate β. This gives us confidence that the range of smoothing parameters we have explored is sufficient i.e. our conclusion that PNKR-3D provides little benefit is not limited by this range.

Fig. A.1. Change in LOSVD recovery error when using PNKR-3D vs PNKR-1D (y-axis) vs the kernel scale length β (x-axes) and power-law index α (colours as shown in legend). The five panels show this for different S/N levels, decreasing from left to right.

Regarding convergence, in contrast to PNKR-1D fits, where we always reached zero active equations, PNKR-3D fits instead stopped at a plateau defined by no decrease in the number of active equations over the last 100 iterations. The number of iterations required to reach the plateau increases with S/N, from 200 at S/N = 20 to 350 for S/N = 200, with small variations (approx 10) depending on the smoothing kernel parameters. For S/N = 200 we ran additional tests with a plateau defined by no decrease in the number of active equations over the last 300 iterations, but found no change.

A.3. BLOSVD-1D

For BLOSVD-1D, we use 15 PCA templates to represent the SSP grid (which accounts for > 99% of the variance of the full grid), and run the NUTS MCMC sampler for 1000 warm-up steps followed by 2000 sampling steps. For all S/N, BLOSVD-1D fits successfully achieve the converge criteria in Sect. 2.2). Run times were 30 minutes at high S/N and 3 minutes at low S/N. These times are much smaller than those for PNKR-1D since we can parallelise over the spaxels.

A.4. BLOSVD-2D

For BLOSVD-2D we use again use 15 PCA templates to represent the SSP grid, and again run the NUTS MCMC sampler for 1000 warm-up steps followed by 2000 sampling steps. In this case, however, when running the NUTS sampler, we found that using including all spaxels led to a prohibitively slow warm-up stage. We therefore limited ourselves to using only a fifth of the spaxels (i.e. 6 out of 30 spaxels per column) for the warm-up stage, took the resulting post warm-up settings and used these to initiate the sampling stage including all spaxels. The initial positions of the MCMC chains are posterior samples from 1D fits from each spaxel.

The CAR hyper-parameters discussed in Sect. 2.2 are set as follows. For the correlation parameter 0 < ρ_CAR < 1 we take a hyper-prior strongly peaked at ρ_CAR = 1 to promote smoothness. Concretely, we take the beta distribution p(ρ_CAR) = Beta(ρ_CAR; 4, 1) which has a single mode at ρ_CAR = 1 and a median at ρ_CAR = 0.84. We encountered poor convergence diagnostics when we used a hyper-prior on the second CAR hyper-parameter σ_CAR. Specifically, when using a uniform prior on log σ_CAR, we encountered many divergent transitions when MCMC chains approached σ_CAR = 0. When truncating the prior below some value σ_CAR > 0, the resulting samples invariably bunched up tightly against the truncation limit. To have more control over these experiments, we therefore abandoned the use of a hyper-prior on σ_CAR and instead tried two fixed values σ_CAR ∈ [0.001, 0.03].

The convergence of the resulting BLOSVD-2D depended on the choice of σ_CAR. All fits converged successfully for σ_CAR = 0.03. For σ_CAR = 0.001 many fits failed, exhibiting poor mixing and small effective sample sizes. Out 30 columns of spaxels, we encountered [30,6,5,1,0] failures for our five S/N bins increasing from 20 to 200 i.e. all columns failed at S/N = 20 while all succeeded for S/N = 200. For the results shown in Sect. 4 we required all columns to converge successfully, therefore for S/N < 200 we show only results with σ_CAR = 0.03.

Appendix B: Quality of fit

Figures B.1 to B.6 visualise the quality of fit for (1) PNKR-1D at S/N = 200, (2) PNKR-1D at S/N = 20, (3) BLOSVD-1D at S/N = 200, (4) BLOSVD-1D at S/N = 20, (5) BLOSVD-2D at S/N = 200, (6) BLOSVD-2D at S/N = 20. In each case, the right panels show two spectra selected at the best fit (blue, bottom row) and worst fit (red, top row) spaxels, with the residuals normalised by σ in the middle row. The images on the left show χ² maps between observed and reconstructed datacubes, displayed as the difference k from the number of wavelength bins N in units of $\sqrt{2N}$. In these units, a perfect fit would centre at 〈k〉 = 0 with a variance 〈k₂〉 = 1 and display no spatial structure. Both BLOSVD-1D and 2D approach this ideal scenario for S/N = 20 (Figs. B.4 and B.6); however, at S/N = 200 the χ² maps show some spatial structure. PNKR-1D reconstructions (Figs. B.1 and B.2) lie further from the observations than BLOSVD-1D with significantly more spatial structure. This is to be expected due to PNKR’s implicit regularisation scheme, which trades fidelity of the data reconstruction for accuracy in the reconstruction of the unknown density. Nevertheless, in all cases, the example spectra in the right panels show that the reconstructed spectra (black lines) do tend to trace the data (coloured dots), giving us confidence that all algorithms tested are functioning as expected.

Fig. B.1. Quality of fit for PNKR-1D at S/N = 200. Left: On-sky map of the χ2 difference between observed and reconstructed datacubes (in units k as described in Appendix B). Top right and bottom right: Observed spectra (coloured dots) compared with reconstructions (black lines). These are shown in blue (red) for two spaxels selected to have the best (worst) fits, whose locations are shown by coloured circles in the left panel. Middle-right: Fit residuals normalised by the uncertainty σ, summarised as a histogram in the small, rightmost panel.

Fig. B.2. Same as Fig. B.1 but for PNKR-1D at S/N = 20.

Fig. B.3. Same as Fig. B.1 but for BLOSVD-1D at S/N = 200.

Fig. B.4. Same as Fig. B.1 but for BLOSVD-1D at S/N = 20.

Fig. B.5. Same as Fig. B.1 but for BLOSVD-2D at S/N = 200.

Fig. B.6. Same as Fig. B.1 but for BLOSVD-2D at S/N = 20.

Appendix C: Definitions

This appendix defines several quantities referenced throughout this work.

C.1. Light-weighted LOSVDs

The distribution f(x, v, z, t) defined in Eq. (1) is mass-weighted, since SSPs S(λ, z, t) are defined as spectra per unit-mass of a stellar population. The equivalent light-weighted distribution f^ℒ is achieved by scaling by the integrated flux of the SSPs, i.e.

$$\begin{array}{c}\hfill f^{\mathcal{L}}(\mathit{x},v,z,t)\propto f(\mathit{x},v,z,t){\displaystyle \int S}(\lambda ,z,t)\mathrm{d}\lambda .\end{array}$$(C.1)

The light-weighted LOSVD at position x is then given by

$$\begin{array}{c}\hfill f^{\mathcal{L}}(v|\mathit{x})=\frac{\int f^{\mathcal{L}}(\mathit{x},v,z,t)\mathrm{d}z\mathrm{d}t}{\int f^{\mathcal{L}}(\mathit{x},v,z,t)\mathrm{d}v\mathrm{d}z\mathrm{d}t}.\end{array}$$(C.2)

Note that BLOSVD adopts a simplifying assumption - Eq. (3) – which imposes that, locally, stellar populations and kinematics are independent, and therefore light-weighted or mass-weighted LOSVDs are, by definition, identical.

C.2. LOSVD recovery error

Given the true light-weighted LOSVDs f_true^ℒ(v|x) and a recovery f_rec^ℒ(v|x), we define the error between them as

$$\begin{array}{c}\hfill \mathrm{Error}=\frac{1}{N_{\mathrm{spax}}}{\displaystyle \int \int \left|\frac{f_{\mathrm{true}}^{\mathcal{L}}(v|\mathit{x})-f_{\mathrm{rec}}^{\mathcal{L}}(v|\mathit{x})}{f_{\mathrm{true}}^{\mathcal{L}}(v|\mathit{x})}\right|}{f}_{\mathrm{true}}^{\mathcal{L}}(v|\mathit{x})\mathrm{d}\mathit{x}\mathrm{d}v.\end{array}$$(C.3)

This can be interpreted as the fractional error, weighted by the true LOSVD, averaged over all spaxels and velocity bins. Note that - by definition in Eq. (C.2) - the true LOSVDs are probability density functions over v at each x, and hence they are appropriately normalised to serve as weighting factors in this definition. In practice, we calculate this error ignoring the two cancelling factors of f_true^ℒ(v|x), to avoid numerical issues, and approximate the integrals as discrete sums over spaxels and velocity bins.

C.3. Local velocity-metallicity relations

The local velocity-metallicity relations shown in Sect. 4.3 are light-weighted mean metallicities conditional on position and velocity - i.e. 𝔼(z|v, x) - evaluated for some fixed choice of position x. To define 𝔼(z|v, x) we start from the light-weighted joint density f^ℒ(x, v, z, t) defined in Eq. (C.1). The light-weighted metallicity distribution conditional on velocity and position is then given by

$$\begin{array}{c}\hfill f^{\mathcal{L}}(z|v,\mathit{x})=\frac{\int f^{\mathcal{L}}(\mathit{x},v,z,t)\mathrm{d}t}{\int f^{\mathcal{L}}(\mathit{x},v,z,t)\mathrm{d}z\mathrm{d}t}.\end{array}$$(C.4)

Finally, the light-weighted mean metallicity conditional on velocity and position is given by

$$\begin{array}{c}\hfill \mathbb{E}(z|v,\mathit{x})={\displaystyle \int z}{f}^{\mathcal{L}}(z|v,\mathit{x})\mathrm{d}z.\end{array}$$(C.5)

C.4. Global velocity-metallicity relation

The global velocity-metallicity relations shown in Sect. 4.3 are light-weighted mean metallicities conditional on velocity. This is calculated in a similar way as the local relations defined in the previous section, with additional integrals over position in the numerator and denominator of Eq. (C.4).

Appendix D: Wing artefacts

Figure D.1 shows the prominence of the flat wing artefacts as a function of S/N for PNKR (1D & 3D) and BLOSVD (1D & 2D). As described in Sect. 4.1, we measure the wing prominence as the median value of recovered LOSVDs for |v_LOS|> 600 km/s (i.e. outside the maximum velocity of the ground truth model) averaged over all spaxels. For all recovery algorithms, wing prominence decreases with S/N in a manner well described by a power law with wing prominence ∝S/N^−0.85. Though the power-law slope is common between all algorithms, the normalisation of the wing prominence at fixed S/N varies between the algorithms by ≈10%. The similarity of all four algorithms in Fig. D.1 suggests that a power-law suppression of LOSVD wings with S/N may be a general outcome of non-parametric kinematic recoveries, however the exact slope will likely vary in different contexts.

Fig. D.1. Prominence of the flat wing artefacts as a function of S/N for our four recovery algorithms (coloured lines). For all algorithms, wing prominence decreases with S/N in a manner well described by a power-law with slope -0.85 (dashed black line).

Appendix E: Age-velocity relations

This appendix shows age-velocity relations equivalent to the metallicity-velocity relations shown in Sect. 4.3. Figure E.1 shows the local relations (equivalent to Fig. 10) while Fig. E.2 shows global relations at a range of S/N (equivalent to Fig. 11). For both cases, the results are qualitatively similar to those seen in in Sect. 4.3 however the tendency for PNKR to push the recovered ages towards the middle of the age range present in the SSP grid is stronger than for metallicity. For example, in Fig. E.1, PNKR correctly infers that the dominant component is youngest, but the distinction between the ages of the two components is much blurrier than the clear metallicity difference seen in Fig. 10.

Fig. E.1. Local age-velocity relations recovered with PNKR-1D at S/N = 200.

Fig. E.2. Global age-velocity relations recovered with PNKR-1D over a range of S/N.

All Figures

	Fig. 1. Toy demonstration of LOSVD recovery with PNKR. Left: Number of active wavelengths per iteration of the algorithm, with iterations 300 and 1000 highlighted. Right: Corresponding LOSVD recoveries. At iteration 300, the recovery is smooth and matches the truth well, while extending to iteration 1000 introduces noise.
In the text

Fig. 2. Demonstration of LOSVD recovery with BLOSVD-1D. Each panel shows BLOSVD-1D recoveries via their posterior median (thick blue line), 99% credible interval (shaded region), and ten randomly selected posterior samples (thin blue lines). For the left panel, the template spectrum contains a single absorption line exactly one pixel wide; here, the recovery tightly encloses the truth. For the right panel, the absorption line is two pixels wide, which introduces a degeneracy between neighbouring velocity bins, evidenced by the zig-zagging posterior samples. This degeneracy results in noisy median recoveries and inflated credible intervals.

In the text

Fig. 3. Our ground-truth galaxy model in age-metallicity space (left) as a light-weighted image (centre) and mean velocity map (right). The model consist of two components. The dominant component (blue contours) is a young metal-rich thin disk with negative LOS velocity. The counter-rotating weaker component (orange contours) is older, more metal poor, and more extended vertically. Successive contours in the central panel show changes in flux by 15%. The grid of crosses in the right panel indicates the spaxels used for illustration in Figs. 4, 7, and 10.

In the text

Fig. 4. Comparison of LOSVDs recovered from 1D spectral fits. True LOSVDs (black) are shown alongside recoveries from PNKR-1D (red) and BLOSVD-1D (blue). The 3 × 3 grids in the sub-panels correspond to the spaxels highlighted in Fig. 3. At S/N 200 (left grid), PNKR recovers the LOSVD well everywhere, while BLOSVD recoveries show spike artefacts in some cases e.g. panels 1, 2 and 9. At S/N 20 (right grid) the recoveries are worse. Both algorithms produce LOSVDs with extended, flat wings which reach the edge of the velocity range. For spaxels 5-9 we see that BLOSVD-1D captures the dip between the two components whereas PNKR-1D oversmooths this feature.

In the text

	Fig. 5. Error in recovered LOSVDs from 1D spectral fits as a function of S/N. Comparing the two methods, PNKR (red) achieves the smallest error at high S/N, while BLOSVD (blue) is best at low S/N. By modifying PNKR to retain all wavelengths throughout the recovery (pink line), we see an improvement at low S/N, with a recovery error similar to BLOSVD, at the expense of larger errors at high S/N.
In the text

Fig. 6. Error in recovered LOSVDs from 2D/3D fits as a function of S/N. The solid lines are repeated from Fig. 5, showing the results of 1D spectral fits with PNKR (red) and BLOSVD (blue). The dashed lines show results when fitting multiple spaxels simultaneously and promoting spatial smoothness in the recoveries. PNKR-3D (dashed red), which fits the full 3D datacube, produces results that are almost indistinguishable from the 1D PNKR fits. In contrast, BLOSVD-2D (dashed green), which fits columns of spaxels simultaneously, achieves a significant improvement over the 1D case for all S/N values.

In the text

Fig. 7. Comparison of LOSVDs recovered from with BLOSVD-1D (blue lines) against BLOSVD-2D (green lines). At S/N 200 (left grid), the spike artefacts seen in the 1D fits are alleviated in the 2D fits (e.g. panels 1, 2, 5, and 9). At S/N 20 (right grid), the 2D fits improve the recovered LOSVD shapes in many cases, particularly for the dominant component with vLOS < 0 (e.g. panels 4, 6, and 9). At high and low S/N, the effect of 2D fitting is weaker for the subdominant component with vLOS > 0.

In the text

Fig. 8. Channel maps showing the spatial distribution of stars in velocity bins labelled on the left. The first column shows true distributions, which are smooth by construction. The remaining columns show reconstructions for S/N = 200 data. The second column shows BLOSVD-1D reconstructions, which show large spaxel-to-spaxel variations. The final two columns show BLOSVD-2D reconstructions for two choices of σCAR. These display vertical stripes arising from spatially regularising over columns of spaxels. The degree of smoothing is greater for the dominant negative velocity component (top row) compared to the positive velocity component (bottom), and also for σCAR = 0.001 (fourth column) vs. σCAR = 0.03 (third). Spatial regularisation has proved most effective in the top right panel, where σCAR = 0.001 is lower than the density present in the channel map i.e. 0.002–0.004, as seen in the colour bar.

In the text

Fig. 9. Uncertainty of LOSVDs recovered with BLOSVD (blue) and BLOSVD-2D (green). The side panels show an example from the central spaxel in Fig. 7: the true LOSVD (black line) is compared against the median recovery (solid line) and 99% credible intervals (shaded regions) for S/N 20 (left panels) and 200 (right). All four side panels share the same y-axis range. For 1D fits with BLOSVD (top left/right), the uncertainties are so large that almost the entire LOSVD body (i.e. within the vertical dashed lines) is consistent with zero; increasing from S/N 20 to 200 reduces uncertainty in the wings, but much less so in the body. BLOSVD-2D, in contrast, gives credible intervals that enclose the truth more tightly: for S/N 20 (bottom left), the dominant negative-velocity component is recovered with high statistical significance, while for S/N 200 (bottom right), this is also true for the weaker component. The central panel shows the uncertainty (measured by the width of the 99% credible intervals) averaged over all spaxels and velocity bins for a range of S/N. BLOSVD-2D shrinks the uncertainties in the LOSVD body at all S/N (solid lines). In the wings (dashed lines), the difference is weaker, but varying S/N affects the uncertainties much more strongly. A horizontal dotted line at y = 10−3 has been added in all five panels to aid comparison.

In the text

Fig. 10. Recovery of metallicity-velocity relations with PNKR. The 3 × 3 grids of sub-panels correspond to the spaxels highlighted in Fig. 3. LOSVDs are coloured by the mean metallicity of stars with a given vLOS, for the truth (left) and PNKR recovery at S/N 200 (right). Overall, PNKR successfully recovers that the negative-velocity component is more metal rich than the positive. This occurs where the two components are well separated in vLOS (e.g. spaxels 7-9), but also where they overlap (spaxel 1). Although the metallicity ordering of the two components is successfully inferred, the absolute values of the recovered metallicities are less extreme than the true values, and some spaxels (e.g. 5 and 8) show an erroneous upturn in recovered metallicity at high positive vLOS.

In the text

Fig. 11. Spatially integrated metallicity-velocity relations. The true relation (black line) is compared against recoveries at varying S/N (coloured lines). The dotted lines show velocities with a negligible contribution to the galaxy (selected via a LOSVD threshold of 0.01) where the metallicity is unconstrained. A sharp drop in metallicity is recovered for S/N = 200, marginally so at S/N = 112, but it washes out for lower S/N.

In the text

Fig. 12. Recoveries of spatially integrated metallicity-velocity relations for nine different ground-truth models. The models vary in the separation of the two components in velocity space illustrated by the LOSVDs above each column (corresponding to spaxel 7 in Fig. 10), and the metallicity difference between the two components (in the top, middle and bottom rows, metallicity is a decreasing, flat or rising function of velocity). All nine recoveries were performed at S/N = 200. The dotted lines show velocities with a negligible contribution, where the metallicity is unconstrained.

In the text

	Fig. A.1. Change in LOSVD recovery error when using PNKR-3D vs PNKR-1D (y-axis) vs the kernel scale length β (x-axes) and power-law index α (colours as shown in legend). The five panels show this for different S/N levels, decreasing from left to right.
In the text

Fig. B.1. Quality of fit for PNKR-1D at S/N = 200. Left: On-sky map of the χ2 difference between observed and reconstructed datacubes (in units k as described in Appendix B). Top right and bottom right: Observed spectra (coloured dots) compared with reconstructions (black lines). These are shown in blue (red) for two spaxels selected to have the best (worst) fits, whose locations are shown by coloured circles in the left panel. Middle-right: Fit residuals normalised by the uncertainty σ, summarised as a histogram in the small, rightmost panel.

In the text

	Fig. B.2. Same as Fig. B.1 but for PNKR-1D at S/N = 20.
In the text

	Fig. B.3. Same as Fig. B.1 but for BLOSVD-1D at S/N = 200.
In the text

	Fig. B.4. Same as Fig. B.1 but for BLOSVD-1D at S/N = 20.
In the text

	Fig. B.5. Same as Fig. B.1 but for BLOSVD-2D at S/N = 200.
In the text

	Fig. B.6. Same as Fig. B.1 but for BLOSVD-2D at S/N = 20.
In the text

	Fig. D.1. Prominence of the flat wing artefacts as a function of S/N for our four recovery algorithms (coloured lines). For all algorithms, wing prominence decreases with S/N in a manner well described by a power-law with slope -0.85 (dashed black line).
In the text

	Fig. E.1. Local age-velocity relations recovered with PNKR-1D at S/N = 200.
In the text

	Fig. E.2. Global age-velocity relations recovered with PNKR-1D over a range of S/N.
In the text