Arbitrary emission lines strengths with low velocity dispersion at a given spectral resolution

[renskesmit]

Hi Adam,

BagPipes is great - thanks for making it so easily accessible! I have an issue that I can't quite figure out though and I hope you can help me. I took a screenshot of a test notebook to make it clear (sorry I would have attached it but it won't let me). When I generate a (very simple) template with emission lines and I let the spectra resolution be very low, there is a weird 'stochasticity' for when the lines actually show up - in the screen shot, a difference of one Angstrom in spectral resolution changes the strength of Lya and OII completely! The problem goes away if I go to higher resolution or if I take a larger velocity dispersion. With other words, if the velocity dispersion is low compared the the spectral resolution of the template, the emission line implementation becomes 'unstable' if you will - do you know why?

In this case I am trying to use BagPipes on simulated JWST/NIRSpec R=100 spectra. These are large files and in order to keep the runtime down I am re-binning quite strongly (for now) so there will be many sources that have low velocity dispersion compared to the resolution. I noticed that Lya basically never shows up in the fitted templates, even if it does in higher resolution templates I generated, which led me to investigate the particular issue above.

I hope you can help.

Thanks, Renske Smit

[ACCarnall]

Hi Renske,

Thanks for getting in touch, it's great to hear you've been getting on well trying out the code!

I suspect this is a different manifestation of another possible bug someone else recently reported to me about the stellar continuum doing strange things when working with low velocity dispersions. I'll look into this properly over the next couple of days and try to come up with an answer/fix.

Cheers, Adam

[ACCarnall]

Hi Renske,

Unfortunately not quite as simple a solution here as I'd hoped, but I do think I understand what's going on.

Firstly, I was able to produce the effect you're seeing without implementing the velocity dispersion option in the code at all (and for absorption lines without including the nebular component), so the main disappearing emission line problem isn't actually connected to this at all. To get a better idea of what's happening, it's best to look at the internal spectral model bagpipes generates, which is then resampled to whatever wavelength grid you specify in spec_wavs. This is stored under model.spectrum_full, with corresponding (rest-frame) wavelengths under model.wavelengths. I've shown this in black below, with the 100AA sampling overplotted in blue, and the 99AA sampling in red.

The basic reason this is happening is that bagpipes resamples the internal (black) model to your chosen wavelength grid using the numpy.interp function, which uses some sort of highly optimised piecewise linear interpolation algorithm. I got quite concerned about this at various points as I was developing the code, mostly because there's no requirement for this algorithm to conserve flux. I even went as far as writing my own flux-conserving spectral resampling algorithm, spectres. Even so, I ended up using numpy.interp in bagpipes, basically because it's faster and makes no difference in the vast majority of circumstances. It looks like you've managed to find a circumstance where that's not the case! It seems that if you have a very sharp change in the function you're interpolating, and you're undersampling, something like an emission line can get completely missed by the interpolation algorithm, at least if it's not close to one of the new sample points.

So, what to do about this. Well, I'm going to have to have a think about whether this kind of thing is likely to be a bigger problem than I thought before you pointed this out. In the short term however, I'd recommend you just re-implement my spectres algorithm in your version of the code. This is pretty simple, you just have to go to lines 387-389 in models/model_galaxy.py and replace the fluxes= command with

fluxes = spectres.spectres(self.spec_wavs,
                           redshifted_wavs,
                           spectrum, fill=0)

You'll also have to download spectres via pip and put the import statement at the top of the file. This gives behaviour that's hopefully a bit more like what you'd expect:

It's probably worth mentioning that I believe there is still a problem with the stellar continuum at low velocity dispersions, which you might also fall foul of. If you're working with low resolution spectra like this then there's probably no risk of you resolving velocity dispersion in your galaxies, so you can probably safely get around this problem by either a.) just not using the veldisp option at all, or b.) setting veldisp to your best guess at the velocity resolution of your instrument. I know some of the JWST grisms have quite strongly varying resolution as a function of wavelength, sadly there's currently no way of incorporating this into models generated by bagpipes, but it is something I'm interested in working on in the future.

I hope that all makes some sort of sense, all the spectral sampling/resolution/velocity dispersion stuff can really become a major headache when doing full spectral fitting. If you'd like to talk things through in more detail I'd be happy to arrange a meeting sometime, just let me know!

Cheers, Adam

[gbrammer]

Hi @ACCarnall ,

I just saw this go by and thought I'd mention that there is a fast cythonized version of a flux-conserving interpolator in grizli. Your spectres also looks nice, and you could consider cythonizing it. There was really a big speedup w.r.t. pure python for the grizli function.

[renskesmit]

Hi Adam,

Thanks again for your detailed answer and providing those solutions! I'll hopefully have time to start trying them out and will let you know if that works!

Cheers, Renske

[ACCarnall]

Thanks Gabe! I had a vague look into this a few years back, but never ended up with a reliably faster version. Do you find people can get cython installed seamlessly via pip, or do you get issues with non-python dependencies?

No problem Renske, hope it works out.

[gbrammer]

I haven't heard of any problems, anyway. It compiles and tests without any trouble in python 3.6, 3.7, 3.8 in the github actions environment. You get some extra speedups with the @cython.boundscheck and @cython.cdivision decorators, but I guess that's a bit risky.

[ACCarnall]

I've now made spectres the default option for resampling in the model_galaxy class from v1.0.2 of the code. Still to think about are how the code initially resamples the raw stellar grids (though I think this is probably fine) and options for cython, numba etc to speed the code up.