Hi!

I work on/with a code that inherits a ton of functionality from specutils Spectrum objects (called muler) and I noticed an unexpected behavior with the Spectrum objects arithmetic recently. When I read in my spectra I immediately barycenter velocity correct and RV shift because I slice high resolution spectra into pretty narrow band regions centered around specific lines to make computation time faster. When I was doing arithmetic I noticed that the wavelengths of two objects I was comparing were shifting in wavelength space which I thought was unusual.

Reading through the code in spectrum.py I noticed that whenever arithmetic is done, ._return_with_redshift is called and shifts the wavelengths of the Spectrum object based on the objects RV (if the spectral axis is not in pixels or undefined units). However, since I BVC and RV shift the spectrum before any arithmetic is done to find the correct line regions this introduces an additional unrequired RV shift that misaligns my objects in wavelength space when I try to compare them.

I guess I was wondering why this is standard behavior for Spectrum objects because it seems like an unusual choice to me to not let the user elect to RV shift the object as a separate function. I was also surprised that this wasn't documented in the specutils Spectrum Arithmetic documentation because I feel like outlining this behavior is important to avoid confusion for users.

I haven't found a way around this to fix it.