Uncertainties in gas kinematics arising from stellar continuum modelling in integral field spectroscopy data: the case of NGC2906 observed with MUSE/VLT

TL;DR

This study evaluates how different stellar continuum subtraction and line fitting methods influence the derived gas kinematic parameters in integral field spectroscopy data of galaxy NGC 2906, highlighting the importance of method choice and error considerations.

Contribution

It systematically compares multiple stellar subtraction and line fitting approaches, identifying the most impactful factors on gas kinematic measurements in IFU data.

Findings

Stellar subtraction method significantly affects gas kinematic results.

[NII] line measurements are more robust against stellar absorption effects.

Adding a systematic error component accounts for method-related uncertainties.

Abstract

We study how the use of several stellar subtraction methods and line fitting approaches can affect the derivation of the main kinematic parameters (velocity and velocity dispersion fields) of the ionized gas component. The target of this work is the nearby galaxy NGC 2906, observed with the MUSE instrument at Very Large Telescope. A sample of twelve spectra is selected from the inner (nucleus) and outer (spiral arms) regions, characterized by different ionization mechanisms. We compare three different methods to subtract the stellar continuum (FIT3D, STARLIGHT and pPXF), combined with one of the following stellar libraries: MILES, STELIB and GRANADA+MILES. The choice of the stellar subtraction method is the most important ingredient affecting the derivation of the gas kinematics, followed by the choice of the stellar library and by the line fitting approach. In our data, typical uncertainties in the observed wavelength and width of the H\alpha and [NII] lines are of the order of <\delta\lambda>_rms \sim 0.1\AA\ and <\delta\sigma>_rms \sim 0.2\AA\ (\sim 5 and 10km/s, respectively). The results obtained from the [NII] line seem to be slightly more robust, as it is less affected by stellar absorption than H\alpha. All methods considered yield statistically consistent measurements once a mean systemic contribution \Delta\bar\lambda=\Delta\bar\sigma=0.2xDelta_{MUSE} is added in quadrature to the line fitting errors, where \Delta_{MUSE} = 1.1\AA\ \sim 50 km/s denotes the instrumental resolution of the MUSE spectra. Although the subtraction of the stellar continuum is critical in order to recover line fluxes, any method (including none) can be used in order to measure the gas kinematics, as long as an additional component of 0.2 x Delta_MUSE is added to the error budget.

Tables

Table 1. Table 1: Set of SSP models and stellar libraries used for the spectral fitting.

Method	Model	Library	Isochrone	IMF	Metallicities	Nt	Ages [Gyr]	FWHM [Å]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
FIT3D	SED@	GRANADA(+MILES)	Pad00	Salp	0.002, 0.008, 0.02, 0.03	156	0.001–13	2
STARLIGHT	CB07	MILES	Pad94	Chab	0.004, 0.008, 0.02, 0.05	66	0.001–18	2.5
STARLIGHT	BC03	STELIB	Pad94	Chab	0.0001, 0.0004, 0.004, 0.008, 0.02, 0.05	150	0.001–18	3.0
pPXF	V10	MILES	Pad00	Salp	0.0004, 0.001, 0.004, 0.008, 0.019, 0.03	156	0.1–18	2.5

Table 2. Table 2: Different cases considered for the line fitting analysis for the wavelength and for the width.

Cases	Wavelength ($\lambda$)	Width ($\sigma$)
[1]	fixed	fixed
[2]	free	fixed
[3]	fixed	free
[4]	free	free

Table 3. Table 3: Mean values and standard deviations for the wavelength and line width of the H α 𝛼 \alpha and [NII] emission lines obtained when combining all line fitting cases and stellar continuum subtraction methods for each individual spaxel.

Spectrum	Offset	log(NII/H$\alpha$)	${\overline{\lambda }}_{H\alpha }\pm \Delta {\overline{\lambda }}_{H\alpha }$	${\overline{\sigma }}_{H\alpha }\pm \Delta {\overline{\sigma }}_{H\alpha }$	${\overline{\lambda }}_{\left[\mathrm{NII}\right]}\pm \Delta {\overline{\lambda }}_{\left[\mathrm{NII}\right]}$	${\overline{\sigma }}_{\left[\mathrm{NII}\right]}\pm \Delta \overline{\sigma }\left[\mathrm{NII}\right]$
name	(′′, ′′)		[Å]	[Å]	[Å]	[Å]
(1)	(2)	(3)	(4)	(5)	(6)	(7)
a	$(0.0,-0.6)$	0.74	6563.60 $\pm$ 0.32	1.84 $\pm$ 0.39	6584.240 $\pm$ 0.097	1.91 $\pm$ 0.19
b	$(+0.6,0.0)$	0.59	6561.97 $\pm$ 0.15	1.80 $\pm$ 0.27	6582.554 $\pm$ 0.048	1.89 $\pm$ 0.13
c	$(0.0,+0.6)$	0.46	6563.35 $\pm$ 0.23	2.02 $\pm$ 0.29	6584.141 $\pm$ 0.094	1.95 $\pm$ 0.17
d	$(0.0,0.0)$	0.39	6563.69 $\pm$ 0.21	2.02 $\pm$ 0.24	6584.04 $\pm$ 0.12	2.16 $\pm$ 0.12
e	$(-0.6,0.0)$	0.20	6565.249 $\pm$ 0.085	1.55 $\pm$ 0.18	6585.810 $\pm$ 0.051	1.64 $\pm$ 0.12
f	$(-3.0,+11.4)$	–0.06	6564.27 $\pm$ 0.11	1.27 $\pm$ 0.21	6584.761 $\pm$ 0.086	1.03 $\pm$ 0.19
g	$(+4.2,+5.2)$	–0.12	6562.75 $\pm$ 0.15	1.17 $\pm$ 0.11	6582.93 $\pm$ 0.31	1.36 $\pm$ 0.15
h	$(+3.6,+5.2)$	–0.33	6563.224 $\pm$ 0.039	1.239 $\pm$ 0.097	6583.776 $\pm$ 0.047	1.17 $\pm$ 0.15
i	$(-1.6,-7.8)$	–0.49	6563.559 $\pm$ 0.026	1.15 $\pm$ 0.13	6584.197 $\pm$ 0.056	1.11 $\pm$ 0.14
j	$(-4.8,-8.8)$	–0.53	6564.460 $\pm$ 0.031	1.15 $\pm$ 0.11	6585.004 $\pm$ 0.053	1.15 $\pm$ 0.13
k	$(-4.5,+8.6)$	–0.61	6564.893 $\pm$ 0.046	1.242 $\pm$ 0.093	6585.33 $\pm$ 0.16	1.27 $\pm$ 0.17
l	$(-4.4,+7.0)$	–0.62	6565.241 $\pm$ 0.031	1.21 $\pm$ 0.12	6585.51 $\pm$ 0.43	1.37 $\pm$ 0.27

Table 4. Table 4: Mean and standard deviation of the mean kinematic results averaged over all individual spectra.

	Case			Method
		Unsubtracted	FIT3D	STARLIGHT (MILES)	STARLIGHT (STELIB)	pPXF
$⟨\delta {\lambda }_{mc}^{H\alpha }⟩$	[1]	0.023 $\pm$ 0.043	0.000 $\pm$ 0.057	0.004 $\pm$ 0.034	–0.021 $\pm$ 0.036	–0.047 $\pm$ 0.049
$\pm$	[2]	0.091 $\pm$ 0.069	0.041 $\pm$ 0.069	0.053 $\pm$ 0.029	–0.013 $\pm$ 0.039	–0.12 $\pm$ 0.13
$\Delta {\lambda }_{mc}^{H\alpha }$	[3]	0.024 $\pm$ 0.047	0.001 $\pm$ 0.055	–0.032 $\pm$ 0.069	–0.024 $\pm$ 0.0402	–0.050 $\pm$ 0.047
[Å]	[4]	0.101 $\pm$ 0.105	0.045 $\pm$ 0.072	0.054 $\pm$ 0.033	–0.018 $\pm$ 0.052	–0.28 $\pm$ 0.14
$⟨\delta {\sigma }_{mc}^{H\alpha }⟩$	[1]	–0.092 $\pm$ 0.037	–0.037 $\pm$ 0.038	–0.078 $\pm$ 0.027	0.038 $\pm$ 0.046	0.238 $\pm$ 0.031
$\pm$	[2]	–0.100 $\pm$ 0.035	–0.045 $\pm$ 0.036	–0.084 $\pm$ 0.028	0.038 $\pm$ 0.045	0.240 $\pm$ 0.033
$\Delta {\sigma }_{mc}^{H\alpha }$	[3]	–0.28 $\pm$ 0.12	–0.035 $\pm$ 0.018	–0.092 $\pm$ 0.037	0.134 $\pm$ 0.093	0.233 $\pm$ 0.041
[Å]	[4]	–0.31 $\pm$ 0.13	–0.044 $\pm$ 0.025	–0.11 $\pm$ 0.037	0.136 $\pm$ 0.097	0.240 $\pm$ 0.055
$⟨\delta {\lambda }_{mc}^{\left[\mathrm{NII}\right]}⟩$	[1]	0.102 $\pm$ 0.058	0.077 $\pm$ 0.071	0.083 $\pm$ 0.061	0.058 $\pm$ 0.055	0.031 $\pm$ 0.076
$\pm$	[2]	–0.0290 $\pm$ 0.0402	–0.013 $\pm$ 0.049	–0.039 $\pm$ 0.037	–0.051 $\pm$ 0.037	–0.0740 $\pm$ 0.0702
$\Delta {\lambda }_{mc}^{\left[\mathrm{NII}\right]}$	[3]	0.104 $\pm$ 0.065	0.078 $\pm$ 0.073	0.047 $\pm$ 0.054	0.055 $\pm$ 0.059	0.029 $\pm$ 0.075
[Å]	[4]	–0.068 $\pm$ 0.058	–0.058 $\pm$ 0.076	–0.073 $\pm$ 0.044	–0.16 $\pm$ 0.16	–0.105 $\pm$ 0.089
$⟨\delta {\sigma }_{mc}^{\left[\mathrm{NII}\right]}⟩$	[1]	–0.121 $\pm$ 0.037	–0.066 $\pm$ 0.028	–0.107 $\pm$ 0.034	0.009 $\pm$ 0.047	0.209 $\pm$ 0.045
$\pm$	[2]	–0.129 $\pm$ 0.036	–0.074 $\pm$ 0.031	–0.113 $\pm$ 0.035	0.009 $\pm$ 0.049	0.211 $\pm$ 0.046
$\Delta {\sigma }_{mc}^{\left[\mathrm{NII}\right]}$	[3]	–0.067 $\pm$ 0.043	–0.076 $\pm$ 0.043	–0.063 $\pm$ 0.025	–0.060 $\pm$ 0.045	0.283 $\pm$ 0.033
[Å]	[4]	–0.026 $\pm$ 0.059	–0.038 $\pm$ 0.083	–0.035 $\pm$ 0.041	–0.041 $\pm$ 0.042	0.294 $\pm$ 0.043

Table 5. Table 5: Mean and root mean square values for each stellar subtraction method.

	Unsubtracted	FIT3D	STARLIGHT (MILES)	STARLIGHT (STELIB)	pPXF
$\delta {\lambda }_m^{H\alpha }$	0.060 $\pm$ 0.079	0.022 $\pm$ 0.070	0.020 $\pm$ 0.052	–0.019 $\pm$ 0.046	–0.083 $\pm$ 0.109
$\delta {\sigma }_m^{H\alpha }$	–0.19 $\pm$ 0.11	–0.040 $\pm$ 0.033	–0.091 $\pm$ 0.034	0.087 $\pm$ 0.085	0.238 $\pm$ 0.045
$\delta {\lambda }_m^{\left[\mathrm{NII}\right]}$	0.027 $\pm$ 0.072	0.021 $\pm$ 0.079	0.004 $\pm$ 0.063	–0.023 $\pm$ 0.108	–0.030 $\pm$ 0.089
$\delta {\sigma }_m^{\left[\mathrm{NII}\right]}$	–0.086 $\pm$ 0.054	–0.063 $\pm$ 0.056	–0.080 $\pm$ 0.041	0.021 $\pm$ 0.052	0.249 $\pm$ 0.051

Table 6. Table 6: Mean and root mean square values for each line fitting case.

	case 1	case 2	case 3	case 4
$\delta {\lambda }_c^{H\alpha }$	–0.008 $\pm$ 0.045	0.011 $\pm$ 0.081	–0.016 $\pm$ 0.052	0.013 $\pm$ 0.091
$\delta {\sigma }_c^{H\alpha }$	0.014 $\pm$ 0.065	0.010 $\pm$ 0.066	–0.007 $\pm$ 0.107	–0.017 $\pm$ 0.116
$\delta {\lambda }_c^{\left[\mathrm{NII}\right]}$	0.070 $\pm$ 0.063	–0.041 $\pm$ 0.048	0.062 $\pm$ 0.064	–0.092 $\pm$ 0.092
$\delta {\sigma }_c^{\left[\mathrm{NII}\right]}$	–0.015 $\pm$ 0.066	–0.013 $\pm$ 0.068	0.003 $\pm$ 0.073	0.031 $\pm$ 0.081

Table 7. Table 7: Mean value (and standard deviations) for the H α 𝛼 \alpha and [NII] lines among the different methods and different cases.

	Unsubtracted	FIT3D	STARLIGHT (MILES)	STARLIGHT (STELIB)	pPXF
(1)	(2)	(3)	(4)	(5)	(6)
[1] ${\lambda }_{H\alpha }$	6563.87 $\pm$ 1.02	6563.85 $\pm$ 0.96	6563.85 $\pm$ 1.00	6563.83 $\pm$ 0.98	6563.81 $\pm$ 1.05
[2] ${\lambda }_{H\alpha }$	6563.94 $\pm$ 1.04	6563.89 $\pm$ 0.93	6563.91 $\pm$ 0.98	6563.84 $\pm$ 0.96	6563.74 $\pm$ 1.13
[3] ${\lambda }_{H\alpha }$	6563.88 $\pm$ 1.01	6563.85 $\pm$ 0.95	6563.82 $\pm$ 1.05	6563.83 $\pm$ 0.98	6563.81 $\pm$ 1.05
[4] ${\lambda }_{H\alpha }$	6563.95 $\pm$ 1.08	6563.89 $\pm$ 0.92	6563.91 $\pm$ 0.99	6563.83 $\pm$ 0.95	6563.73 $\pm$ 1.13
[1] ${\sigma }_{H\alpha }$	1.38 $\pm$ 0.35	1.43 $\pm$ 0.36	1.39 $\pm$ 0.35	1.51 $\pm$ 0.45	1.71 $\pm$ 0.40
[2] ${\sigma }_{H\alpha }$	1.37 $\pm$ 0.35	1.43 $\pm$ 0.36	1.39 $\pm$ 0.35	1.51 $\pm$ 0.47	1.71 $\pm$ 0.40
[3] ${\sigma }_{H\alpha }$	1.19 $\pm$ 0.18	1.44 $\pm$ 0.34	1.38 $\pm$ 0.29	1.61 $\pm$ 0.53	1.71 $\pm$ 0.40
[4] ${\sigma }_{H\alpha }$	1.17 $\pm$ 0.16	1.43 $\pm$ 0.34	1.36 $\pm$ 0.30	1.61 $\pm$ 0.54	1.71 $\pm$ 0.42
[1] ${\lambda }_{\left[NII\right]}$	6584.46 $\pm$ 1.01	6584.44 $\pm$ 0.96	6584.44 $\pm$ 1.01	6584.42 $\pm$ 0.99	6584.39 $\pm$ 1.05
[2] ${\lambda }_{\left[NII\right]}$	6584.33 $\pm$ 1.03	6584.34 $\pm$ 1.01	6584.32 $\pm$ 1.03	6584.31 $\pm$ 1.02	6584.28 $\pm$ 0.99
[3] ${\lambda }_{\left[NII\right]}$	6584.46 $\pm$ 1.01	6584.44 $\pm$ 0.96	6584.41 $\pm$ 1.06	6584.41 $\pm$ 0.98	6584.39 $\pm$ 1.05
[4] ${\lambda }_{\left[NII\right]}$	6584.29 $\pm$ 0.99	6584.30 $\pm$ 0.96	6584.28 $\pm$ 1.00	6584.20 $\pm$ 0.96	6584.25 $\pm$ 0.98
[1] ${\lambda }_{\left[NII\right]}$	1.37 $\pm$ 0.35	1.43 $\pm$ 0.36	1.39 $\pm$ 0.35	1.51 $\pm$ 0.47	1.71 $\pm$ 0.39
[2] ${\lambda }_{\left[NII\right]}$	1.37 $\pm$ 0.35	1.43 $\pm$ 0.36	1.39 $\pm$ 0.35	1.51 $\pm$ 0.47	1.71 $\pm$ 0.40
[3] ${\lambda }_{\left[NII\right]}$	1.43 $\pm$ 0.44	1.42 $\pm$ 0.42	1.44 $\pm$ 0.40	1.44 $\pm$ 0.45	1.78 $\pm$ 0.36
[4] ${\lambda }_{\left[NII\right]}$	1.48 $\pm$ 0.44	1.46 $\pm$ 0.43	1.47 $\pm$ 0.40	1.46 $\pm$ 0.44	1.79 $\pm$ 0.36

Table 8. Table 8: Mean values (and standard deviations) for the H α 𝛼 \alpha and [NII] lines when considering the four cases of line fitting for each method.

	Unsubtracted	FIT3D	STARLIGHT (MILES)	STARLIGHT (STELIB)	pPXF
(1)	(2)	(3)	(4)	(5)	(6)
${\lambda }_{H\alpha }$	6563.91 $\pm$ 0.50	6563.87 $\pm$ 0.46	6563.87 $\pm$ 0.49	6563.83 $\pm$ 0.47	6563.77 $\pm$ 0.53
${\sigma }_{H\alpha }$	1.28 $\pm$ 0.14	1.43 $\pm$ 0.17	1.38 $\pm$ 0.16	1.56 $\pm$ 0.23	1.71 $\pm$ 0.20
${\lambda }_{\left[NII\right]}$	6584.39 $\pm$ 0.49	6584.38 $\pm$ 0.47	6584.36 $\pm$ 0.50	6584.33 $\pm$ 0.48	6584.33 $\pm$ 0.50
${\sigma }_{\left[NII\right]}$	1.43 $\pm$ 0.19	1.46 $\pm$ 0.19	1.44 $\pm$ 0.18	1.51 $\pm$ 0.22	1.76 $\pm$ 0.19

Table 9. Table 10: Reduced chi–square χ r ​ e ​ d 2 subscript superscript 𝜒 2 𝑟 𝑒 𝑑 \chi^{2}_{red} values for the different spaxels and models.

Spectrum	${\chi }_{red}^2$	${\chi }_{red}^2$	${\chi }_{red}^2$	${\chi }_{red}^2$	${\chi }_{red}^2$	Best fit	Secondary
name	(Unsubtracted)	(FIT3d)	(STARLIGHT/MILES)	(STARLIGHT/STELIB)	(pPXF)	model	best model
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
a	2.55	1.14	1.73	1.58	1.91	FIT3d	STELIB
b	2.64	1.43	1.25	1.86	1.74	STARLIGHT/MILES	FIT3d
c	1.53	0.73	0.80	0.84	1.04	FIT3d	MILES
d	3.05	1.36	2.30	1.65	2.40	FIT3d	STELIB
e	2.09	1.55	1.51	2.11	2.28	STARLIGHT/MILES	FIT3d
f	1.03	0.87	0.98	1.07	0.92	FIT3d	pPXF
g	1.09	1.14	1.06	1.02	1.21	STARLIGHT/STELIB	MILES
h	0.79	0.85	0.63	0.83	0.81	STARLIGHT/MILES	Unsubtracted
i	1.00	0.93	0.98	1.08	1.09	FIT3d	MILES
j	0.90	0.84	0.86	0.91	0.74	pPXF	FIT3d
k	0.90	1.00	0.85	0.84	0.67	pPXF	STELIB
l	0.77	0.81	0.99	0.86	0.50	pPXF	Unsubtracted

Equations

$${\overline{\lambda }}_i=\frac{1}{N_m}\frac{1}{N_c}\sum _{m=1}^{N_m}\sum _{c=1}^{N_c}{\lambda }_{i,mc}\text{leavevmode};\text{leavevmode}{\overline{\sigma }}_i=\frac{1}{N_m}\frac{1}{N_c}\sum _{m=1}^{N_m}\sum _{c=1}^{N_c}{\sigma }_{i,mc}$$

$$\delta {\lambda }_{i,mc}={\lambda }_{i,mc}-{\overline{\lambda }}_i\text{leavevmode};\text{leavevmode}\text{leavevmode}\delta {\sigma }_{i,mc}={\sigma }_{i,mc}-{\overline{\sigma }}_i.$$

$$\Delta {\overline{\lambda }}_i=\sqrt{\frac{{\sum }_{mc}\delta {\lambda }_{i,mc}^2}{N_m{N}_c}}\text{leavevmode};\text{leavevmode}\text{leavevmode}\Delta {\overline{\sigma }}_i=\sqrt{\frac{{\sum }_{mc}\delta {\sigma }_{i,mc}^2}{N_m{N}_c}}$$

$$⟨\delta {\lambda }_{mc}⟩=\frac{{\sum }_i\delta {\lambda }_{i,mc}}{N_i}$$

$$⟨\delta {\sigma }_{mc}⟩=\frac{{\sum }_i\delta {\sigma }_{i,mc}}{N_i}$$

$$⟨\delta {\lambda }_c⟩=\frac{{\sum }_m\delta {\lambda }_{mc}}{N_m}$$

$$⟨\delta {\sigma }_c⟩=\frac{{\sum }_m\delta {\sigma }_{mc}}{N_m}$$

$$⟨\delta {\lambda }_m⟩=\frac{{\sum }_c\delta {\lambda }_{mc}}{N_c}$$

$$⟨\delta {\sigma }_m⟩=\frac{{\sum }_c\delta {\sigma }_{mc}}{N_c}$$

$$\Sigma {\lambda }_{i,mc}=\frac{\delta {\lambda }_{i,mc}}{\Delta {\lambda }_{i,mc}}\text{leavevmode};\text{leavevmode}\text{leavevmode}\Sigma {\sigma }_{i,mc}=\frac{\delta {\sigma }_{i,mc}}{\Delta {\sigma }_{i,mc}}.$$

$$\Sigma ({\lambda }_{i,mc}{)}^{Corr}=\frac{\delta {\lambda }_{i,mc}}{\Delta {\lambda }_{i,mc}^{Corr}};\Sigma ({\sigma }_{i,mc}{)}^{Corr}=\frac{\delta {\sigma }_{i,mc}}{\Delta {\sigma }_{i,mc}^{Corr}}.$$

$$\Delta {\lambda }_{i,mc}^{Corr}=\sqrt{\Delta {\lambda }_{i,mc}^2+\Delta {\overline{\lambda }}_i^2};\Delta {\sigma }_{i,mc}^{Corr}=\sqrt{\Delta {\sigma }_{i,mc}^2+\Delta {\overline{\sigma }}_i^2}.$$

$$\Delta \overline{\lambda }=\Delta \overline{\sigma }=CF\times \frac{{\mathtt{FWHM}}_{instr}^{\mathtt{MUSE}}}{2.354}=CF\times {\Delta }_{\mathtt{MUSE}}$$

$$CF=\frac{\Delta \overline{\lambda }}{{\Delta }_{\mathtt{MUSE}}}\lesssim \frac{0.2}{1.1}=0.18\sim \mathrm{0.2.}$$

$${\chi }_{red}^2=\sum _{i=1}^{N_{\lambda }}{\left(\frac{O_i-M_i}{{\sigma }_i}\right)}^2/N_{\lambda }$$

$$M_i=M_{i,\star }+\sum _{k=1}^{N_c}{M}_k^{i,gas}/N_c=M_{i,\star }+M_{i,gas}.$$

Full text

*Aims. * We study how the use of several stellar subtraction methods and line fitting approaches can affect the derivation of the main kinematic parameters (velocity and velocity dispersion fields) of the ionized gas component.

*Methods. * The target of this work is the nearby galaxy NGC 2906, observed with the MUSE instrument at Very Large Telescope (VLT). A sample of twelve spectra is selected from the inner (nucleus) and outer (spiral arms) regions, characterized by different ionization mechanisms. We compare three different methods to subtract the stellar continuum (FIT3D, STARLIGHT and pPXF), combined with one of the following stellar libraries: MILES, STELIB and GRANADA+MILES.

Results. * The choice of the stellar subtraction method is the most important ingredient affecting the derivation of the gas kinematics, followed by the choice of the stellar library and by the line fitting approach. In our data, typical uncertainties in the observed wavelength and width of the H$\alpha$ and [NII] lines are of the order of $\langle\delta\lambda\rangle_{\rm rms}\sim 0.1$ Å and $\langle\delta\sigma\rangle_{\rm rms}\sim 0.2$ Å (i.e. $\sim 5$ and $10$ km s-1*, respectively). The results obtained from the [NII] line seem to be slightly more robust, as it is less affected by stellar absorption than H$\alpha$. All methods considered yield statistically consistent measurements once a mean systemic contribution $\Delta\bar{\lambda}=\Delta\bar{\sigma}=0.2\,\Delta_{\tt MUSE}$ is added in quadrature to the line fitting errors, where $\Delta_{\tt MUSE}=1.1$ Å $\sim 50$ km s*-1* denotes the instrumental resolution of the MUSE spectra.

*Conclusions. * Although the subtraction of the stellar continuum is critical in order to recover line fluxes, any method (including none) can be used in order to measure the gas kinematics, as long as an additional component $\Delta\bar{\lambda}=\Delta\bar{\sigma}=0.2\,\Delta_{\tt MUSE}$ is added to the error budget.

1 Introduction

The kinematic characterization of different galaxy populations is a key observational input to distinguish between different evolutionary paradigms, such as the ‘major merger’ (Tacconi et al., 2006, 2008; Kartaltepe et al., 2012) and ‘steady cold gas accretion’ (Kereš et al., 2005; Ocvirk et al., 2008; Dekel et al., 2009b; Ceverino et al., 2010) scenarios, since it allows one to determine the fraction of rotating disks to mergers at different cosmic epochs (see e.g., Genzel et al., 2001; Tacconi et al., 2008; Dekel et al., 2009a; Förster Schreiber et al., 2009; Lemoine-Busserolle & Lamareille, 2010; Lemoine-Busserolle et al., 2010; Bournaud et al., 2011; Epinat, 2011; Glazebrook, 2013; Bellocchi et al., 2016). Kinematics is also important in order to study the physical processes that govern the formation and evolution of the galaxies, providing a powerful diagnostic to infer the main source of dynamic support (Puech et al., 2007; Epinat et al., 2009; Cappellari et al., 2013; Zhu et al., 2018), to distinguish between relaxed virialized systems and merger events (Flores et al., 2006; Shapiro et al., 2008; Bellocchi et al., 2012; Barrera-Ballesteros et al., 2015; Bellocchi et al., 2016), to infer fundamental galaxy quantities like the dynamical mass (Colina et al., 2005; Bellocchi et al., 2013; Aquino-Ortíz et al., 2018) and also to detect and characterize radial motions associated with feedback mechanisms, like massive gas outflows (Shapiro et al., 2009; Rupke & Veilleux, 2013; López-Cobá et al., 2017a, b; Maiolino et al., 2017).

With the aid of integral field spectroscopy (IFS), we can resolve the kinematical status and internal processes at work within galaxies and more clearly understand the role of the dominant mechanisms involved at early epochs of galaxy formation. This technique provides much more information of the observed galaxies than was previously available through classical long–slit spectroscopy: on the one hand, the full two–dimensional spatial information (i.e., maps, rather than radial profiles) can be reconstructed at high resolution in regions of bright emission; on the other hand, one of the advantages of IFS is that large areas may be combined in order to properly characterize weak signals, such as in the galaxy outskirts. The observed spectrum at each spaxel (or area) includes gas in different phases, traced by distinctive features at different wavelengths (e.g., H$\alpha$ emission allows one to study the warm ionized emission, NaD absorption traces the neutral gas, the CO emission helps to trace the cold molecular clouds, etc.; see Bekeraité et al., 2016; Leung et al., 2018; Levy et al., 2018) as well as the underlying stellar population (a smooth continuum with absorption features. To properly characterize their physical properties, we need to apply methods that are able to separate the contribution of gas and stars to the total emission. To this aim, in the last ten years, many routines have been developed to subtract the stellar continuum. Some of them like STARLIGHT (Cid Fernandes et al., 2005), FIT3D (Sánchez et al., 2016c), pPXF (Cappellari & Emsellem, 2004), STECKMAP (Ocvirk et al., 2006), sedfit (Walcher et al., 2006) and PyPARADISE (Husemann et al. in prep.) are widely used nowadays. These routines allow to study in detail the underlying stellar population(s) of the source and derive interesting parameter/properties such as the total stellar mass, the star formation rate (e.g., Pérez et al., 2013; González Delgado et al., 2017), the age and metallicities, as well as the kinematics of the stars.

In this paper we want to investigate how the choice of different stellar subtraction methods affects the derivation of the gas kinematics, traced by the ionized component through the H$\alpha$ and [NII] lines. We also analyze how the choice of different line fitting approaches influence the derivation of the kinematic parameters (mean velocity and velocity dispersion). The target of this study is a nearby galaxy observed with MUSE/VLT, from which we obtain high resolution spectra, very suitable to our aims (a related discussion has been recently carried out by Belfiore et al. (2019) in the context of the MaNGA Data Analysis Pipeline). We describe the observations in Section 2. In Section 3 we present the selection of the spectra, a brief description of the methods involved for the stellar continuum subtraction, and the line fitting approaches considered. Section 4 is devoted to the discussion of the results obtained by different methods, and our main conclusions are summarized in Section 5. Appendix A presents some details on the line fitting analysis along with complementary kinematic results and the quality of the fits is discussed in Appendix B. Throughout the paper we consider H0 = 70 km s*-1* Mpc*-1*, $\Omega_{M}$ = 0.3 and $\Omega_{\Lambda}$ = 0.7.

2 Observations

NGC 2906111This galaxy was previously observed with the integral field spectrograph PMAS/PPak mounted on the 3.5m telescope at the Calar Alto Observatory. It is part of the Third Public Data Release (DR3) of the Calar Alto Legacy Integral Field Area (CALIFA) survey (Sánchez et al., 2016a). is a nearby (z = 0.007138, d = 30 Mpc) ‘normal’ star–forming galaxy which was observed during the science verification of the MUSE instrument (Bacon et al., 2004), and included in the pilot study of the All–weather MUse Supernova Integral–field of Nearby Galaxies (AMUSING; Galbany et al., 2016) given that it hosts the type IIn supernova (SN) 2005ip (Boles et al., 2005; Modjaz et al., 2005). According to its morphology, this object is classified as Scd galaxy. In Fig. 1 the continuum and the H$\alpha$222The H$\alpha$ flux emission is actually the flux peak of the line having removed the continuum emission computed in a region close to the H$\alpha$–[NII] complex. emissions of this object are shown.

The MUSE instrument provides a FoV of 1 arcmin2 with spaxels of 0.2${}^{\prime\prime}\times 0.2^{\prime\prime}$ in Wide Field Mode (WFM). At this redshift, the arcsec to kpc conversion gives 0.147 kpc arcsec*-1*. In this configuration the data cube is composed by 320 $\times$ 322 spaxels with 3680 elements of spectral sampling. The wavelength coverage ranges between 4750–9350 Å  with a mean resolution of R$\sim$3000333See the MUSE spectral resolution (R) as a function of the wavelength at:

[https://www.eso.org/sci/facilities/paranal/instruments/

muse/inst.html](https://www.eso.org/sci/facilities/paranal/instruments/muse/inst.html). At the H$\alpha$ wavelength, it corresponds to a spectral resolution R $\sim$ 2500, deriving a FWHM of 2.6 Å. and sampling of 1.25 Å. This object has been observed with a seeing (FWHM) of 0.88*′′* ($\sim$0.13 kpc). With the MUSE instrument it is possible to combine good seeing conditions with large FoV and a spatial sampling sufficient to properly sample the point spread function (PSF). Thus, it is the suitable instrument to study in detail the kinematics of our object.

3 Analysis

In this section we first present how the spectra have been selected, then describe the codes used to subtract the stellar continuum emission from the observed unsubtracted data as well as the stellar libraries involved, and finally discuss the different approaches that have been applied to fit the H$\alpha$–[NII] line complex.

3.1 Selected spectra

This work focuses on the analysis of twelve spectra. Half of them are located in the nuclear region, while the remaining ones are placed in different parts of the galaxy. For practical reasons, we refer to the different spaxels using letters (a, b, ..., k, l), as shown in Figs. 1,  2 and  4, ordered according to their (decreasing) log([NII]$\lambda$6583/H$\alpha$) flux ratios444The flux ratios have been computed applying the most general line fitting approach (case [4]; see Sect. 3.3)..

The spatial distribution of the selected spaxels is shown in Fig. 1. The intensity peak in the continuum image (spectrum d) has been considered the origin of our system of reference, and the offsets of all spectra with respect to that origin are listed in Tab 3. The minimum separation between the spectra is 0.6*′′, reached in the region close to the nucleus, up to a maximum separation of 11.8′′* from the centre.

We have tried to ensure that these spectra are representative of a wide range of physical conditions and signal–to–noise (S/N) ratios. The original spectra, extracted from from the MUSE datacube, are shown on the left column of Fig. 2, whereas the other columns display the results obtained after subtracting the stellar continuum with each of the four methods described below. Although the precise line fluxes and equivalent widths are rather sensitive to the adopted method, it is evident from Fig. 2 that the values of the $\log([{\rm NII}]\lambda 6583/{\rm H}\alpha)$ ratio illustrate the whole range of values typically found in a Baldwin et al. (1981) diagram (see also Tab. 3).

3.2 Stellar continuum subtraction methods

In order to characterize nebular emission, the stellar continuum emission needs to be subtracted. The H$\alpha$ and H$\beta$ emission line profiles may be significantly affected from Balmer absorption especially in evolved stellar populations (e.g., Sarzi et al., 2007; Cresci et al., 2015; Gomes et al., 2016).

This is clearly visible in NGC 2906, where the H$\alpha$ emission in the central spectra is quite absorbed (from a to d spectra in Fig. 2), while for the outer spectra this effect does not seem to be so relevant (from f to l spectra in Figs. 2 and 3).

Here we remove the stellar contribution using three different software programs that are publicly available and have been widely used in the literature:

FIT3D555Specifically, the pipeline Pipe3D was used which is based on the FIT3D code.

http://www.astroscu.unam.mx/~sfsanchez/FIT3D/ (Sánchez et al., 2016b, c). 2. 2.

STARLIGHT666http://www.starlight.ufsc.br (Cid Fernandes et al., 2004, 2005, 2013). 3. 3.

Penalized Pixel–Fitting777http://www-astro.physics.ox.ac.uk/~mxc/software/ (pPXF; Cappellari & Emsellem, 2004; Cappellari, 2017).

These spectral synthesis codes combine the spectra from a base of simple stellar populations of various ages and metallicities in order to match an observed spectrum. The search for an optimal model and the coefficients of the population vector associated with the base elements also allows to derive interesting outputs, such as the star formation and chemical enrichment histories of a galaxy or region, as well as to recover non–linear parameters, such as the amount of dust extinction and the main kinematic parameters (velocity shift and dispersion) of both stars and gas. The reader is referred to the corresponding articles for an in–depth discussion of each method.

In any case, the characteristics of the stellar libraries adopted as a basis are as important as the fitting algorithm. The choice of the stellar library generally depends on the resolution of the data. In our case, MUSE has a nominal spectral resolution FWHM $\sim$ 2.6888The intrinsic spectral resolution is slightly lower than the theoretical one (FWHM $\sim$ 2.2 Å). Å, and we combined each code with one of the following stellar libraries in order to achieve a similar (or lower) resolution:

STELIB: empirical library from Bruzual & Charlot (2003, ; hereafter BC03); 2. 2.

MILES: ‘Medium-resolution Isaac Newton Telescope Library of Empirical Spectra’ from Vazdekis et al. (2010); 3. 3.

GRANADA & MILES: a combination of the theoretical stellar libraries from González Delgado et al. (2005) and González Delgado & Cid Fernandes (2010) with Vazdekis et al. (2010).

The STELIB999The STELIB stellar library can be found at http://webast.ast.obs-mip.fr/stelib. library (Le Borgne et al., 2003) consists of an homogeneous library of 249 stellar spectra in the visible range (3200 to 9500 Å), with an intermediate spectral resolution ($\leq$3 Å) and sampling (1 Å). Only 187 over 249 stars have measured metallicity and can be used to compute the predicted spectra. This library includes stars of various spectral types and luminosity classes, spanning a relatively wide range in metallicity (Z = 0.0001 – 0.05).

The MILES101010The MILES stellar library can be found at [http://www.iac.es/proyecto/miles/pages/webtools/tune-

ssp-models.php](http://www.iac.es/proyecto/miles/pages/webtools/tune-ssp-models.php). stellar library (Vazdekis et al., 2010) consists of stars spanning a large range in atmospheric parameters. The 985 spectra were obtained at the 2.5m Isaac Newton Telescope (INT) and cover the range 3525–7500 Å (Sánchez-Blázquez et al., 2006) at 2.51 Å (FWHM) spectral resolution (Falcón-Barroso et al., 2011). These models cover ages between $\sim$0.1 and $\sim$18 Gyr and metallicities between 0.0004 – 0.03. Vazdekis/Miles models are based on the previous models from Vazdekis (1999) and Vazdekis et al. (2003) and they use Padova2000 isochrones. The Padova isochrones (1994 and 2000) are presented in Girardi et al. (2002)111111Available at http://pleiadi.pd.astro.it/ based on the (solar scaled mixture) tracks from Bertelli et al. (1994) and Girardi et al. (2000), respectively.

GRANADA121212The GRANADA library is a high resolution library. & MILES stellar library (gsd156) comprises 156 templates that cover 39 stellar ages (1 Myr to 13 Gyr), and 4 metallicities (Z = 0.002, 0.008, 0.02 and 0.03). These templates have been extracted from a combination of the synthetic stellar spectra from the GRANADA ((Martins et al., 2005)) and the SSP libraries provided by the MILES project. This SSP–library uses the Salpeter (1955) Initial Mass Function (IMF) and Girardi et al. (2000) and Geneva131313The Geneva isochrones have been computed with the isochrone program presented in Meynet (1995) and they follow the prescriptions quoted in Cerviño et al. (2001) from the evolutionary tracks from Schaller et al. (1992), Charbonnel et al. (1993), Schaerer et al. (1993b) and Schaerer et al. (1993a). tracks. This library is described in detail in Cid Fernandes et al. (2013).

The ‘code–stellar library’ combinations we used are:

the FIT3D code with the GRANADA and MILES libraries; 2. 2.

the STARLIGHT code with the MILES library (see Galbany et al., 2016) and references therein for more information141414This model is considered a ‘modified’ version of the BC03 models (i.e., CB07; see Bruzual (2007) for further details).); 3. 3.

the STARLIGHT code with the STELIB library; 4. 4.

the pPXF code with the MILES library.

In Tab. 1 we summarize all combinations ‘code–stellar library’ used in this analysis, including details on the isochrones and the age of the stellar populations considered.

3.3 Line fitting cases

In this work we focus on the H\alpha$$\lambda6563Å and the forbidden [NII]\lambda$$\lambda6548, 6583Å emission lines151515We actually do not discuss the results of [NII]$\lambda$6548Å since this line is three times fainter than [NII]$\lambda$6583Å (according to atomic physics) and more subjected to higher uncertainties in the fit. Furthermore, in most cases, it would share the same kinematics as the [NII]$\lambda$6583Å line.. The routine MPFITEXPR (implemented in IDL code by Markwardt, 2009) is quite commonly used to derive the kinematic parameters (e.g., Arribas et al., 2008; Bellocchi et al., 2012, 2013; Cazzoli et al., 2014; Cairós et al., 2015; Bellocchi et al., 2016; Poggianti et al., 2017). This algorithm allows to fit the observed lines of the individual spectra to Gaussian profiles and the continuum emission as a line with a certain slope. It derives the best set of lines that match the available data. In case of adjusting multiple lines, the line flux ratios can be fixed according to the atomic physics (as in the case of the [NII] doublet). This routine gives in output the wavelength (i.e., centroid $\lambda$), the width ($\sigma$) and flux of the line. It also allows one to combine in different ways the $\lambda$ and $\sigma$ parameters keeping them fixed or let free to vary. Thus, when fitting the H$\alpha$–[NII] complex, there are four cases that can be taken into account:

$\lambda$ & $\sigma$ FIXED: the wavelengths of the different lines are fixed according to the atomic physics and their widths are constrained to be equal for all lines (e.g., Bellocchi et al., 2013, 2016); 2. 2.

$\lambda$ FREE, $\sigma$ FIXED: the wavelengths of the different lines are free to vary but their widths are constrained to be equal for all lines; 3. 3.

$\lambda$ FIXED, $\sigma$ FREE: the wavelengths of the different lines are fixed according to the atomic physics but their widths are left free to vary; 4. 4.

$\lambda$ & $\sigma$ FREE: both wavelengths and widths are free to vary for all lines.

All these cases are summarized in Tab. 2. In the next section the results of our analysis are shown and discussed. It is worth mentioning that in this work the MUSE instrumental profile ($\sigma_{INS}$) has not been subtracted161616To have a rough estimate of the instrumental profile we apply the line fitting case 1 to the [OI]$\lambda$6300 sky line in a relatively large region of the galaxy, where the [OI] emission is present. We derive a typical (median) value for the wavelength of $\lambda_{INS}\leavevmode\nobreak\ \sim$ 6299.99 Å and $\sigma_{INS}\sim$ 0.88 Å, with the respective mean values of 6299.99 ($\pm$ 0.91) Å and 0.96 ($\pm$ 0.48) Å..

4 Results

In this section we present and discuss the kinematic results obtained when different stellar subtraction methods and line fitting approaches are considered. The raw observed data (hereafter, unsubtracted) are also included in the analysis in order to study the case in which the stellar emission is not removed at all and thus understand the importance of removing (or not) the stellar component from the data for the derivation of the kinematics of the gas component.

In Sect. 4.1 we focus on the results derived for each individual spectrum to better study the dispersion related to the choice of a specific stellar subtraction method as well as the one associated to the use of a particular line fitting approach. In Sect. 4.2 we then analyze the mean trends obtained when considering all spectra averaged according to the number of methods or the number of cases considered.

4.1 Mean and dispersion for individual spectra

In this section the derivation of the kinematic results of each individual spectrum ($\lambda_{i,mc}$, $\sigma_{i,mc}$) is presented, where the index $i$ runs over our 12 individual spectra, and a combination of twenty values per spectrum is analyzed (i.e., 5 methods $\times$ 4 cases; $N_{mc}=N_{m}\,N_{c}=20$). We remind that, according to our selection criterion, the [NII] contribution is very significant in the first six spectra (from spectrum [a] to [f]), while the H$\alpha$ flux emission dominates in the remaining six spectra (from spectrum [g] to [l]).

We first determined the wavelength $\lambda_{i,mc}$ and the velocity dispersion $\sigma_{i,mc}$ for both lines in each spectrum (Fig. 4) as derived when applying different line fitting approaches (‘cases’) and continuum–subtraction methods. The more asymmetric shape of the H$\alpha$ profile in the first spectra, characterized by lower S/N as a consequence of the H$\alpha$ absorption, reflects in larger dispersions when deriving the kinematic values when using different values. Conversely, the quality of the fit of the [NII] line is also affected by the corresponding S/N ratio, although the differences are not as marked as for the Balmer line.

Considering that each measurement provides a wavelength $\lambda_{i,mc}$ and a velocity dispersion $\sigma_{i,mc}$, we defined the mean values for each spectrum among all twenty combinations as

[TABLE]

and we characterize the behavior of the different methods and cases in terms of the offset from the appropriate mean value:

[TABLE]

As the true kinematic parameters of every individual spectrum are unknown, we cannot determine the optimal combination of subtraction method and fitting approach (in each particular case), but we may quantify the uncertainties associated to such choice in a statistical sense, computing the standard deviations

[TABLE]

arising from our $N_{mc}=20$ combinations.

The mean values of the wavelength $\overline{\lambda}_{i}$ and the velocity dispersion $\overline{\sigma}_{i}$ obtained for each spectrum, along with the measurement dispersions (uncertainty ranges) $\Delta\overline{\lambda}_{i}$ and $\Delta\overline{\sigma}_{i}$, are listed in Tab. 3 and shown in Fig. 4. The latter, as well as the fractional uncertainties ($\Delta\overline{\lambda}_{i}$/$\overline{\lambda}_{i}$ and $\Delta\overline{\sigma}_{i}$/$\overline{\sigma}_{i}$) are plotted in Fig. 5. We find that, for H$\alpha$, $\Delta\overline{\lambda}_{i}$ decreases from 0.3 Å  down to 0.03 Å when moving from spectrum [a] to [l], and a (much milder) decreasing trend is also found for the velocity dispersion uncertainty $\Delta\overline{\sigma}_{i}$, which varies from $\sim$0.4 down to $\sim$0.1. The respective fractional uncertainties for the wavelength parameter, $\Delta\overline{\lambda}_{i}$/$\overline{\lambda}_{i}$, thus cover the range of values between 4.9$\times$10*-5* Å down to 3.9$\times$10*-6* Å from spectrum [a] to [l], while those obtained for the velocity dispersion, $\Delta\overline{\sigma}_{i}$/$\overline{\sigma}_{i}$, varies from $\sim$ 0.21 down to $\sim$0.08. The opposite trend is derived for the [NII] line, reflecting the line intensity ratios, although more stable and constant values are obtained for both line wavelength and width for different spaxels.

It is important to bear in mind that, on the one hand, we are quantifying the dispersion of the measured parameters with respect to the average value, which may be different from the true solution. On the other hand, assuming that the estimated central wavelengths and line widths inferred from the various codes and combinations of stellar population synthesis (SPS) models are statistically independent is conservative, in the sense that any off–diagonal terms in the covariance matrix would likely reduce the quoted systematic uncertainties. For these reasons, our estimate of the dispersion between codes may be actually regarded as a upper limit to the true systematic uncertainty associated to the choice of a specific method.

4.2 Trends for different methods and cases

Figure 6 shows the mean and standard deviation

[TABLE]

averaged over the $N_{i}=12$ spaxels, for each emission line, stellar continuum subtraction method, and line fitting case. Numerical values are reported on Tab. 4.

In Tabs. 5 and 6 the results from further averaging over the $N_{m}=5$ methods

[TABLE]

for each line fitting case, as well as over the $N_{c}=4$ cases

[TABLE]

at fixed stellar subtraction method are considered, respectively. These results are shown in Fig. 7.

Our results show clear systematic trends depending on the stellar continuum subtraction algorithm. The measurements obtained from the FIT3d method and those derived after applying the STARLIGHT (MILES) method are relatively similar, and both prescriptions tend to yield slightly larger wavelengths and lower widths (with respect to the mean values) for either line, regardless of the line fitting case. On the other hand, higher line widths are found when the STARLIGHT (STELIB) and pPXF algorithms are used, with a mild bias towards lower wavelengths.

In particular, for the H$\alpha$ line the wavelength deviations follow a decreasing trend, from higher to lower wavelengths (from 0.06 Å to –0.083 Å) when moving from the Unsubtracted to the pPXF methods. On the other hand, the width deviations increase from a mean value of –0.19 Å for the Unsubtracted method up to $\sim$ 0.24 Å for the pPXF one. When taking into account the [NII] line, all but pPXF method derived similar wavelength and width deviations with respect to the total mean values. However, they still follow a similar trend (i.e., decreasing or increasing if we consider the wavelength or width deviations, respectively) than that shown for the H$\alpha$ line but covering a smaller range of values (see Tab. 5). It is important to note, once again, that this does not necessarily imply that they are closer to the unknown true solution.

Focusing on the results derived from the STARLIGHT code with different stellar libraries (i.e., MILES, STELIB) and comparing with the measurements obtained with pPXF based on the MILES library, it is evident that the details of the minimization algorithm are at least as important as the choice of spectral basis. Starting from the same data, each program converges to a different solution. Subtle differences in the residual spectra are apparent in Fig. 2, and it is not surprising that these differences propagate into the final kinematic parameters. In our case, the difference when using the same method with different stellar library (STARLIGHT using MILES or STELIB) or the same stellar library in different stellar subtraction method (STARLIGHT and pPXF using MILES) is at least a factor of $\times$2, deriving a larger gap among the kinematic values ($\lambda$, $\sigma$) when different methods are considered.

Regarding the line fitting approach, we found a certain similarity between the cases that share the same kinematics for the velocity dispersion parameter: i.e., when $\sigma$ is fixed (cases [1] and [2]) and when $\sigma$ is free to vary (cases [3] and [4]). We derived the same trends for the wavelength and the velocity dispersion of both H$\alpha$ and [NII] emission lines. The magnitude of the uncertainties is similar to those stemming from the stellar subtraction method, but the effect is much more random rather than systematic.

In general, the two lines yield consistent measurements, although the differences in the recovered kinematics arising from using different methods are somewhat stronger for the Balmer line. For [NII], all methods result in fairly similar values, with the only possible exception of the velocity dispersion obtained when pPXF is used to subtract the stellar continuum. Velocity dispersion is actually more sensitive to the adopted prescription than the mean velocity derived from the line centroid, causing the larger spread already observed in Figs 4 and 5.

The kinematic results obtained from fitting the H$\alpha$–[NII] complex (see Fig. 7, bottom) show a more robust determination of the global kinematic properties with respect to the adopted line fitting approach (right panel), as the strongest (and thus more stable) line tends to dominate the result, and random deviations in each individual line tend to compensate each other on average (see also Belfiore et al., 2019). However, the systematic differences arising from the continuum subtraction method (left panel) persist, as they are present in both emission lines.

5 Discussion

In order to gauge whether these uncertainties are statistically compatible with the errors returned by the MPFITEXPR routine, we have normalized (i.e., standardized) the deviation of each individual measurement from the mean ($\delta\lambda_{i,mc}$, $\delta\sigma_{i,mc}$) using the line fitting errors derived from the output covariance matrix ($\Delta\lambda_{i,mc}$, $\Delta\sigma_{i,mc}$):

[TABLE]

The normalized deviations give us an estimate of the ‘number of $\sigma$’ that each result departs from the mean value, assuming that the line fitting errors should be representative of the expected standard deviation. As can be readily seen in Fig.8, the dispersion of these variables is often well above unity, as expected from Gaussian statistics, for both H$\alpha$ and [NII] lines.

In the next step we modified the Eq. 10 including along with the line fitting errors an extra correction term accounting for the contribution coming from the choice of a specific stellar subtraction methods ($\Delta\overline{\lambda}_{i}$ and $\Delta\overline{\sigma}_{i}$) in the computation of the standard deviation. The new (i.e., corrected) formulas are then modified:

[TABLE]

where

[TABLE]

For both lines we derive a mean dispersion of $\Delta\bar{\lambda}\sim$ 0.13 Å for the wavelength while $\Delta\bar{\sigma}\sim$ 0.2 Å for the velocity dispersion. In order to consider a general value for both parameters and for both lines we take into account the highest dispersion value derived for the wavelength and velocity dispersion parameters ($\sim$0.2 Å), as derived for the velocity dispersion of the H$\alpha$ line.

As can be seen in in Fig. 9, applying the same mean value for $\Delta\bar{\lambda}$ and $\Delta\bar{\sigma}$ we finally derived smaller dispersions (of order unity) for both axes and for both lines. Within the (corrected) 1-$\sigma$ error bar, the results of all methods can be considered statistically consistent with each other.

We then propose the following prescription to compute the correction factor (CF) that needs to be considered to take into account the uncertainty associated to the choice of a specific stellar subtraction method. The formula can be related to the instrumental profile of the instrument, $\Delta_{\tt MUSE}$:

[TABLE]

In the case of the MUSE data the (theoretical) instrumental resolution is of the order of 1.1 Å (FWHM=2.6 Å). Considering the mean uncertainty value of 0.2 Å, the correction factor (CF) can be derived as:

[TABLE]

In this particular case, the (observed, computed) instrumental resolution $\Delta_{\tt MUSE}$ is of the order of 0.9 Å, deriving a CF of the order of $\sim$0.22. To check the validity of this equation the same analysis should be applied to other local sources observed with MUSE, as well as other instruments, in order to test whether our simple prescription can be extended to other instrumental resolutions.

6 Summary and Conclusions

We have compared the kinematic results obtained when applying three different stellar subtraction methods to a sample of twelve spectra in the nearby normal star–forming galaxy NGC2906, observed with MUSE instrument on the VLT. The spectra have been selected according to their different ionization phases, located in different parts of the galaxy, from the nuclear regions to the outer spiral arms.

In order to study and compare the kinematics of the gas component, as that derived from applying the stellar subtraction methods, with that derived from the unsubtracted (gas and stars) one, we have considered the following stellar subtraction methods: FIT3D, STARLIGHT and pPXF. Each of them has been combined with a specific stellar library: MILES for the STARLIGHT and pPXF methods, STELIB for the STARLIGHT method and a combination of two stellar libraries, GRANADA & MILES for the FIT3D one. From their combination we obtained four different data sets of the gaseous ionized component. To better characterize the dispersion when deriving the kinematic parameters we also include different line fitting approaches. In particular, we have considered all possible ways to combine the wavelength and the velocity dispersion parameters allowing them to be free to vary or letting them fixed according to the atomic physics.

The main results of the present study can be summarized as follows:

• $\bullet$

When comparing the kinematic deviations (for both H$\alpha$ and [NII] lines) derived for the raw unsubtracted data (gas and star; no stellar subtraction is applied) with those derived when applying different stellar subtraction methods for the whole sample of spectra we see that the kinematic results obtained with FIT3D are very similar to that derived in the STARLIGHT (MILES) case whereas high kinematic deviations are found when using the STARLIGHT (STELIB) and pPXF;

• $\bullet$

When averaging over all cases, the velocity dispersion $\sigma$ is the parameter most affected by the choice of the different combinations method--line fitting, whose deviations range between ${}^{+0.3}_{-0.3}$ Å, while the wavelength $\lambda$ parameter is less affected (<$${}^{+0.1}_{-0.2} Å or ${}^{+5}_{-9}$ km s*-1*).

• $\bullet$

When averaging over all methods, the kinematic deviations obtained for the H$\alpha$ and [NII] lines ($\langle\delta\lambda_{m}\rangle$, $\langle\delta\sigma_{m}\rangle$) show *point–*symmetric distribution for the cases [1]--[4] and the cases [2]--[3].

• $\bullet$

When analyzing the kinematic results for each individual spectrum we derive larger dispersion (error bars) in the kinematic parameters when the line taken into account show a low S/N ratio: this happened for both H$\alpha$ and [NII] lines, deriving a typical value of $\lesssim$ 0.2 Å. They represent the systemic error as due to the selection of a specific stellar subtraction method.

The implication of these results allow us to draw the following conclusions:

In order to properly derive the kinematics of the ionized gas component the subtraction of the stellar continuum is not crucial, although it is when the real flux of the line has to be recovered; however, the determination of the H$\alpha$ width is quite affected from the choice of the stellar subtraction method. This factor has to be considered in presence of outflow; 2. 2.

Any method (including none) can be used to measure the gas kinematics as long as an additional term ($\Delta\overline{\lambda}$ = $\Delta\overline{\sigma}$ $\sim$ 0.2 $\Delta_{\tt MUSE}$) is added to the error budget. Further analysis with other local sources observed with MUSE, as well as other instruments, should be considered to establish the validity of this equation; 3. 3.

The kinematic results obtained from combining the line fitting results of the H$\alpha$–[NII] complex show a robust determination of the global kinematic properties, as the strongest (and thus more stable) line tends to dominate the result and the most extreme deviations, which could be derived when individual lines are considered, are thus diluted. This is supported by the fact that the [NII] line is not affected by stellar absorption.

Acknowledgements.

We thank the anonymous referee for useful comments and suggestions which helped us to improve the quality and presentation of the manuscript. This paper is based on observations carried out at the European Southern Observatory, Paranal (Chile), Program 60.A-9100(B) with MUSE at Very Large Telescope (VLT). This research made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. EB acknowledges support from the Spanish Ministerio de Economía y Competitividad through the grant AYA2016-79724-C4-1-P as well as the exchange programme ‘Study of Emission–Line Galaxies with Integral–Field Spectroscopy’ (SELGIFS, FP7–PEOPLE–2013–IRSES–612701), funded by the EU through the IRSES scheme. YA is also supported by the Ramón y Cajal programme (RyC–2011–09461; Mineco, Spain). HIM wants to acknowledge the Mexican CONACyT grant (Ciencia Básica) 285080 and PAPIIT–DGAPA–IA104118 (UNAM) project for funding the develop of the Pipe3D code. SFS thanks Conacyt CB–285080, FC–2016–01–1916 and PAPIIT–DGAPA–IN100519 projects.

Appendix A Complementary line fitting results

In this Appendix we present complementary kinematic results which have not been shown in the main text. In Fig. 10 we present the line fitting results obtained when considering the case 2,3 and 4 applied to the H$\alpha$–[NII] complex to the spectrum ‘a’ as an example. This kind of analysis has been applied to the whole sample of spectra to derive the kinematic results discussed in the text.

In Fig. 11 we show the kinematic results ($\lambda_{rest}$, $\sigma_{obs}$) obtained for all spectra when applying the four cases for the five methods. When no stellar subtraction method (Unsubtracted) is applied to the H$\alpha$ line, we derive the smallest value of the velocity dispersions (i.e., $\sigma\sim$1.3 Å) followed by that obtained when applying the STARLIGHT (MILES) method. When considering the pPXF routine we get the highest velocity dispersion (i.e., $\sigma\sim$1.7 Å). The [NII] velocity dispersion values cover a similar range (1.4 – 1.8 Å) and show an analogous trend than that derived for the Balmer line.

We also present in Tabs. 7, 9 and 9 the main results obtained when combining different line fitting approaches with different stellar subtraction methods. In particular, in Tab. 7 the respective mean values for the wavelength and width of the two emission lines are shown for all different cases and methods. In Tabs. 9 and 9 we complement the information given in the previous table considering the mean values of the four cases for each method and vice versa (i.e., mean values of the five methods for each case).

Appendix B Chi–square results

In this section we present the results obtained for each spaxel when computing the chi–square ($\chi^{2}$) values in order to test which stellar continuum subtraction code more faithfully explains the raw observed data (within the errors).

The goodness of the fits has been evaluated by means of the reduced $\chi^{2}$ ($\chi^{2}_{red}$), where the number of degrees of freedom are the number of data points used in the fit, following the equation:

[TABLE]

where $N_{\lambda}$ is the number of observed data points, $O_{i}$ and $M_{i}$ are the raw observed flux data and the total flux model values for the $i$–th point, and $\sigma_{i}$ is the corresponding observed error.

We constrain the derivation of the $\chi^{2}_{red}$ values in the range defined by the lines of our interest, which involves the H$\alpha$–[NII] complex (rest frame wavelength range from 6540 Å to 6600 Å).

In particular, the total model emission in the wavelength range of our interest ($M_{i}$) is given by the combination of the stellar continuum model ($M_{i,\star}$) and the pure gas model ($M_{i,gas}$) emissions. The pure gas emission ($M_{i,gas}$) has been derived as the average of the models obtained when applying the four line fitting cases ($N_{c}$ = 4, see text for details). Thus, each model can be described as follows:

[TABLE]

In Figs. 12 and 13 we compare the data $O_{i}$ with the different models $M_{i}$, derived as described in Eq. 16, for each spaxel. The $\chi$ values as a function of the wavelength ($\lambda$) are also shown for each model.

The reduced chi–square $\chi^{2}_{red}$ values for all spaxels are reported in Tab. 10 and they are also shown in Fig. 14. The best solution (lowest $\chi^{2}_{red}$) as well as the secondary best $\chi^{2}_{red}$ results are highlighted in Tab. 10.

Our results suggest that for the majority of the spectra, especially those characterized by a log([NII]/H$\alpha$) ratio in between –0.49 and $\sim$0.0 (low H$\alpha$/[NII] ratio; spectra from a to f), the best solution is achieved when using the FIT3d method. The pPXF method seems to give good fits only when spectra with higher H$\alpha$/[NII] ratio are considered (i.e., j, k, l). In between these two ranges (i.e., spectra g, h, i) FIT3d and STARLIGHT (MILES and STELIB) methods derive the best fits.

However, the Unsubtracted model does not give good solutions in any case, deriving for the innermost spectra (i.e., from [a] to [d] with low H$\alpha$/[NII] ratio) the most discrepant results (a factor of $\times\sim$2) among all methods considered in this analysis. When considering spectra with higher H$\alpha$/[NII] ratio, the Unsubtracted model can also give reasonable secondary fits in a few cases, although for the majority of the spectra the best candidates for a secondary best $\chi^{2}_{red}$ are FIT3D and STARLIGHT (MILES and STELIB).