Covariance-regularized reconstruction of data cubes in integral field spectroscopy and application to MaNGA data

TL;DR

This paper introduces a new covariance-regularized resampling method for integral field spectroscopy data that achieves near-intrinsic resolution and low pixel correlations, improving upon Shepard's method.

Contribution

A novel resampling technique that enhances resolution and reduces pixel correlations in irregularly sampled integral field data.

Findings

Achieves ~16% better resolution than Shepard's method.

Produces images with resolution close to the instrument's intrinsic limit.

Low statistical correlations between neighboring pixels.

Abstract

Integral field spectroscopy can map astronomical objects spatially and spectroscopically. Due to instrumental and atmospheric effects, it is common for integral field instruments to yield a sampling of the sky image that is both irregular and wavelength-dependent. Most subsequent analysis procedures require a regular, wavelength independent sampling (for example a fixed rectangular grid), and thus an initial step of fundamental importance is to resample the data onto a new grid. The best possible resampling would produce a well-sampled image, with a resolution equal to that imposed by the intrinsic spatial resolution of the instrument, telescope, and atmosphere, and with no statistical correlations between neighboring pixels. A standard method in the field to produce a regular set of samples from an irregular set of samples is Shepard's method, but Shepard's method typically yields images with a degraded resolution and large statistical correlations between pixels. Here we introduce a new method, which improves on Shepard's method in both these respects. We apply this method to data from the Mapping Nearby Galaxies at Apache Point Observatory survey, part of Sloan Digital Sky Survey IV, demonstrating a full-width half maximum close to that of the intrinsic resolution (and ~ 16% better than Shepard's method) and low statistical correlations between pixels. These results nearly achieve the ideal resampling. This method can have broader applications to other integral field data sets and to other astronomical data sets (such as dithered images) with irregular sampling.

Tables

Table 1. Table 1: Sharpness measures of PSF.

	pixel scale per pixel	Kernel	Shepard’s	CRR
FWHM(arcsec)	0.5	1.905	2.316	1.973
Strehl ratio	0.5	1	0.664	0.996
FWHM(arcsec)	0.75	1.905	2.317	1.973
Strehl ratio	0.75	1	0.664	0.996

Table 2. Table 2: Correlation coefficient between pixels with different separations (quadratic mean of values between the central pixel and its surrounding pixels), for Shepard’s method and CRR.

	pixel scale arcsec/pixel	1 pixel	2 pixels
Shepard’s	0.5	0.876	0.703
Shepard’s	0.75	0.767	0.465
CRR	0.5	0.303	0.094
CRR	0.75	0.050	0.024

Table 3. Table 3: Radial profile measurements for real galaxies. FWHM is in units of arcsec, while the central flux is in units of 10 − 17 superscript 10 17 10^{-17} erg/s/cm / 2 Å {}^{2}/\mbox{\AA} /arcsec 2 .

	method	H$\alpha$ 6563Å		O[III] 5007Å
		FWHM	center	FWHM	center
8718-12701	CRR	2.511	79.5	2.181	42.9
	Shepard	2.729	62.7	2.644	29.0
8549-12701	CRR	2.193	405	2.184	201.6
	Shepard	2.640	278	2.591	202
8482-12704	CRR	2.283	181	2.264	37.8
	Shepard	2.729	128	2.729	26.2

Equations

$$S[j]=\sum _{i=1}^{N}W[j,i]f[i]$$

$$W[j,i]=\frac{b[i]}{W_0[j]}\exp \left(-\frac{r[i,j{]}^2}{2{\sigma }_0^2}\right),$$

$$W_0[j]=\sum _{i}b[i]\exp \left(-\frac{r[i,j{]}^2}{2{\sigma }_0^2}\right),$$

$$m[i]=\sum A[i,j]F[j],$$

$$A[i,j]=K_i(x_i-X_j,y_i-Y_j)$$

$${\chi }^2=\sum _i\frac{{\left(f[i]-m[i]\right)}^2}{\sigma [i{]}^2}.$$

$$\vec{F}=({\mathbf{A}}^{\mathrm{T}}{\mathbf{N}}^{-1}\mathbf{A}{)}^{-1}{\mathbf{A}}^{\mathrm{T}}{\mathbf{N}}^{-1}\vec{f},$$

$$\mathbf{C}=⟨\Delta \vec{F}\Delta {\vec{F}}^{\mathrm{T}}⟩=({\mathbf{A}}^{\mathrm{T}}{\mathbf{N}}^{-1}\mathbf{A}{)}^{-1}.$$

$${\mathbf{C}}^{-1}=\mathbf{PD}{\mathbf{P}}^{-1}=(\mathbf{P}{\mathbf{D}}^{\frac{1}{2}}{\mathbf{P}}^{-1})(\mathbf{P}{\mathbf{D}}^{\frac{1}{2}}{\mathbf{P}}^{-1})=\mathbf{QQ},$$

$${\mathbf{N}}^{-\frac{1}{2}}\mathbf{A}=\mathbf{U\Sigma }{\mathbf{V}}^{\mathrm{T}}.$$

$$\vec{F}=\mathbf{V}{\mathbf{\Sigma }}^{-1}{\mathbf{U}}^T{\mathbf{N}}^{-\frac{1}{2}}\vec{f}$$

$${\mathbf{C}}^{-1}=\mathbf{V\Sigma \Sigma }{\mathbf{V}}^{\mathrm{T}}$$

$$\mathbf{Q}=\mathbf{V\Sigma }{\mathbf{V}}^{\mathrm{T}}.$$

$$R[i,j]=\frac{1}{{\sum }_{j}Q[i,j]}Q[i,j],$$

$${\mathbf{C}}_G=\mathbf{RC}{\mathbf{R}}^{\mathrm{T}}$$

$$C_G[i,i]={\left(\sum _{j}Q[i,j]\right)}^{-2}$$

$$\vec{G}$$

$$=$$

$$\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}$$

$${\chi }^2=\sum _i\frac{[f[i]-m[i]){]}^2}{\sigma [i{]}^2}+{\lambda }^2\sum _{j}F[j{]}^2$$

$$\vec{F}$$

$$=$$

$${\Sigma }_{\ast ,ii}^{-1}=\frac{{\Sigma }_{ii}}{{\Sigma }_{ii}^2+{\lambda }^2}$$

$$\vec{G}=\mathbf{RV}{\mathbf{\Sigma }}_{\ast }^{-1}{\mathbf{U}}^{\mathrm{T}}{\mathbf{N}}^{-\frac{1}{2}}\vec{f}$$

$$\vec{G}=\mathbf{W}\vec{f},$$

$$\mathbf{W}=\mathbf{RV}{\mathbf{\Sigma }}_{\ast }^{-1}{\mathbf{U}}^{\mathrm{T}}{\tilde{\mathbf{N}}}^{-\frac{1}{2}}$$

$${\mathbf{C}}_G=⟨\mathbf{W}\Delta \vec{f}\Delta {\vec{f}}^{\mathrm{T}}{\mathbf{W}}^{\mathrm{T}}⟩=\mathbf{WN}{\mathbf{W}}^{\mathrm{T}}.$$

$$\tilde{\mathbf{C}}=\mathbf{W}\tilde{\mathbf{N}}{\mathbf{W}}^T.$$

$$\frac{{\sigma }_2}{{\sigma }_1}$$

$$\frac{A_2}{A_1}$$

$$\text{FWHM}=1.05{\text{ FWHM}}_1=2.473{\sigma }_1$$

$$S/N=\frac{f_{\max }\cdot s}{\sqrt{f_{\max }\cdot s}}=\sqrt{f_{\max }\cdot s}$$

Full text

††pubyear: 2015††pagerange: Covariance-regularized reconstruction of data cubes in integral field spectroscopy and application to MaNGA data–Covariance-regularized reconstruction of data cubes in integral field spectroscopy and application to MaNGA data

1 Introduction

Integral field spectroscopy (IFS) yields a rich set of information about extended objects such as galaxies. Modern facilities for performing IFS include the Multi Unit Spectroscopic Explorer (MUSE; Laurent et al, 2006), the Keck Cosmic Web Imager (KCWI; Morrissey et al, 2012), and several efforts to create large samples of nearby galaxies, such as the Calar Alto Legacy Integral Field Area Survey (CALIFA; Sánchez et al, 2012), the Sydney-AAO Multi-object IFS (SAMI; Croom et al, 2012) and Mapping Nearby Galaxies at Apache Point Observatory (MaNGA; Bundy et al, 2014).

Here, we consider the MaNGA survey, one of three core programs in the fourth-generation Sloan Digital Sky Survey (SDSS-IV; Blanton et al, 2017). MaNGA uses integral field units (IFUs) consisting of optical fiber bundles to obtain spectra across the face of a sample of 10,000 low redshift galaxies, making it possible to map the spectroscopic properties of galaxies and to interpret this spectroscopy in terms of two-dimensional maps of stellar age, gas phase and stellar phase elemental abundances, star formation histories, and kinematics.

The MaNGA fiber bundles are arranged in a hexagonal grid with a separation of 151 microns. Each fiber has a cladding and buffer in its outer annulus and has an active core size of 120 microns, corresponding to 2 arcsec in the Sloan Foundation Telescope focal plane. This configuration results in an effective filling factor for the active cores of 56% (Law et al, 2015). With a typical seeing full-width half-maximum (FWHM) of 1.5 arcsec at Apache Point Observatory (APO), a single observation undersamples the sky image considerably. Thus, to increase the sampling of the point spread function (PSF) and to avoid sampling irregularities, MaNGA uses dithered observations. Ideally, each observation consists of three dithered exposures, performing a set of subsequent exposures on each side of an equilateral triangle. This pattern of observations leads to a finer, but still hexagonal, pattern on the sky at the guiding wavelength. In the ideal case when all observations are taken at the same hour angle, the dithers produce a hexagonal pattern at all wavelengths, with an overall wavelength-dependent shift of the hexagonal pattern on the sky due to chromatic differential refraction by the atmosphere. If the observations are taken at different hour angles, the set of dithers will not produce a perfect hexagonal pattern of samples at all wavelengths. Thus, generally speaking, the MaNGA data produce an irregular sampling of fluxes on the sky that varies with wavelength.

Most analysis techniques are designed to handle regular, usually rectangular, grids. Except for analyses involving a full “forward modeling” of the MaNGA data, resampling the fluxes onto such a regular grid is important. The resulting product is referred to as a data cube. These data cubes are the primary product of the MaNGA data reduction pipeline (DRP; Law et al 2016). The DRP extracts the spectra and calibrates them, producing row-stacked spectra (RSS) — one spectrum per fiber all interpolated onto the same wavelength grid. At each wavelength, MaNGA then uses a modified version of Shepard’s method (Shepard 1968) to resample the RSS data for that wavelength onto a rectangular grid. This method, which is a flux-conserving variation of Shepard’s interpolation method, is widely adopted because of its simplicity and robustness (Yang et al, 2004), for example in the CALIFA pipeline (e.g. Sánchez et al, 2012).

In this paper we reconsider this choice of resampling method and propose an alternative. Before beginning, we ask: what is the best performance we should expect from a resampling algorithm? Here we define this “best performance” as being equivalent to actually sampling the real sky with the real instrument at the locations of every grid point in the desired resampling. This process would produce an image with a resolution equivalent to the intrinsic spatial resolution of the instrument, telescope, and atmosphere, and zero statistical correlations between the errors in each pixel. This definition means we are not attempting to deconvolve the intrinsic resolution — just to resample with as little loss of resolution or introduction of statistical correlations as possible.

Relative to this ideal, Shepard’s method suffers two problems. First, it broadens the final reconstructed PSF substantially. In the cases we consider in this paper, the FWHM of the final PSF can be 20% broader than the intrinsic resolution. Second, the errors in pixels in the reconstruction are highly correlated with each other. In the cases we consider in this paper, for Shepard’s method pixels separated by 1 arcsec can have correlation coefficients between their errors as high as 0.70, complicating any correct analysis of the images that accounts for errors. These deficiencies drive us to find a way to improve the performance of the resampling.

To do so, we make use of the techniques developed in a different context by Bolton et al (2010), which they term “spectroperfectionism.” This name alludes to the fact that they are performing a spectroscopic extraction. However, we prefer the name “covariance-regularized reconstruction” as more generally applicable, since the approach of Bolton et al (2010) is appropriate for a number of other contexts beyond spectroscopic extraction or (as applied here) integral field spectroscopy image reconstruction.

We will make some important distinctions between different terms in this paper to avoid confusion:

• •

Kernel: The response $K(x,y)$ to a delta-function source of the atmosphere, telesecope, and instrument combined. In the specific case of MaNGA, the kernel is the atmospheric seeing, convolved with the telescope’s optical response, convolved with the fiber profile (which we approximate as a 2 arcsec diameter top hat).

• •

Kernel-convolved image: The actual or model image on the sky convolved with the kernel.

• •

Shepard’s image: The resampling resulting from Shepard’s method. Referred to as a vector of pixel values $\vec{S}$.

• •

Deconvolved reconstruction: A model of the image on the sky with the kernel deconvolved. Referred to as a vector of pixel values $\vec{F}$ (never used directly in the final result).

• •

Covariance-Regularized Reconstruction (CRR): The resampling resulting from our method. Referred to as a vector of pixel values $\vec{G}$.

• •

Point spread function (PSF): In our usage, the PSF will refer to the response to a delta-function on the sky of the output image from the analysis, either Shepard’s image or the CRR. In an ideal reconstruction, the PSF would be identical to the kernel values at each pixel (using the delta-function location as the center of the kernel function).

This article is organized as follows. Section 2 describes the resampling methods we consider here. Section 3 presents a series of tests of this method for a simulation. Section 4 presents a demonstration of the method for real data. Section 5 summarizes and discusses the results.

2 Image Resampling Methods

2.1 Shepard’s method, the MaNGA DRP standard

A specialized version of Shepard’s method is described in Sánchez et al (2012) in the context of CALIFA and Law et al (2016) in the context of MaNGA. The input data from the IFU survey consists of the flux intensity $f[i]$ and the variance of the error $N[i]=\sigma[i]^{2}=\langle\Delta f[i]^{2}\rangle$ in each fiber. For a system with $N_{\rm fiber}$ fibers, with which we have taken $N_{\rm exp}$ exposures, for each wavelength channel there are $N=N_{\text{fiber}}\times N_{\text{exp}}$ values, each corresponding to a different location on the sky. We will refer below to each such observation as a fiber-exposure. For the MaNGA DRP, the output grid is rectangular with a pixel size 0.5 arcsec per pixel, and we write the total number of pixels as $M$.

The transformation from intensities at the irregularly located fiber locations to the Shepard’s image is:

[TABLE]

where $\vec{S}$ is Shepard’s image, and the $M\times N$ matrix $W[j,i]$ is the weight of each fiber location $i$ contributing to the output grid point $j$. In Shepard’s method, the weight function is a circularly symmetric Gaussian that depends on the distance $r[i,j]$ between the fiber location $i$ and the grid point $j$:

[TABLE]

for $r[i,j]<r_{\rm lim}$, and zero otherwise. $\sigma_{0}$ defines the width of the Gaussian function, and for the MaNGA DRP is set to 0.7 arcsec. The MaNGA DRP takes $r_{\mathrm{lim}}=1.6$ arcsec. The normalization parameter $W_{0}[j]$ is defined as the sum of the $N$ weights for each output grid point $j$, to guarantee the conservation of flux:

[TABLE]

over all pixels $i$ for which $r[i,j]<r_{\rm lim}$. In Equations 2 and 3, $b[i]$ is zero if the inverse variance $N[i]^{-1}=0$ and is unity otherwise.

2.2 Covariance-regularized reconstruction

Our method can be written in the same form as Shepard’s method (Equation 1), but with a different choice of $W[j,i]$. The choice is motivated in the following manner, which explicitly designs the weights so that the final reconstruction remains consistent with the samples and so that its covariance matrix has small off-diagonal entries. We refer to an image using our method as a covariance-regularized reconstruction (CRR).

Consider fitting a linear model to reconstruct the observables $f[i]$. Our model consists of a set of delta functions with fluxes $F[j]$, distributed on a regular grid (the same grid we want to use for the reconstruction). Each fiber-exposure $i$ observes this function convolved with the kernel appropriate for that particular observation, $K_{i}(x,y)$. In the case of MaNGA, the kernel is the convolution of the atmospheric seeing, the telescope optical response, and the fiber. This model for the observables can be written as

[TABLE]

where:

[TABLE]

is the value of the kernel at the separation between the model pixel $j$ at $(X_{j},Y_{j})$ and the fiber-exposure $i$ at $(x_{i},y_{i})$. The kernel is a known property of the observations based on the estimated seeing at the time of each observation.

We can fit for the model parameters $F[j]$ by minimizing the $\chi^{2}$ error:

[TABLE]

Because the model is linear in $F[j]$, the solution can be written:

[TABLE]

where $\mathbf{N}$ is the data covariance matrix, which in our case is diagonal with diagonal entries $\sigma[i]^{2}$.

Therefore, we can estimate $\vec{F}$, which is analogous to a deconvolved image based on the fiber samples, because it is the sky image before convolution with the kernel. As one expects for a deconvolution, the error covariance for $\vec{F}$ is highly non-diagonal — there are strong correlations and anticorrelations in the errors of neighboring pixels. The covariance matrix is

[TABLE]

Using $\vec{F}$ for science is very undesirable because of these correlations, which would complicate any error analysis but also lead to large fluctuations among the values of $\vec{F}$.

A common technique is to regularize the values of $\vec{F}$, either under a Tikhonov regularization, a maximum entropy criterion, or something else (e.g. Warren & Dye 2003). However, here we take a different approach, which is to regularize the covariance such that it is diagonal, which turns out to be similar to reconvolving $\vec{F}$ to a resolution similar to the kernel.

The covariance matrix can be whitened through a linear transformation of $\vec{F}$ that can be found by taking a square root of the inverse covariance matrix $\mathbf{C}^{-1}=\mathbf{A}^{\mathrm{T}}\mathbf{N}^{-1}\mathbf{A}$. For this symmetric and positive definite matrix, we can take its square root by finding its eigensystem:

[TABLE]

where $\mathbf{D}$ here is a diagonal matrix of eigenvalues and $\mathbf{P}$ is the matrix of eigenvectors. For $\mathbf{D}^{\frac{1}{2}}$, the positive root is always chosen.

Another path to finding the same matrix $\mathbf{Q}$ is through the singular value decomposition (SVD):

[TABLE]

where $\mathbf{\Sigma}$ is diagonal, $\mathbf{V}$ is orthogonal ($\mathbf{V}\mathbf{V}^{\mathrm{T}}=\mathbf{1}$), and $\mathbf{U}$ is close to orthogonal ($\mathbf{U}\mathbf{U}^{\mathrm{T}}$ is diagonal with ones for dimensions $i$ with $\Sigma_{i}\neq 0$ or zeros for dimensions $i$ with $\Sigma_{i}=0$). In this case:

[TABLE]

if $\Sigma_{i}\neq 0$ for all $i$. We can also use the standard Moore-Penrose inverse technique and set $\Sigma_{i}^{-1}=0$ for $\Sigma_{i}=0$, which allows us to handle truly degenerate cases smoothly. We can show that

[TABLE]

and therefore,

[TABLE]

Thus, we can find $\mathbf{Q}$ without ever constructing the full covariance matrix or its inverse, or explicitly finding its eigenvectors, which may be useful in cases when the covariance matrix is ill-conditioned.

The linear transformation of $\vec{F}$ that we will use is not quite $\mathbf{Q}$, because we want the transformation to conserve the total flux in $\vec{F}$. We can achieve this goal by normalizing over each row to give the transformation matrix:

[TABLE]

$\mathbf{R}$ is then a linear transformation under which $\mathbf{C}$ can be diagonalized:

[TABLE]

with the entries of the diagonal matrix $\mathbf{C}_{G}$ given by

[TABLE]

(In detail, the resulting covariance matrix is not precisely diagonal if any $\Sigma_{i}=0$). This result means that if we define:

[TABLE]

Then $\mathbf{C}_{G}$ is the covariance matrix of $\mathbf{G}$. Multiplying by the matrix $\mathbf{R}$ turns out to be similar to a reconvolution of $\vec{F}$ with the kernel, although it is not strictly speaking a convolution.

$\vec{G}$ will be our reconstruction, which by design should have close to a diagonal covariance matrix. We will show later that in the context of MaNGA it also has a sharper PSF than Shepard’s method. Like Shepard’s method it is just a linear combination of the input fluxes so can be written in the same form as Equation 1.

There are some adjustments we will make to this method. The first adjustment is that we will apply a regularization term to handle singular values. This adjustment will make the procedure more numerically robust but have almost no impact on the results. We can think of this regularization term in terms of an adjustment to the function we are minimizing in Equation 6 with a quadratic term, the simplest version of Tikhonov regularization:

[TABLE]

The solution $\vec{F}$ is altered under this regularization to:

[TABLE]

where $\mathbf{\Sigma_{\ast}}$ is diagonal and:

[TABLE]

and the solution $\vec{G}$ becomes:

[TABLE]

The second adjustment is that we will not use the actual flux noise vector $\vec{\sigma}$ estimated from the data. If we did so, then there would be a correlation between the flux in a fiber-exposure (which determines the true noise vector) and the values in the corresponding column of $\mathbf{W}$. This correlation will lead to a dependence of the PSF in the CRR on the signal that will greatly complicate the interpretation of the image. In addition, it is usually the case that the noise is estimated from the signal and using it as a weight in the fit will therefore be biased. For these reasons, in the calculation, we use $\mathbf{\tilde{N}}$, which equals unity where $\mathbf{N}^{-1}$ is not zero, and zero where $\mathbf{N}^{-1}$ is zero.

A third adjustment we make is to slightly alter the kernel to remove contributions for which the radius is larger than 4 arcsec. This choice makes little difference in the final result but makes the calculation considerably faster.

The fourth adjustment we make is to remove from the SVD calculation any model pixels $j$ that are further than some distance $r_{\rm lim}$ (which we here set to 1.6 arcsec, the same as DRP’s choice) from any fiber — basically, any pixels not well-constrained by the data — which lowers the condition number of $\mathbf{A}$ and makes the SVD more stable.

The final method then can be written as follows:

[TABLE]

where:

[TABLE]

where $\mathbf{V}$, $\mathbf{U}$, and $\Sigma$ are from the SVD of $\mathbf{A}$, $\Sigma_{\ast}$ is defined in Equation 24, and $\mathbf{R}$ is defined by Equations 13 and 14. The covariance matrix of $\vec{G}$ is then:

[TABLE]

Here, $\mathbf{N}$ is the diagonal variance matrix from the individual fibers. This covariance matrix is not guaranteed to be diagonal, but because of how we have constructed the weights it will prove to be much closer to diagonal than Shepard’s method.

When operating on the resulting data cube for science, it is important to identify unreliable spaxels about which there is little information provided by the fiber data. The MaNGA pipelines use two masks, NOCOV to indicate that there is no information about the pixel, and LOWCOV to indicate that there is a low amount of information about the pixel. These masks are based on the fiber-level maskbits, which record the effect of hot pixels, cosmic ray hits, and other effects that can make the data from a fiber unusable at some or all wavelengths. Our approach to defining these masks for the data cube spaxels is different from the standard MaNGA pipeline. The MaNGA pipeline counts the total contribution of fibers to each pixel, based on the ${\mathbf{W}}$ matrix for Shepard’s method, to identify poorly constrained spaxels. We use a different approach, which is to use the diagonal elements of the covariance matrix when we assume constant, unit noise in each fiber that is not masked:

[TABLE]

This covariance matrix is not the same as the actual covariance ${\mathbf{C}}_{G}$. We set spaxels with a variance more than twice the median variance to LOWCOV, so that they may be ignored in scientific analysis. This procedure appropriately masks spaxels too near the edge of the fiber bundle or that are affected by bad fibers.

There are several free parameters in the CRR method as applied to MaNGA, which we now discuss. First, there is the shape of the kernel, which is determined by our estimate of the observing conditions at each exposure. Second, there are the conditions for dropping edge pixels. Third, there is the regularization parameter $\lambda$. Fourth, there is the pixel scale of CRR. The first three parameters prove to have no significant affect on the results for MaNGA, as we show later in this paper. We will explore below the effect of the reconstruction pixel scale.

Some general aspects of the method, compared to Shepard’s method, are worth noting before describing those tests. Although our method is also just a linear combination of the fluxes, the weights are not determined by a stationary function as they are for Shepard’s method. They are also not restricted to be nonnegative. These properties, particularly the latter, are essential to reducing the off-diagonal covariances and producing a PSF close to the kernel resolution.

Although we motivated the method based on the model fit to the parameters $\vec{F}$ expressed in Equation 6 and the diagonalization of their covariance matrix $\mathbf{C}$, we never need to explicitly determine either $\vec{F}$ or $\mathbf{C}$. We will nevertheless calculate these quantities below in order to demonstrate their properties.

Our approach is mathematically identical to the approach proposed by Bolton et al (2010) in the different context of spectroscopic extraction. Our notation differs somewhat from theirs. Specifically, their $\vec{p}$ is our $\vec{f}$, their $\vec{f}$ is our $\vec{F}$, and their $\tilde{f}$ is our $\vec{G}$.

3 Tests on Simulated Data

In this section, we use simulations to characterize the performance of Shepard’s method, which is what the MaNGA DRP uses, and our CRR method.

3.1 General information

MaNGA uses a hexagonal pattern of fibers to detect flux as a function of position on the sky. The hexagonal pattern is dithered in position between different exposures, and varies on the sky with wavelength. Our goal is to resample the fluxes onto a rectilinear grid. In Figure 1, the blue points show the fiber locations at $\lambda=5500$Å for the plate-IFU 8720-1901, shown as a function of $X$ and $Y$ position in the focal plane (corresponding to RA and Dec). The red points show the pixels we are using in the reconstruction. As explained above, we exclude pixels when their distance to all the fibers is larger than 1.6 arcsec. In this figure, we show both 0.5 arcsec/pixel and 0.75 arcsec/pixel scales. We will examine the effect of difference choices of pixel scales on our reconstruction later.

We will simulate observations for a point source to test the methods. Since the methods we test are linear, the point source responses at different locations can fully characterize the performance. The flux for each fiber is sampled from the kernel function $K$ as the response of a point source at some location (in most cases, we use $X=Y=0$).

Panel (a) of Figure 2 shows our adopted kernel for wavelength $\lambda=5500\mbox{\AA}$. We use a double Gaussian function to simulate the atmospheric seeing, and convolve with the fiber profile (a 2 arcsec top-hat) to generate the kernel function. The normalized double Gaussian is defined by the standard deviation $\sigma_{1}$ for the inner Gaussian, the ratio $\sigma_{2}/\sigma_{1}$ of the outer to inner Gaussian standard deviation, and the ratio $A_{2}/A_{1}$ of the central values of the outer and inner Gaussians. We use a typical pair of ratios that approximates the inner parts of a Moffat-like atmospheric PSF (Law et al 2015; Jim Gunn, private communication):

[TABLE]

The FWHM of the resulting model PSF is:

[TABLE]

The MaNGA data reduction reports the FWHM of atmospheric seeing for each exposure at the guider wavelength $\lambda_{0}=5400\mbox{\AA}$. For other wavelengths, we will assume the seeing varies as $\lambda^{-1/5}$ (Yan et al, 2016).

For each simulation, the fiber locations are provided by a MaNGA observation with several exposures, for which positions accounting for a variety of observational effects (bundle metrology, dithering, chromatic and field differential refraction, etc.) are stored in the MaNGA data products, for example those shown in Figure 1. Normally the number of MaNGA exposures ranges from 6–21 in order to achieve uniformity and the required $S/N$ ratio.

Each pixel in the output image is a weighted sum of the values sampled for each fiber-exposure. Figure 3 shows an example for the central pixel. For Shepard’s weights (right panel), the weights for each fiber-exposure are just a decreasing function of distance from the pixel. For the CRR, the weights are both positive and negative and have an oscillatory nature, similar to that found in sinc-interpolation methods.

Whereas the fluxes sampled from the kernel are noiseless, the fluxes of the actual observations are not. The main sources of noise are Poisson noise in the number of electrons due to the object, sky, and dark current, plus the read noise from the amplifiers. For the bulk of locations and wavelengths in MaNGA, the noise is object-dominated. When we include simulated noise in our tests, we concentrate on this regime, so that the noise is proportional to the square root of the flux. We characterize the signal-to-noise ratio of the simulations based on that of the fiber with the maximum simulated (noiseless) flux $f_{\rm max}$. For a chosen signal-to-noise ratio we then define a scale factor $s$ converting flux to number of photons $N=fs$, so that:

[TABLE]

Then for each fiber we apply Poisson noise based on the resulting $N$ for each fiber.

3.2 Nominal case

We start our tests with the nominal case at 0.5 arcsec/pixel, in order to compare directly at the MaNGA DRP pixel scale. Figure 2 shows the response of a point source in the center, averaged over the $g$-band wavelength range. Panel (b), (c), and (d) show, respectively, MaNGA’s DRP datacube (which uses Shepard’s method), the CRR image, and our implementation of Shepard’s image, all at 0.5 arcsec/pixel. Panel (e) and (f) show the CRR image and Shepard’s image at 0.75 arcsec/pixel. As we said previously, panel (b) and (d) are almost the same because they use the same method. We show the DRP result in (b) to confirm that we are analyzing the data consistently.

We compare the performance of the methods quantitatively by measuring the FWHM and also a pseudo-Strehl ratio for the resulting PSF. We define the pseudo-Strehl ratio as the ratio of the peak image data cube slice intensity compared to the center of the kernel-convolved point source image. This is the ratio of the peak intensity to the largest possible given the resolution of the instrument. We show here a single typical wavelength slice at $\lambda=5500\mbox{\AA}$ with seeing FWHM around 1.19 arcsec, but we find similar results throughout the whole spectrum and various plate-IFUs.

Figure 4 shows the profile of the intensity of the resulting images. To determine a FWHM reliably, we fit a model to the pixel values. In particular, we use the kernel function $K$ as our model, which is parametrized by a choice of seeing. We use least-squares to find the best-fitting seeing value, and infer the FWHM and the pseudo-Strehl ratio from the corresponding kernel. In analyzing the MaNGA DRP results, Law et al (2016) used a Gaussian fit instead; we find that this is adequate for a broad PSF but not for the CRR image PSF. The results for plate-IFU 8720-1901, including the FWHM and the Strehl ratio for the kernel, for Shepard’s method, and for the CRR method, are listed in Table 1.

For consistency, given a pair of Strehl ratios or FWHM values $v_{1}$ and $v_{2}$ for the Shepard and CRR images, we define the improvement in this quantity as $\frac{2|v_{2}-v_{1}|}{v_{2}+v_{1}}$. The radial profile shows that there is a 16.0% percent improvement in the FWHM between Shepard’s method and CRR, which turns out to be a typical level of improvement. Denothing the kernel FWHM as $v_{0}$, the maximum possible improvement of the FWHM would be $\frac{2(v_{2}-v_{0})}{v_{2}+v_{0}}=19.5\%$ at most. The pseudo-Strehl ratio increase is around 40%. Our results are very close to the best possible performance without performing a deconvolution (i.e. very close to the kernel). We note in passing that we were surprised that the kernel had a FWHM slightly less than 2 arcsec, given that the fiber top-hat FWHM is precisely 2 arcsec; however, this result is correct — slightly blurring the top-hat with a Gaussian slightly reduces the FWHM.

Next we compare the covariance between pixels for both methods. The covariance between pixels of the final image is expressed in Equation 28, for both our method or Shepard’s method, depending only on which weights $\mathbf{W}$ are used. In order to factor out the magnitude of the covariance along the diagonal, we will examine the correlation matrix between pixels, defined in the usual way as $\rho_{jk}=C_{jk}/\sqrt{C_{jj}C_{kk}}$.

Figure 5 shows the resulting correlation matrix between pixels. We omit masked pixels, which are mostly at the corners. The image scaling is arcsinh to best examine the off-diagonal components. Shepard’s output image has a broad covariance matrix, while our method is nearly diagonal, as expected from the previous section. We show pixel scales of 0.5 and 0.75 arcsec/pixel, and for the CRR covariances the difference in the level of correlation between these cases is clear. We consider the choice of pixel size and its effect on covariance more quantitatively in the next subsection.

3.3 Choice of pixel size

Here we discuss the choice of pixel size. The primary considerations are the sharpness of the resulting PSF (as quantified by FWHM and pseudo-Strehl ratio), the covariance between pixels, and the sampling of the image. The results shown in this section for plate-IFU 8720-1901 and $\lambda=5500$ Å are representative of what we find at other wavelengths and in other plate-IFUs.

The sharpness of the PSF differs very little between the two pixel scales, as Table 1 demonstrates.

The correlation coefficients do tend to depend on pixel scale, as shown in Table 2, based on the central pixel of plate-IFU 8720-1901 at 5500 Å. In this table, we list the quadratic means of the correlation coefficients between the central pixel and the pixels which are separated from it by 1 pixel or 2 pixels. We examine the correlations more fully in Figure 6, which shows the full correlation between the central pixel and all others for both methods with the two pixel scales. Each image corresponds to a row of the correlation matrix. The CRR image outperforms the Shepard’s image in either case. The CRR image has some ringing in its correlation matrix that is greatly reduced at a pixel scale of 0.75 arcsec/pixel. Statistical independence is desirable so that the fluxes in the cube can be used to fit models and to propagate errors without tracking a broad covariance matrix, which is complex and burdensome.

The final major consideration for pixel size is whether the pixels provide sufficient sampling of the image. At 0.75 arcsec/pixel, the FWHM of the images are about 2.5 pixels, which is above the usual rough guide for critical sampling. Although for unusual PSFs of the sort produced by this instrument there may yet be some poorly sampled power at 0.75 arcsec/pixel, sampling more densely leads to an increase in covariance between pixels, as shown above.

These considerations of sampling and covariance lead us to 0.75 arcsec/pixel as the best choice for our application. This pixel scale provides decent sampling without inducing excess correlations between pixels or ringing behavior.

3.4 PSF as a function of source position

The PSF response for both CRR and Shepard’s is sensitive to the relative position on the sky of the point source and the fibers. We therefore need to characterize the PSF response across the face of IFU and verify that we satisfy MaNGA’s requirement that the PSF FWHM vary by less than 10% across the IFU. In Figure 7 we consider the FWHM homogeneity and in Figure 8 we consider the pseudo-Strehl ratio homogeneity. Below, we also discuss the homogeneity of the axis ratio $b/a$ of the PSF.

Each figure uses a grid of locations of the point source, using the 0.75 arcsec/pixel CRR and the 0.5 arcsec/pixel Shepard’s method results (to match MaNGA’s implementation). We quantify the variation of each quantity relative to its median:

[TABLE]

Here, $\sigma$ is the FWHM or the pseudo-Strehl ratio, and $\sigma_{0}$ is its median value. As noted above, the CRR image is consistently higher resolution than the Shepard’s image for point sources regardless of where they are relative to the fibers. In addition, in this case the CRR image also shows less fractional variation in resolution than Shepard’s method, especially in the FWHM. In the example plate-IFU 8720-1901, the fractional standard deviation of the FWHM is 0.033 for Shepard’s method, compared to 0.023 for our method. Meanwhile, the fractional standard deviation of the pseudo-Strehl ratio is 0.056 in Shepard’s method, compared to 0.048 in our method. Thus, the CRR image slightly outperforms Shepard’s image in terms of PSF homogeneity for this fiber bundle. This comparison depends on the location of the fiber-exposures; for the other fiber bundles we have tested, the CRR images always give similar or better homogeneity than Shepard’s images for the FWHM and pseudo-Strehl ratio.

However, for the PSF roundness we find that Shepard’s method outperforms CRR. We quantify the roundness with the axis ratio $b/a$ of the PSF. For each point source location within a radial distance to the center equalling to the FWHM, we fit the PSF with a 2D Gaussian. The parameters of this fit yield the minor-to-major axis ratio. Following (Law et al, 2015), we quantify the performance with $(b/a)_{99}$, which is the axis ratio for which 99% of the point source locations yield a PSF with $b/a>(b/a)_{99}$. This quantity is a conservative lower limit on the PSF roundness and its variation. In order to compare more directly to Law et al (2015), we implemented our test with constant seeing. In our analysis, Shepard’s image has $(b/a)_{99}=0.91$, whereas the CRR image has $(b/a)_{99}=0.85$.

We were unable to reproduce the results of (Law et al, 2015) for Shepard’s method, who found a substantially higher value of $(b/a)_{99}=0.96$ using a very similar methodology to ours. This is partly because the coordinates we adopt is from an actual observational fiber bundle. If we use a three-point dithered hexagonal grid with seeing FWHM=1.4 arcsec to do the simulation, we got $(b/a)_{99}=0.942$ with Shepard’s reconstruction method. The homogeneity statistics parameters are consistently better than with the real observational coordinates. But we also found Shepard’s method outperforms our result in the homogeneity.

3.5 Adding noise

The simulated data we test on above is noiseless. However, we need to test our performance in the presence of realistic noise. The sources of noise in MaNGA are outlined in Section 3.1, but as noted there we will only consider the effect of Poisson noise due to the object signal, which is the most problematic (increasing with the signal rather than staying constant) and which is usually dominant.

Figure 9 shows a point source with noise, quantified by the signal-to-noise ratio of peak flux in the reconstruction image, for both the CRR and Shepard’s image. Obviously the image exhibits some noise. For Shepard’s image, the result looks smoother, due to the larger off-diagonal covariances — i.e. the pixel-to-pixel fluctuations are reduced because neighboring pixels are correlated.

Figure 10 examines the FWHM and pseudo-Strehl ratios under a range of noise conditions. The CRR image remains sharper and brighter in the center regardless of the noise level.

3.6 Effect of an inaccurate PSF model

We build our analysis based on the assumption that the kernel is the same as our assumed one — that is that our double Gaussian model for the seeing with a given width, convolved with the fiber profile, is correct. However, the seeing width is an estimate and its FWHM can deviate from reality as much as 20%; other aspects of the kernel model may be incorrect as well.

Shepard’s method does not use any information about the kernel, and thus its behavior is independent of how accurately we know the kernel. However, our CRR method depends on a kernel model to determine the weights. Therefore we need to test whether the CRR image varies under reasonable assumptions about the inaccuracy of our kernel model.

To test the behavior of the algorithm under these conditions, we varied the FWHM of the assumed seeing from the actual value in our simulations, and compared the CRR image with the case that the assumed seeing was the same as the actual seeing. The results are shown in Figure 11. We varied the seeing by as much as 30% from the observational estimate for this plate of 1.19 arcsec. The FWHM and pseudo-Strehl ratio for this case are 1.973 arcsec and 0.996, when the seeing is correctly estimated. For a 10% error in our assumed seeing, the FWHM of CRR image changes by around 0.30%; for a 20% error in our assumed seeing, the FWHM of CRR image changes by around 0.40%. The fractional variation of the Strehl ratio is also similarly small. Therefore our results are relatively insensitive to whether or not we know the kernel exactly. We have checked that this variation is small independent of the actual seeing in the simulation (e.g. if the actual seeing is poor rather than the relatively good value of 1.19 arcsec).

3.7 Regularization parameter

One of the free parameters in the method is the regularization parameter $\lambda$. Here we examine the effect of this parameter on the results. Figure 12 shows results without regularization in the upper panels and with regularization in the lower panels, for comparison. The left panels show the deconvolved reconstruction $\vec{F}$ in the upper left panel (using 0.75 arcsec/pixel) for a simulated point source, showing significant ringing, with large correlations and anti-correlations between pixels. The upper center panel quantifies the covariance between the central pixel of $\vec{F}$ and all the other pixels (the central row of the covariance matrix), showing large anticorrelations between the central pixel and its immediate neighbors.

These correlations and anticorrelations are large and can introduce numerical instability to our analysis, since $\vec{F}$ results from the last four factors in Equation 27. However, with appropriate regularization, the behavior of $F$ (and therefore those last three factors) is much more constrained. While improving the numerical stability, this regularization will not affect the extraction of our final reconstruction from the fiber data, due to the first factor $\mathbf{R}$ in Equation 27.

The right panels show the contribution of a particular fiber to all the pixels on the grid, equivalent to a column of the $\mathbf{W}$ matrix. Unlike the Shepard’s method weights, the contributions are not all positive and exhibit some ringing. This behavior is characteristic of accurate interpolation kernels, such as sinc-interpolation, so it is expected. These weights remain identical even under regularization, meaning that for the values of $\lambda$ we use (or smaller), it does not affect the final results.

3.8 Consistency across wavelength

The locations of fibers and the kernel shape are a function of wavelength, which causes a necessary variation with wavelength of the weights $\mathbf{W}$ in our method. Here we examine whether the CRR method produces spectra that are consistent across wavelength under this variation.

We simulate a constant $f_{\lambda}$ point source for all wavelengths and use the meta-data (positions, atmosphere conditions) for plate-IFU 8720-1901. We compare the output spectra for CRR and Shepard’s method in Figure 13 for the total flux within the bundle. In Figure 14, we perform the same comparison for just the spectrum in the central pixel.

Figure 13 shows that the variation of the sum of all the pixels in the simulation results is within about 3% of constant. Once the overall difference in amplitude of $\sim 1\%$ between the CRR and Shepard’s method are accounted for, the relative variation is consistent between the methods to about 0.5%. The variation in both methods is due to not having enough dithers of exposures, which inevitably leads to a variation in the flux due to the changing position of the point source relative to the fibers as a function of wavelength. Artificially adding in more dithers can reduce this variation for the simulations (though of course we cannot do that for the observations!).

Figure 14 shows the variation of the central pixel intensity across wavelength. The variation of intensity is 6% for the CRR result, compared to 3% for Shepard’s result. Therefore, for an individual pixel there is a slightly larger spectrophotometric inconsistency in the CRR method than in Shepard’s method.

For the real case, besides the non-constant flux and extended fiber configuration, we need to consider the bad fibers identified by the MaNGA DRP for each plate-IFU, leading to the “low coverage” or “no coverage” masks. These fiber masks occur because of cosmic ray events, bad flat fields, CCD defects, broken fibers, or 2D extraction problems (Law et al, 2016). They will result in zero weights in the $\mathbf{W}$ matrix. But we have verified that they do not affect the spectrophotometric consistency as a function of wavelength, because of the normalization over the contributions from all the good fibers.

We applied our method on 140 fiber bundles in order to demonstrate that the example plate-IFU (8720-1901) is typical. In each case, we simulate a point source at the center of the cube. We compare the PSF FWHM and central intensity between CRR and Shepard’s image at several different wavelengths. As shown in Figure 15 and Figure 16, the reconstruction in our method is generally narrower than Shepard’s result and with higher pseudo-Strehl ratio.

4 Demonstration with actual data

The tests above regard a simulated point source. In this section, we apply the method to several plate-IFUs in MaNGA as a demonstration and to verify that it is working as expected. We will defer a thorough battery of tests of real data to the next phase of this work.

Figure 17 shows four test galaxies. As the $irg$-band color images in the left column show, two are elliptical galaxies and two are spiral galaxies. The first two rows (plate-IFUs 8720-1901 and 8143-6101) are AGN, and the second two (plate-IFUs 8247-6101 and 9183-12702) are non-AGN. We use these different type of galaxies just to show the variety of possible cases.

The remaining columns, from left to right, show the H$\alpha$ emission in the CRR method and in Shepard’s method, the [O III] 5008 emission in both methods, and the mean continuum of the 5300–6000 Å region in both methods. The H$\alpha$ and [O III] line emission fluxes are estimated by subtracting a continuum estimate using side bands around each line, and an unweighted integration of the line flux in a fixed rest frame wavelength range around each line. This relatively crude method is adequate to characterize the image quality, and we defer to the next phase of this work a more careful analysis with the MaNGA Data Analysis Pipeline (Westfall et al 2019).

In each image, the sharper nature of CRR image is clear. This clarity is most dramatically shown in the H$\alpha$ emission in the bottom row and in the [O III] emission in the second row. It is also clear that the off-diagonal covariance in the CRR images is lower, manifesting as a noisier-looking image. For example, in the [O III] images in the bottom two rows, the regions without significant emission show uniform white noise in the CRR image, but show the characteristic mottling of correlated noise in Shepard’s method.

We also check the profile as a function of radius from the center for real images. Considering that the morphology of the galaxy will affect the radial profile, it is best to choose point-like sources for the comparison. In Figure 18, we select three targets identified as AGN by Rembold et al (2017). Using the same methods we use for the PSF fitting, we measure the FWHM and central flux. Table 3 shows the results. Only the comparison of FWHM and central intensity for different methods has meaning, not the absolute value of these quantities, since the galaxies are more extended than a PSF and may have more irregular morphologies than a PSF. In addition, we expect for extended sources the change in their size to be smaller than for a PSF. The CRR image reveals a 28.3% brighter center and a 16.6% smaller spatial size on average based on these three cases, showing that the method is working as we expect from the simulation tests in Section 3.2.

Finally, we check the full spectra of the CRR and Shepard’s image results for real data. Figure 19 shows the ratio of the spectrum from CRR to the spectrum from Shepard’s image. We consider the full spectrum summed over all pixels, and the intensity of the central pixel, as labeled, for plate-IFU 8720-1901. The thick solid curve overlaid on each spectral ratio is the running mean for the ratio smoothing over 100 wavelength slices or 84Å. For the sum of all pixels of the running mean, the ratio of Shepard’s image and our reconstruction is 0.9729 on average with 0.0033 for the standard deviation, meaning it is almost constant in the running mean. The central intensity of our reconstruction is 29.9% brighter than Shepard’s result, with a standard deviation of 0.0439. These results indicate that we have not introduced major wavelength dependent artifacts relative to what might exist in Shepard’s method.

5 Discussion

Shepard’s method is one of the most commonly used techniques to interpolate irregularly sampled data onto a regular grid (see Franke et al 1980; Dell’Accio et al 2016). The largest IFU surveys to date are all using a variant of it. MaNGA and CALIFA both used the form of Shepard’s method described here, and SAMI uses a differently motivated technique, but one which is very similar to Shepard’s method with a different kernel.

However, Shepard’s method produces images whose flux errors are correlated with one another; i.e. they have a very non-diagonal covariance matrix (Law et al 2016). In addition, as we found in this investigation, it unnecessarily broadens the PSF of the resulting image.

Other techniques exist in the statistical literature for interpolating from irregularly distribution samples, for example radial basis function techniques, Wiener interpolation, “kriging,” and Gaussian processes (Krige 1951; Wiener 1964; Schaback 1995; Hartkamp et al 1999; Press 2007; Rasmussen et al 2006). However, these methods are not designed to produce a consistent and tight PSF, and generally lead to highly off-diagonal covariances.

Motivated by a desire to avoid off-diagonal covariances, we examined the techniques of Bolton et al (2010), and in this paper have adapted them to the imaging context. The result is a technique which successfully reduces the off-diagonal covariances to a very small level and also provides a final image PSF which is better than other methods by around 16% in the FWHM, for the MaNGA example we consider.

The off-diagonal covariances produced by other methods can heavily affect the subsequent analysis. For simple measurements such as aperture fluxes, the correct propagation of errors becomes cumbersome and complex. For more complicated measurements such as maximum likelihood model fitting, even determining the best fit parameters depends on accounting for the covariance accurately. The MaNGA data, and IFU data generally, is often used for such measurements. Our reconstruction represents a way to simplify these measurements up front.

One may wonder why it is possible for us to obtain near-zero off-diagonal covariances, considering that dithered samples are spaced every 1.44” while our pixel scale is 0.75”. First, it is worth noting that the covariance between adjacent neighbors is not entirely eliminated, but remains at the few percent level for this pixel size. Second, the differences in seeing and even slight differences in coordinates for different exposures allow the method to create nearly independent data points in neighboring pixels.

The improved image PSF is substantial. It is equivalent to building an instrument with 30–40% more fibers. It provides a greater ability to resolve structures within galaxies and to measure gradients accurately, as well as to find fainter point-like features (e.g. AGN).

A critical aspect of obtaining diagonal covariances is that some of the weights are negative, whereas Shepard’s method has all positive weights (see Figure 3). All positive weights will always produce correlated errors. Although at first glance the fact that a positive fiber flux can contribute negatively to a pixel flux may be non-intuitive, this same feature exists in all accurate image interpolation techniques on a grid (e.g. those based on a sinc kernel) and is not in itself a cause for concern.

We plan to apply this method to the entire MaNGA sample and to test the MaNGA Data Analysis Pipeline code on the revised cubes. Since our method is still experimental (we have only applied it to around 100 or so cubes) we cannot determine yet if it will yield a practical improvement to the results, but we view its prospects as promising.

So far as we can determine, although our technique was derived from that of Bolton et al (2010), this application is entirely new, and represents a new method for scattered image interpolation. Since its main function is to control the covariance matrix of the result, rather than to regularize the result by enforcing a notion of smoothness, we refer to it as covariance-regularized reconstruction. It is relevant when the input samples are noisy and irregular, when result of the interpolation is meant to be a specific grid of values, and when there is a natural resolution (in our case, the kernel) in the sampled image that is meant to be preserved. These conditions apply in other IFU data sets, as well as to a number of ground-based and space-based imaging data sets. Thus, this method may provide an alternative and improved method to analyze those data sets.

Acknowledgements

We thank Yacine Ali-Haïmoud, Stephen Bailey, Adam Bolton, Cristina Mondino, and David W. Hogg for useful discussions regarding this work. MRB was supported in part by National Science Foundation grant NSF-AST-1615997.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.