Estimating dust attenuation from galactic spectra. I. methodology and tests

We develop a method to estimate the dust attenuation curve of galaxies from full spectral fitting of their optical spectra. Motivated from previous studies, we separate the small-scale features from the large-scale spectral shape, by performing a moving average method to both the observed spectrum and the simple stellar population model spectra. The intrinsic dust-free model spectrum is then derived by fitting the observed ratio of the small-scale to large-scale (S/L) components with the S/L ratios of the SSP models. The selective dust attenuation curve is then determined by comparing the observed spectrum with the dust-free model spectrum. One important advantage of this method is that the estimated dust attenuation curve is independent of the shape of theoretical dust attenuation curves. We have done a series of tests on a set of mock spectra covering wide ranges of stellar age and metallicity. We show that our method is able to recover the input dust attenuation curve accurately, although the accuracy depends slightly on signal-to-noise ratio of the spectra. We have applied our method to a number of edge-on galaxies with obvious dust lanes from the ongoing MaNGA survey, deriving their dust attenuation curves and $E(B-V)$ maps, as well as dust-free images in $g$, $r$, and $i$ bands. These galaxies show obvious dust lane features in their original images, which largely disappear after we have corrected the effect of dust extinction. The vertical brightness profiles of these galaxies become axis-symmetric and can well be fitted by a simple model proposed for the disk vertical structure.


INTRODUCTION
The observed spectrum of a galaxy is a combination of several components: a continuum, absorption and emission lines. The continuum and absorption lines are both dominated by starlight, thus usually referred to as the stellar component of the spectrum, while the emissionline component is produced in Hii regions around hot stars, or emission-line regions of active nuclei, or both. All these components, however, are modified by the attenuation of dust grains distributed in the inter-stellar space. Dust attenuation can affect galaxy spectra over a wide range of wavelengths, from ultraviolet (UV), optical to infrared, by absorbing short-wavelength photons in UV/optical and re-emitting photons in the infrared, and the absorption is stronger in shorter wavelength. Consequently, dust attenuation can cause changes in the overall shape of a galaxy spectrum. Such attenuation has to be taken into account before one can measure the different components of an observed spectrum reliably.
Various schemes have been used to measure dust attenuation. Observations of individual stars can directly probe dust attenuation along lines of sight. But this method is limited to Milky Way and very nearby galaxies, such as the Magellanic Clouds (e.g., Cardelli et al. 1989;Fitzpatrick 1999;Gordon et al. 2003). Far-infrared (FIR) observations provide direct measurements of dust attenuation for more distant galaxies because the emission by dust dominates the spectral energy distribution (SED) of FIR, although such observations can be made only through space telescopes. When both infrared (IR) and UV photometry are available, dust attenua-tion may be estimated by the IR-to-UV luminosity ratio, L IR /L UV , known as IRX, as described in (e.g., Meurer et al. 1999;Gordon et al. 2000). In the absence of IR data, the slope of the UV continuum spectrum β can also be used as an alternative estimator (e.g., Calzetti et al. 1994;Meurer et al. 1999). When multi-band photometry covering a wide range of wavelengths is available, a commonly used method is to fit the observed SED with stellar population synthesis models, which provide estimates of a variety of stellar population parameters including dust attenuation (e.g., Noll et al. 2009;Johnson et al. 2013;Chevallard & Charlot 2016). The spectrum of a galaxy should, in principle, contain more detailed information about the physical properties of the galaxy. For star-forming regions which are usually dusty, the attenuation of various lines is commonly estimated from the Balmer decrement, by comparing the observed H α to H β line ratio (H α /H β ) with the intrinsic one predicted by atomic physics applied to a given environment (Osterbrock & Ferland 2006). Dust attenuation in the stellar component, however, cannot be estimated in the same way due to at least two factors. First, the attenuation by star-forming regions and that by stars can be quite different, with the former expected to be stronger in most cases. Second, the Balmer decrement is not measurable in non-star forming regions where emission lines are weak. Therefore, dust attenuation in the stellar component is usually estimated through full spectral fitting. For instance, a simple approach is to match a dust-reddened galactic spectrum with its un-reddened counterpart that is produce by a similar stellar population (e.g., Calzetti et al. 2000;Wild et al. 2011;Reddy et al. 2015;Battisti et al. 2017a,b). Alternatively and more commonly, the estimate of the stellar dust attenuation is made by fitting the full spectrum with a stellar population synthesis model. A number of spectral synthesis fitting codes are now publicly available and are commonly adopted for such purposes. These include PPXF developed by Cappellari & Emsellem 2004;Cappellari 2017, STARLIGHT by Cid Fernandes et al. 2005, STECKMAP by Ocvirk et al. 2006a,b, andVESPA by Tojeiro et al. 2007. In the conventional fitting method, one usually makes an assumption about the shape of dust attenuation curve and treats the amount of attenuation as an free parameter to be obtained from the fitting. Obviously, the fitting result will depend on the theoretical dust attenuation curve chosen. In addition, it is known that dust attenuation has significant degeneracy with stellar age in such fitting.
In order to overcome some of the problems in the conventional method, Wilkinson et al. (2015Wilkinson et al. ( , 2017 have developed a code of full spectral fitting (Firefly), us-ing a new approach to break the degeneracy between dust and other stellar population properties. In this method, dust attenuation is assumed to affect the largescale shape of an observed spectrum but little the spectral shapes on small scales. Their solution is to first apply a high-pass filter (HPF) to remove the large-scale features in both the observed and model spectra, before fitting the filtered observed spectrum with the filtered model spectra.
In this paper we adopt a method similar to that of Wilkinson et al. (2015Wilkinson et al. ( , 2017 to fit the spatially resolved spectra of edge-on disk galaxies selected from SDSS IV MaNGA (Bundy et al. 2015;Blanton et al. 2017) to investigate the dust extinction maps of these galaxies. To do this, we separate small-scale spectral features from large-scale variations using a moving box average rather than the Fourier transform adopted in Wilkinson et al. (2015Wilkinson et al. ( , 2017. Since at a given wavelength dust attenuation affects both the small-scale and large-scale components in a similar way, the ratio between them is expected to be independent of dust extinction, and so can be used to constrain the underlying stellar population. The paper is organized as follows. §2 describes our method. In §3, we test the performance of our method with mock spectra. We then apply our method to the MaNGA data in §4 and compare our results with those obtained earlier. Finally, we summarize and discuss in §5.

THE METHOD
A simple method to extract information from an observed galaxy spectrum is to fit it to a linear combination of simple stellar populations (SSPs). Each SSP is a single, coeval population of stars with a given metallicity and abundance pattern. The spectral synthesis of a SSP, therefore, consists of three components: the evolution of individual stars in the form of isochrones; a library of stellar spectra; and an initial mass function (IMF). Mathematically, the spectrum of a SSP of metallicity Z at the age t can be written as (1) where f star is the spectrum of a star of age t, metallicity Z and initial mass M in the spectral library, and Φ(M ) is the IMF (e.g., Conroy 2013). The dependence of the effective temperature, T eff , and the surface gravity, g, on the initial stellar mass for given t and Z is determined by the stellar evolution model adopted. The lower limit of integration, m 0 , is typically taken to be the hydrogen burning limit, 0.08M , and the upper limit, m(t), is the mass of the most massive stars that can survive to the age t, as determined by the stellar evolution model.
There are several popular stellar population synthesis codes available (e.g., Leitherer et al. 1999;Bruzual & Charlot 2003;Maraston 2005;Vazdekis et al. 2010). In this paper, we use the simple stellar populations given by Bruzual & Charlot (2003, hereafter BC03). BC03 provides a large sample of SSPs, covering 221 ages from t = 0 years to t = 2.0 × 10 10 years, and six metallicities from Z = 0.0001 to Z = 0.05 (note that the solar metallicity Z = 0.02) at a spectral resolution of 3Å. A total of 1326 SSPs are provided by BC03.
Once we have a series of SSPs with different ages and metallicities, the observed spectrum, F obs (λ), can be fitted with a linear combination of the SSPs together with a model of dust extinction: where N * is the number of templates (SSPs in our case), f j SSP is the spectrum of the j th SSP, x j is the weight of the j th SSP, and A(λ) is the dust attenuation curve. So defined, F a (λ) is the model spectrum that takes into account dust extinction, and F (λ) is the dust free model spectrum.

Spectral decomposition
In our analysis, we first decompose a spectrum into two components, one small-scale component, S, and one large-scale component, L. Roughly speaking, the L component can be considered as the continuum shape of the spectrum, and the S component as the composition of absorption and emission line features and noise. Wilkinson et al. (2015Wilkinson et al. ( , 2017 adopted a similar spectral decomposition, using Fourier transforms to separate smalland large-scale components. In our analysis, we adopt a moving average method. Specifically, the large-scale component is obtained through where ∆λ is the size of the wavelength window. The small-scale component is simply defined as Using equation (2), we can write where an approximation is made in the second line that ∆λ is small and A(λ) is smooth so that A varies little over ∆λ. It is then easy to see that where F L (λ) is the large-scale component of the intrinsic spectrum F (λ). Similarly, equation (4) can be written as where F S (λ) is the small-scale component of the intrinsic spectrum. From equations (6) and (7), we see that the ratio between the small-scale and large-scale components is dust free: where R(λ) is the ratio of the intrinsic spectrum. It must be realized that the above relations are derived under the assumption that all the stellar populations have the same extinction given by A(λ). In reality, dust extinction may be more complicated and, in particular, may be different for different stellar populations. Along this line, Charlot & Fall (2000) proposed a twocomponent dust model, in which old stars are assumed to be mixed uniformly with a diffuse dust distribution while young stars are assumed to be embedded in birth clouds of stars that have larger dust optical depth than the diffuse component. This model is motivated by the empirical relation between stellar and nebular dust extinctions (Calzetti et al. 2000).
More generally, suppose that the j th stellar population (not necessarily a SSP) has an extinction curve of A j (λ), defined by where τ j (λ) is the dust optical depth, which is different for different stellar populations. The ratio between the small-scale and large-scale components can then be written as where F j S and F j L are the small-scale and large-scale spectrum of the j th population. This equation reduces to equation (8), as expected, if τ j (λ) is the same for all 'j' (which is the case if the dust distribution is a screen in front of the galaxy, as assumed above), or the intrinsic spectra of different stellar populations are the same. Consider another simple case in which τ j (λ) 1. In this case, one can write where and .
Note that both τ S (λ) and τ L (λ) are weighted averages of τ j (λ). We have estimated T S (λ)/T L (λ) using realistic spectra and extinction curves, and found the ratio to depend only weakly on λ for τ j < 1, so that R a (λ) is again independent of dust attenuation. However, for τ j > 1, the ratio R a (λ) can be affected significantly by the differences in dust attenuation among different stellar populations. We thus conclude that, R a (λ) can well reflect the properties of the underlying intrinsic spectrum under certain assumptions. Throughout this paper we assume R a (λ) ≈ R(λ). We will come back and consider more general cases in the future.

The fitting procedure
In our modeling of galaxy spectra, we fit the ratio R a (λ) obtained from an observed spectrum with the corresponding ratio of the model spectrum. To this end, we decompose each model template to obtain its small-scale and large-scale components as defined by equations (3) and (4). For a given set of fitting coefficients, {x j }, we obtain the model prediction of the ratio, where N * is the number of templates, f j S (λ) and f j L (λ) are the small-scale and large-scale components of the j th template, f j (λ), respectively. In the equation, the constant f 0 is the normalization of the model template, and is the same for all the f j S (λ) and f j L (λ). The predicted ratio R m (λ) is then compared with the observed ratio R a (λ) to determine {x j }. The fitting is carried out by using MPFIT, which is a non-linear Least-squares Fitting code in IDL (Markwardt 2009). We should point out that, the fact that the normalization factor f 0 is cancelled out in equation (14) indicates that this factor cannot be directly determined from this fitting procedure. We will come back to this point later.
In principle, one can use all the SSPs from the BC03 library, or a random subset, as the model templates for the above fitting process. SSPs with different ages and metallicities have different small-scale features. Similarly, the small-scale to large-scale (S/L) component ratios of different SSPs also have different features. This is shown in Figure 1, where the first 10 SSPs have different ages but the same metallicity (solar metallicity), while the last 5 have the same age (5 Gyr) but different metallicities. The left-hand panel shows the original SSPs and the right-hand panel shows the corresponding S/L ratio spectra. It is obvious that different SSPs have different absorption line strengths and other broader features. It is these differences that allow us to derive constraints on the stellar populations from a galaxy spectrum.
In order to speed up the fitting process, we construct a small set of model templates by applying the technique of Principal Component Analysis (PCA, Deeming 1964) to the BC03 SSP spectra. PCA can effectively reduce the size of the template library, as the principal features of the library can be described by the first few eigenspectra. For instance, Li et al. (2005) obtained galactic eigen-spectra using PCA and showed that the first nine eigen-spectra already provided the base to model the stellar spectra of the galaxies in Sloan Digital Sky Survey Data Release 1 (SDSS DR1, Abazajian et al. 2003). In our analysis, we apply the PCA to the 1326 BC03 SSPs, and we adopt the first 14 eigen-spectra resulted from the PCA as the model templates for our spectral fitting. The cumulative contribution to variance by the first 14 eigen-spectra is 99.985%, indicating that they contain nearly all the information of the original 1326 SSPs. Figure 2 displays the 14 eigen-spectra (the left panel), as well as the corresponding S/L ratio spectra (the right panel). We should point out that the number of eigen-spectra is determined so as to include as much as possible information of the original SSPs, and simultaneously to reduce as much as possible the computing time. We have repeated our tests to be presented below, by adopting a different number of eigen-spectra, and found that the spectral fitting is little affected as long as the majority of the original information of the whole SSP library (> 99%) is contained in the adopted eigen-spectra.
Once the coefficients, x j , are obtained from the fitting, the best fitting spectrum, F fit , which is expected to not contain the effect of dust extinction, can be reconstructed by Note that the normalization factor f 0 is unknown, and it is set to be arbitrary for the moment. This means that the fitting procedure described above gives only the shape of the dust-free spectrum. The shape of the dust attenuation curve can then be obtained by comparing the best fitting spectrum with the observed spectrum: Conventionally, A(λ) is written as where A(λ) is the dust attenuation at the wavelength of λ. In practice, we normalize both F a (λ) and F fit (λ) at a given wavelength, and we have where A V is the dust attenuation at λ V , and are the dust-free and the observed spectrum normalized at λ V , respectively. Note that the factor f 0 is implicitly contained in both F fit (λ) and F fit (λ V ), and so is cancelled out due to the normalization. Therefore, what we obtain from this fitting procedure is the relative dust attenuation curve, i.e. A(λ) − A V , which is the dust attenuation as a function of wavelength relative to the attenuation at λ V . Equation (19) shows that our method provides a direct measurement of the relative dust attenuation curve, with no need to assume a functional form for the curve. In practice, however, it is may still be desirable to use  Wavelength (Å) RMSE=0.008 Figure 3. Examples of the fitting procedure for measuring relative dust attenuation curves and E(B − V ). Panels from left to right correspond to three mock spectra with different signal-to-noise ratios (SNRs): 5, 10 and 20. The mock spectra are constructed from a same intrinsic dust-free spectrum, and are reddened with a Calzetti extinction curve assuming a color excess of E(B − V ) = 0.2. The top panels display the mock spectra, and panels in the second row show the small-scale and large-scale components of each spectrum. In the third row, the black line in each panel is the ratio of small-scale to large-scale component (S/L), while and the red line plots the S/L of the best-fitting model spectrum. In the fourth row, the black line is the recovered relative attenuation curve and the red line is the input attenuation curve. The E(B − V ) estimated from the measured attenuation curves is 0.221, 0.194, and 0.205 respectively, as indicated. In the bottom panels, the black lines show the ratios of the measured attenuation curve (smoothed for clarity) to the input one. The root square mean error (RMSE) of these curves is 2.9%, 1.7% and 0.8%, also indicated.
a parametric form to represent the curve. As shown below, our measurements of A(λ) − A(V ) can well be represented by a second-order polynomial in most (if not all) cases, The parameters a 1 and a 2 are to be obtained by fitting the above function to the A(λ) − A V measured by using equation (19). Given the measured A(λ) − A V and setting λ = λ B = 4400Å, we therefore can get the (B − V ) color excess where A B is the dust attenuation at the wavelength, λ B , corresponding to the B-band. Similarly, the selective at- can also be obtained. In practice, dust extinction is sometimes described by the total attenuation curve, defined as where R V is the value of the total attenuation curve in the V -band, R V is also known as the ratio of the total to selective extinction in the V -band. This equation indicates that, if one were to obtain the total (absolute) dust attenuation at given λ, it is necessary to know the value of R V (or more generally speaking, the total attenuation at any given wavelength which falls in the wavelength range considered). For the Milky Way, the typical value of R V is 3.1, but it can vary between 2 and 6 (Cardelli et al. 1989). The R V of the Calzetti dust curve is 4.05 with scatter of 0.8 (Calzetti et al. 2000). Figure 3 shows some examples to demonstrate the step-by-step application of our method to measure the relative dust attenuation curve, A(λ) − A V , as well as the color excess, E(B − V ). The three panels in the first row display three mock spectra that correspond to E(B − V ) = 0.2 but have different Gaussian signal-tonoise ratios SNR=5, 10, and 20, as indicated on the top of the figure. Here a Calzetti dust curve (Calzetti et al. 2000) is adopted to create the mock spectra, and the procedure to construct the mock spectra will be detailed in the next section. In the second row, the moving average filter (see Eqn. 3 and 4) is applied to decompose  each spectrum into the large-scale and small-scale components, which are plotted separately in each panel. The third row shows the small-to large-scale spectral ratio in black lines, and the best-fitting ratio (see Eqn. 14) in red lines, as obtained from the non-linear fitting described above. The coefficients obtained from the fitting are then used to derive the best-fitting model spectrum (see Eqn. 15), which is expected to be dust-free. The relative dust attenuation curves, given by the ratio between the original spectra and the best-fitting models (see Eqn. 19), are shown as the black lines in the fourth row. The input attenuation curve is repeated in the three panels as a red line for comparison. The E(B − V ) values obtained from our measured attenuation curves (see Eqn. 22) are 0.221, 0.194 and 0.205 for the three SNRs, respectively, as indicated in each panel, all very close to the input value, E(B − V ) = 0.2. In the bottom panels, we show the ratios of the measured attenuation curve to the input attenuation curve, which are smoothed for clarity. As one can see, our method recovers the input attenuation curve very well. The rms deviations of the measured attenuation curves around the input one are 2.9%, 1.7% and 0.8% for the three SNRs, respectively.

Examples
3. TEST WITH MOCK SPECTRA

The mock spectra
In order to test the reliability of our method, we have created a series of synthetic spectra from the BC03 SSPs. We first pick 150 SSPs out of the 1326 BC03 SSPs, with 25 different ages and 6 different metallicities. For each metallicity, the 25 SSPs are chosen so as to cover the full range of age from 0.001 Gyr to 18 Gyr, with approximately equal intervals in the logarithm of the age. We then randomly select one of the 25 SSPs for every metallicity, and the six selected spectra of different metallicites and ages are then normalized at 5500Å and combined with random weights. We repeat this step 1000 times, thus creating 1000 synthetic mock spectra with a wide coverage of age and metallicity, as shown in Fig. 4. In the figure, the age-metallicity space is roughly divided into four regions as indicated by the vertical and horizontal lines: (a) young age and high metallicity, (b) old age and high metallicity, (c) young age and low metallicity, (d) old age and low metallicity. Each of the mock spectra is reddened with a Calzetti dust curve (Calzetti et al. 2000), but assuming four different color excesses: E(B − V ) = 0.05, 0.15, 0.3, and 0.5. Finally, a Gaussian noise with SNR=5, 10, 20, or 30 is added to each spectrum. This procedure results in a total of 16000 mock spectra, which are used to test our method. We would like to point out that some of the mock spectra may not be physically meaningful, e.g. those with extremely old ages but high metallicities and those with young ages but low metallicities. We do not exclude these spectra from our test, because the inclusion of them does not affect our test results, as we will show below.

The choice of eigen-spectra number
As mentioned, we adopt the first 14 eigen-spectra resulted from the PCA analysis as the model templates for our spectral fitting. With the mock spectra generated above, we have examined the potential effect of using different numbers of eigen-spectra on the measurement of E(B − V ). For a given number of eigen-spectra and the mock spectra with a given SNR, we calculate the standard deviation, σ, of the difference between the output and input E(B − V ) values, and we use this quantity to indicate the goodness of our method in recovering the true E(B − V ). The result of this analysis is shown in Fig. 5, which plots σ as function of the number of eigen-spectra for four different SNRs.
The figure shows that, for mock spectra with SNR> 10, the value of σ decreases as one adopts more and more eigen-spectra in the fitting, but σ does not decrease any more when the number of eigen-spectra exceeds ∼ 15. The value of σ keeps roughly constant for SNR=10 and even increases for SNR=5. In spectra with low SNRs, features on small scales such as stellar absorption lines can be well dominated by noise, and so the use of too many eigen-spectra actually leads to poor fits. Overall, the figure indicates that the use of 8-15 eigen-spectra appears to be able to achieve a compromise among different SNRs. We opt for the first 14 eigen-spectra considering mainly the fact that the sigma starts to increase with 15 eigen-spectra at all SNRs. Nevertheless, we note that the variation of σ with the number of eigen-spectra used is in general modest in comparison to that produced by different SNRs.

The choice of filtering window size
As mentioned before, filtering to separate a spectrum into small-and large-scale components plays an important role in our method, and it is critical to find a suitable choice of the wavelength window size for the filtering. For this purpose, we have adopted window sizes ranging from 160Å to 720Å when filtering the mock spectra. Fig. 6 shows σ versus the wavelength window size, for mock spectra of different E(B − V ) and SNRs. Defined in the same way as above, the σ is an indicator of the goodness of our method. The four colorful curves in each panel correspond to the four regions divided in Fig. 4, while the black curve is for all the spectra with the given E(B − V ) and SNR. Generally, as expected, the value of σ decreases as SNR increases. In addition, for low values of SNR (5 and 10), σ tends to decrease as the wavelength window size increases, but the trend reverses when SNR and E(B − V ) are high.
The results shown in Fig. 6 suggest that the wavelength window size should be chosen to be about 500Å     for SNR smaller than 10 and about 300Å for higher SNR. This dependence of smoothing window size on SNR is understandable. A larger window size can reduce the effects of noise in the low SNR spectra more effectively, while a smaller window size is needed to retain more real features in the spectra with high SNR. Spectra with older stellar ages (green and yellow lines) also trend to have lower σ than younger ones (red and blue lines), indicating that our method works better for older stellar populations. However, the differences in the results among the four different regions are not large in comparison with the effects of the SNR and the wavelength window size. We have chosen a wavelength window size of 300Å for all the cases, according to the analysis presented in the previous subsection. For low-SNR spectra with SNR=5 and 10, we repeat the test using a larger window size of 500Å, and plot the results as the orange histograms in the top two rows. As can be seen, the results change little with varying wavelength window sizes.

Test results
The standard deviation, σ, decreases as SNR increases, but its dependence on the input E(B − V ) is quite weak. The distribution of ∆E(B−V ) in most cases is roughly a Gaussian centered at ∆E(B − V ) ∼ 0, indicating that the model inference is unbiased. The only exception is the case shown in the first panel, where the distribution is skewed to positive values of ∆E(B − V ). In this case, the inputs of E(B − V ) and SNR are both small. We note that increasing the wavelength window size for the moving average filter does not make a significant improvement of the E(B − V ) measurement in such low-SNR and low-E(B − V ) cases. For reference, the mean and σ of ∆E(B − V ) are listed in Table 1 for mock spectra generated with a number of SNRs.
We also show the results for "unreasonable" mock spectra as the red histograms in individual panels. These spectra are located in regions b and c in Fig. 4: they either have old ages and high metallicities, or young ages and low metallicities. The distributions of ∆E(B − V ) for these spectra are almost the same as those of the full samples at given E(B − V ) and SNR, indicating that the inclusion of such spectra does not introduce biases into our test results.
As described in §2.2, our method can recover the relative attenuation curve. As a demonstration, the upper row of Fig. 8 shows the median A λ − A V curve at given E(B − V ) and SNR. Panels from left to right are for different SNR bins, and the different colors in each panel are for different E(B − V ) bins. For each case the corresponding Calzetti curve is plotted as a dotted line. The differences between the median curves output from our method and the input Calzetti curve are shown in the lower row. As one can see, our method indeed recovers the input curve very well in all cases and at all wavelengths. The largest difference is ∼ 0.03 magnitude, and this occurs at lowest SNR (SNR=5) and at both short (< 5000Å) and long wavelengths (> 6500Å). At higher SNRs the difference is generally very small, less than 0.01 magnitude at all wavelengths. There is a slight but noticeable jump at ∼ 6200Å in all cases. This is caused by the fact that the Calzetti curve is a piecewise function with two sub-functions separated at a fixed wavelength, while the attenuation curves output from our method are smooth, thus not necessarily following any predefined functional forms.
To conclude, the results of the test demonstrate that the input E(B − V ) are well recovered by our method for various cases. The systematic bias in ∆E(B − V ) is quite small, indicating that potential uncertainties, such as that produced by the dust-age degeneracy, is well taken care of by our method.

The MaNGA data
Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey (Bundy et al. 2015;Yan et al. 2016a,b;Law et al. 2016;Wake et al. 2017) is one of the three core programs of SDSS-IV (Sloan Digital Sky Survey IV, Blanton et al. 2017). MaNGA aims to obtain spatially resolved spectroscopy of about 10,000 nearby galaxies. In the outskirts of galaxies, MaNGA can reach a SNR of 4-8 (perÅ per 2 fiber) at 23 AB mag arcsec −2 in the r-band. The wavelength coverage is between 3,600 and 10,300Å with a spectral resolution R ∼ 2000 (Drory et al. 2015). In this paper, we select galaxy samples from the Sloan Digital Sky Survey Data Release 14 1 (SDSS DR14, Abolfathi et al. 2018), which includes integral field spectroscopy (IFS) data from MaNGA for 2812 galaxies, including ancillary targets and ∼ 50 repeated observations. We visually inspect the SDSS gri colorcomposite images of all the galaxies in this release, and select 15 edge-on galaxies that have obvious dust lane features. Galaxies of this kind are much more affected by dust attenuation than those viewed from face-on, and thus provide more demanding tests of our method. In the rest of this paper, we will apply our method to the 15 galaxies to demonstrate that our method can effectively measure the spatially-resolved E(B−V ) maps and relative dust attenuation curves, even for these difficult cases.

Dust extinction maps and profiles
In the current version of our code, velocity dispersion of the observed spectrum needs to be specified so that model templates can be adapted to the observed spectral resolution. We first use pPXF (Cappellari & Emsellem 2004;Cappellari 2017) to measure stellar kinematics and correct the broadening effect of velocity dispersion. We then apply our method to all the spatial pixels in the MaNGA datacubes whose continuum SNRs are greater than 5. Fig. 9 displays the two-dimensional maps obtained from our method for the 15 galaxies in our sample. The first panel in each row shows the SDSS gri color-composite image of the galaxy in question, and the second panel is the E(B − V ) map obtained from our method. As one can see, the dust lanes in the original images are well reproduced in the E(B − V ) maps, in terms of both the spatial distribution and the relative strength.
Once the E(B − V ) map is obtained for a galaxy, we can correct the effect of dust extinction on the gri image. Since the FWHM (full width at half maximum) of the point-spread function (PSF) for MaNGA is ∼ 2.5 , while that of the SDSS gri image is ∼ 1.4 , we convolve the SDSS image in g, r and i bands with a Gaussian kernel to match the MaNGA PSF, before correcting the dust extinction in the images according to the measured attenuation curve. The dust-corrected image in the three bands are then combined to form a new gri color-composite image, which is expected to be dust free. This corrected image is shown in the third panel for each galaxy. The dust lanes in the original images are no longer seen in the dust-corrected images. This can be seen more clearly in the fourth panel of each row, where the red and blue curves show, respectively, the original and corrected brightness profiles in the gband, as measured along the arrow indicated in the third panel. Compared to the corrected brightness profile, the original one is not only lower in amplitude, but also asymmetric in shape due to the strong extinction in dust lane regions. In almost all cases, the prominent features . There are five panels in each row. The first one is SDSS gri color-composite image with the MaNGA footprint in magenta. The second one is the E(B − V ) map derived from our method. The third one is dust corrected gri color-composite image. The red cross shows the center of the integral field. In the fourth panel, the red line is reddened brightness profile while the blue line is un-reddened brightness profile along the direction of arrow in the third panel. The un-reddened brightness profile is fitted by a model that the black dash line shows. n is the key parameter of this model. The last panel shows the selective attenuation curves measured from our method. There are five panels in each row. The first one is SDSS gri color-composite image with the MaNGA footprint in magenta. The second one is the E(B − V ) map derived from our method. The third one is dust corrected gri colorcomposite image. The red cross shows the center of the integral field. In the fourth panel, the red line is reddened brightness profile while the blue line is un-reddened brightness profile along the direction of arrow in the third panel. The un-reddened brightness profile is fitted by a model that the black dash line shows. n is the key parameter of this model. The last panel shows the selective attenuation curves measured from our method. in the vertical surface brightness profiles produced by the dust absorption are no longer present in the dustcorrected profiles. The corrected profiles are roughly symmetric with respect to the peak, as is expected for axis-symmetric thin disks. We attempt to fit the dust corrected brightness profiles with a vertical brightness profile model. Since galaxy disks are not infinitesimally thin, the threedimensional luminosity density of the disk is typically written as The first part of this equation describes the surface brightness as a function of the radius R, which is assumed to have an exponential profile, with R d being the scale-length of the disk. The function f (z) describes the surface brightness distribution in the vertical (z) direction (see §2.3.3 of Mo et al. (2010) for a review). A commonly adopted fitting function of f (z) is where n is a parameter controlling the shape of profile near z = 0, and z d is the scale-height of the disk. In particular, n = 1 corresponds to a self-gravitating isothermal sheet while n = ∞ corresponds to an exponential profile. Note that f n (z) decreases exponentially at large |z|, but a larger value of n gives a steeper profile near the mid-plane. When fitting the observed profile, the effect of seeing on the observed data must be included. This is done by convolving the model profile with a Gaussian kernel that matches the MaNGA PSF. The best-fit model profile is plotted for each galaxy as the black dash line in the fourth panel of Fig. 9. For reference, we also indicate the values of n in the corresponding panels. It is clear that the corrected profiles are well described by the model. Fig. 10 shows the distribution of 2/n for the 15 galaxies in comparison with that obtained by de Grijs et al. (1997) from fitting the K-band images of 24 nearly edgeon spiral galaxies. The two distributions are qualitatively consistent with each other, suggesting again that our method is able to correct dust extinctions for nearly edge-on galaxies. Unfortunately, no quantitative comparison can be made, as the two samples are small and have different selections.

Dust attenuation curves
As described in §2.2, one important advantage of our method is that the relative dust attenuation as a function of wavelength (i.e. the relative attenuation curve) can be directly obtained without the assumption of a functional form for the shape of the attenuation curve or the adoption of a theoretical dust model. The selective attenuation curves measured from our method are plotted in the rightmost panel of each row in Fig. 9. The grey region covers the range spanned by all the individual spaxels in a given galaxy. The blue dashed line is the median, while the pink region shows the standard deviation of the spaxels around the median. For comparison, the commonly-adopted dust curves of Calzetti et al. (2000) and Cardelli et al. (1989) (CCM) are plotted as the black dashed and yellow dashed lines, respectively. Over the wavelength range probed here, the median dust curves measured with our method lie between the Calzetti and CCM curves for almost all the sample galaxies. However, there are spaxels where the recovered dust curves deviate significantly from both empirical curves. The deviation occurs largely at λ > 6, 000Å, where we see large scatter among the individual spaxels. This indicates that the dust curves at long wavelengths may vary from region to region, and probably also from galaxy to galaxy.
In order to better understand the variation of the attenuation curves, we further examine the measured selective attenuation curves in regions of different E(B − V ) and those located at different radii from the galactic center. The results are plotted in Fig. 11 where we show the difference of the median of our selective attenuation curves relative to either the Calzetti curve (upper pan- Wavelength (Å) Figure 11. The first row shows the differences between the medians of all relative dust attenuation curves and Calzetti curve in each E(B − V ) and radius bins. 'Re' is effective radius. The second row shows the differences between the median curves and the CCM curve. els) or the CCM curve (lower panels). Panels from left to right are for different ranges of E(B − V ) as indicated. In each panel the different lines are for different radial intervals, also indicated. Overall, the measured attenuation curves deviate from both empirical models to varying degrees, depending on both E(B − V ) and radial distance. The deviation is the largest in the rightmost panels with E(B − V ) > 0.3, where ∆(A λ − A V ) is close to zero at intermediate wavelengths and increases at both shorter and longer wavelengths, up to ∼ 0.15. At wavelengths λ 6000Å, the attenuation at a fixed wavelength increases from the inner to outer regions, with a difference of ∼ 0.1 magnitude between the region within 0.5R e and that beyond R e . These trends are also seen in the left two columns, but the amplitude of the deviation/difference becomes smaller with decreasing E(B − V ). In the leftmost panels where E(B − V ) < 0.15, the difference is quite small, less than ∼ 0.03 magnitude, comparable to the typical offset between the input and output found in the test based on the mock spectra (see Fig. 8). When comparing the two empirical models, we find that the Calzetti curve is closer to our measured curves at λ 5500Å. At longer wavelengths, the Calzetti curve works better at smaller radii, while the CCM curve works better at larger radii. These results suggest that different regions of a galaxy are described by different dust attenuation curves.

SUMMARY AND DISCUSSION
We have developed a method to estimate the relative dust attenuation curve, A(λ) − A V (see equation 19), which is the dust attenuation as a function of wavelength relative to that at the fixed wavelength λ V . For an observed spectrum of a galaxy or a specific region within a galaxy, our method first decomposes the spectrum into a small-scale component, S, and a large-scale component, L, by adopting a moving box average method. Assuming all the stellar populations that contribute to the observed spectrum have the same extinction given by A(λ), we are able to show that the ratio of the two components, R(λ) = F S (λ)/F L (λ), is free of dust attenuation. The observed R(λ) can be modeled with a given theoretical library of simple stellar populations (SSPs), without the need for modeling the effect of dust attenuation in the fitting. The dust-free intrinsic spectrum is then reconstructed from the SSPs according to the corresponding coefficients that form the best-fit to R(λ). Finally, the relative dust attenuation curve is obtained by comparing the observed spectrum with the 'dust free' model spectrum.
We have performed extensive tests of our method on a set of mock spectra covering wide ranges in stellar age and metallicity, as well as E(B − V ) and spectral signalto-noise ratios (SNRs). These tests show that both the E(B −V ) and the relative dust attenuation curve can be well recovered as long as E(B−V ) > 0.05 and SN R > 5. At lower E(B − V ) and smaller SNRs, on average, our method tends to overestimate E(B − V ) by 0.03 magnitude (see Fig. 7). We have applied this method to the integral field spectroscopy (IFS) of 15 edge-on galaxies from the ongoing MaNGA survey, obtaining both the two-dimensional map and radial profile of E(B − V ) for each galaxy, as well as the spatially-resolved relative dust attenuation curve. Using the E(B − V ) maps we have corrected the effect of dust extinction on the SDSS gri image of these galaxies. The dust lanes in the original images become invisible in the images reconstructed from the extinction-corrected spectra, and the vertical brightness profiles in a given band become almost symmetric and can be well described by a simple model proposed for the disk vertical structure.
Although small in number, the MaNGA galaxies already show variations in the slope of dust attenuation curves, which are the most significant at wavelengths longer than ∼ 6000Å and in the outskirt of galaxies. For instance, at ∼ 8000Å and when E(B − V ) exceeds 0.3, A λ − A V increases from the inner region within R < 0.5R e to the outer region (R > R e ) by ∼ 0.1 mag (see Fig. 11). Therefore, in the optical band as covered by the MaNGA spectra, the variation of dust attenuation curve from galaxy to galaxy and from region to region within a galaxy are quite small.
We have compared the estimated dust attenuation curves of the 15 MaNGA galaxies with empirical models of Calzetti et al. (2000) and (Cardelli et al. 1989, CCM), which have been commonly adopted in previous studies respectively for star-burst galaxies and Milky Way-like galaxies. Fig. 9 and 11 show that both models can reasonably well, but not perfectly describe the dust attenuation curves of these galaxies. Overall, the deviation from both models to our data is less than 0.15 mag. In the inner regions of these galaxies where R < 0.5R e , the Calzetti model works better than the CCM model, with ∆(A λ − A V ) < 0.05 at all E(B − V ) values probed and all wavelengths in the optical. In the outer regions beyond the effective radius (R > R e ), the CCM model appears to work better at wavelengths longer than ∼ 5500Å, while both models present relatively large deviation from the data, particularly at large E(B − V ).
Although neither models can perfectly describe the estimated dust attenuation curves, the deviations are nevertheless quite small, less than ∼ 0.1 mag in most cases. Therefore, our conclusion is that both the Calzetti and CCM models can well be adopted for modelling the optical attenuation of galaxies. This result is probably not unexpected given both the high similarity of the different attenuation models in the optical and the relatively weak dependence of the attenuation on wavelength in the narrow optical band. It is known that the various attenuation models differ more significantly in the ultraviolet. In particular, the CCM model presents a bump at around 2175Å, which is absent in the Calzetti curve. Therefore, previous observational studies of dust attenuation curves have mainly relied on broad-band spectral energy distribution (SED) including UV bands (e.g., Salmon et al. 2016;Battisti et al. 2017a,b;Salim et al. 2018;Tress et al. 2018;Decleir et al. 2019;Salim & Boquien 2019). In many cases, these studies have also found evidence in support of variations in the slope of attenuation curves and/or the strength of the UV bump. In a recent theoretical paper (Narayanan et al. 2018), a model for the origin of such variations was developed, in which the variation in the attenuation curve slope depends primarily on complexities in the star-dust geometry, while the bump strength is primarily influenced by the fraction of unobscured O and B stars. In principle, if applied to galactic spectra covering the UV band, our method should be able to provide an independent way of quantifying the variations in the attenuation curve and the UV bump. Next-generation large spectroscopic surveys such as the Prime Focus Spectrograph surveys (PFS; Takada et al. 2014) will obtain high-quality spectra for hundreds of thousands of galaxies at 1 z 2, for which the rest-frame UV bump will be well covered in the observed spectrum. We expect to apply our method to those galaxies, thus directly deriving the selective dust attenuation curves for high-z galaxies.
As mentioned, our method of measuring the relative dust attenuation curves was motivated from Wilkinson et al. (2015), which originally proposed the idea of decomposing an observed spectrum into small-and largescale components before performing full spectral fitting. In our work, we have made a step forward by introducing two important improvements. First, we determine the best-fit stellar spectrum by fitting the ratio of the small to large-scale components, R λ , instead of the filtered spectrum (i.e. the small-scale component only) as in Wilkinson et al. (2015). As can be seen from Eqn. 7, the small-scale component alone is not fully free of dust, although the effect of dust attenuation is expected to be much weaker than the effect on the large-scale component. In contrast, as derived in § 2.1, the small-scale to large-scale component ratio, R λ is essentially dust free (see Eqn. 8) provided that the dust extinction of all the stellar populations in the probed region can be described by a common curve of A(λ). This is actually the assumption commonly-adopted in previous studies, where the dust is treated as a screen in front of the whole galaxy or a given region of the galaxy, so that all the stellar populations in the galaxy/region have the same extinction. In reality this cannot be always true. However, we prove in § 2.1 that our method still works in the case of multiple dust components, as long as the dust optical depth τ (λ) < 1 for all the stellar populations and at all wavelengths considered.
The second improvement is that, for the purpose of decomposing the total spectrum into the small and largescale components, we adopt the moving average filter (Eqn. 3) instead of a high pass filter (HPF) as adopted in Wilkinson et al. (2015). In that work, the decomposition was done by applying the HPF to the Fourier transform of the total spectrum, thus removing the large-scale features through the use of an empirically chosen window function. This approach is not applicable to our method, as we need to obtain the small-scale and largescale components simultaneously so as to have their ratio. Therefore, we apply the moving average filter to obtain the large-scale component in the first place, which is then subtracted from the total spectrum to give the small-scale component. As we have shown in § 3.3, the moving average filtering can be conveniently and reliably implemented for different spectra over wide ranges of SNRs and E(B − V ), as long as a filtering window size of 300-500Å is adopted.