Radial velocity from Gaussian fit not as expected

[Knusper]

According to https://lime-stable.readthedocs.io/en/latest/outputs/outputs1_measurements.html the v_r output parameter contains the radial velocity (note that there are brackets missing in the definition).

However, the following puzzles me, since I don't get the same radial velocity when using the manual formula and the lime output (the pickle dictionary with spectrum, variance, and wavelength axis is attached).

from astropy import constants as cons
import numpy as np
import pickle
import lime

c = cons.c.to('km/s').value

with open('em_line_fit_debug.pickle', 'rb') as f:
    in_dic = pickle.load(f)

xax0 = in_dic['xax0']  # wavelength axis (rest-frame / air wavelengths)
ln_spec_in = in_dic['ln_spec_in']  # input spectrum (continuum subtracted)
ln_spec_var = in_dic['ln_spec_var']  # variance

# lime bands (for test just Balmer-series & [OIII] doublet
particle_list = ('H1', 'O3')
bands = lime.line_bands(wave_intvl=(xax0[0], xax0[-1]),
                        particle_list=particle_list,
                        units_wave='Angstrom', decimals=None, vacuum=False)

# lime spectrum
lspec = lime.Spectrum(input_wave=xax0, input_flux=ln_spec_in, input_err=np.sqrt(ln_spec_var),
                      units_flux='1e-20*FLAM', redshift=0)
lspec.fit.continuum(degree_list=[3,5,7], emis_threshold=[3,3,3], plot_steps=False)
lspec.line_detection(bands, sigma_threshold=5, plot_steps=False)
lspec.fit.frame(bands=bands)

# check Balmer alpha line fit
l = lspec.frame.loc[['O3_5007A']]

# radial velocity from LiME
v_r_lime = l.v_r.iloc[0]
v_r_manual = (l.center.iloc[0] / l.wavelength.iloc[0] - 1 ) * c 

print('\n')
print('v_r_lime:' + str(v_r_lime))
print('v_r_manual:' + str(v_r_manual))

The output is:

Line fitting progress:
[==========] 100% of 4 lines (H1_6563A)

v_r_lime:-18.850342740006468
v_r_manual:6.74700738846691

Now I would have expected that v_r_manual and v_r_lime are the same - v_r_manual is in this example also the value that appears correct given a previous result on this spectrum.

What is going on here?

(Btw. is there a way to turn of the "line fitting progress" message...)

em_line_fit_debug.pickle.gz

[Vital-Fernandez]

Thank you @Knusper very much for your comment. I have checked your code and there is indeed a difference.

This has to do with the reference wavelength used to calculate the radial velocity. In the case of a single lines, it is the peak wavelength. In your code this would be:

v_r_manual = (l.peak_wavelength.iloc[0] / l.wavelength.iloc[0] - 1 ) * c

where l.peak_wavelength is the wavelength of the maximum flux pixel on the line band, which depending on the spectral resolution may be different from the central wavelength.

In complex blended profiles, specially in IFU data sets, the intensity of the peak may change considerably due to ionization conditions rather than kinematic... which may make it misleading to keep track of the kinematics using the peak wavelength... so Instead I use the theoretical leading line wavelength times (1 + z_input).

I would rather have the same logic for every line type... but I believe most users using longslist spectra prefer to use the peak wavelength and remove the redshift dependence...

Do you think the documentation should be more clear? It maybe good to add an sketch....

[Knusper]

Thank you @Knusper very much for your comment. I have checked your code and there is indeed a difference.

Thanks for checking the code and I excuse very much for the late reply.

This has to do with the reference wavelength used to calculate the radial velocity. In the case of a single lines, it is the peak wavelength. In your code this would be:

v_r_manual = (l.peak_wavelength.iloc[0] / l.wavelength.iloc[0] - 1 ) * c

where l.peak_wavelength is the wavelength of the maximum flux pixel on the line band, which depending on the spectral resolution may be different from the central wavelength.

I don't find peak_wavelength in the documentation. You mean peak_wave?

In complex blended profiles, specially in IFU data sets, the intensity of the peak may change considerably due to ionization conditions rather than kinematic... which may make it misleading to keep track of the kinematics using the peak wavelength...

I think that even for a single line the peak is a very coarse tracer of the kinematics. This can never be as accurate as the v_r inferred from a fit.

so Instead I use the theoretical leading line wavelength times (1 + z_input).

Just to clarify, by theoretical leading line wavelength you refer to the first component (i.e. the one that should have the lowest dispersion in the fitting configuration).

I would rather have the same logic for every line type... but I believe most users using longslist spectra prefer to use the peak wavelength and remove the redshift dependence...

I doubt that anybody prefers the peak wavelength. Moreover, the documentation states that v_r is a Gaussian property.

Do you think the documentation should be more clear? It maybe good to add an sketch....

While I can compute the desired v_r manually with the code as posted in my report, I think it should always be the one that is derived from the Gaussian fit. If, for whatever reason, a v_r based on the integrated property peak_wave is desired, then this should either be an optional parameter, or the user should calculate it as given above.