task: metadata: install: pip install -e .[test]
patch: diff --git a/CHANGES.rst b/CHANGES.rst
index 82cd03d9..165477ae 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -7,6 +7,9 @@ New Features
 Bug Fixes
 ^^^^^^^^^
 
+- Arithmetic operations on ``Spectrum1D`` objects now preserve spectral axis values that
+  were updated by setting redshift or radial velocity. [#1158]
+
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index 7702b58f..5942a139 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -680,6 +680,10 @@ class Spectrum1D(OneDSpectrumMixin, NDCube, NDIOMixin, NDArithmeticMixin):
     def radial_velocity(self, val):
         self.shift_spectrum_to(radial_velocity=val)
 
+    def _return_with_redshift(self, result):
+        result.shift_spectrum_to(redshift=self.redshift)
+        return result
+
     def __add__(self, other):
         if not isinstance(other, (NDCube, u.Quantity)):
             try:
@@ -687,7 +691,7 @@ class Spectrum1D(OneDSpectrumMixin, NDCube, NDIOMixin, NDArithmeticMixin):
             except TypeError:
                 return NotImplemented
 
-        return self.add(other)
+        return self._return_with_redshift(self.add(other))
 
     def __sub__(self, other):
         if not isinstance(other, NDCube):
@@ -696,25 +700,25 @@ class Spectrum1D(OneDSpectrumMixin, NDCube, NDIOMixin, NDArithmeticMixin):
             except TypeError:
                 return NotImplemented
 
-        return self.subtract(other)
+        return self._return_with_redshift(self.subtract(other))
 
     def __mul__(self, other):
         if not isinstance(other, NDCube):
             other = u.Quantity(other)
 
-        return self.multiply(other)
+        return self._return_with_redshift(self.multiply(other))
 
     def __div__(self, other):
         if not isinstance(other, NDCube):
             other = u.Quantity(other)
 
-        return self.divide(other)
+        return self._return_with_redshift(self.divide(other))
 
     def __truediv__(self, other):
         if not isinstance(other, NDCube):
             other = u.Quantity(other)
 
-        return self.divide(other)
+        return self._return_with_redshift(self.divide(other))
 
     __radd__ = __add__
 

problem_statement: Dividing fluxes of  a redshifted spectrum
I load a spectrum and then I shift the spectral_axis to a given redshift value

spec = Spectrum.read('myspectrum.fits')
spec.shift_spectrum_to(redshift=2.0)

then I want to divide the fluxes by a given value (i.e 1+z)

spec /= 3.0

Now the fluxes are divided by 3, but the spectral_axis is back to the rest frame values (and the redshift is back to 0.0).

Is this a bug? Or an intended behaviour?

Thank you for any advice

repo: astropy/specutils
test_cmd: pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
test_patch: diff --git a/specutils/tests/test_arithmetic.py b/specutils/tests/test_arithmetic.py
index b432609b..5f1a698f 100644
--- a/specutils/tests/test_arithmetic.py
+++ b/specutils/tests/test_arithmetic.py
@@ -131,3 +131,12 @@ def test_with_constants(simulated_spectra):
     r_sub_result = 2 - spec
     l_sub_result = -1 * (spec - 2)
     assert_quantity_allclose(r_sub_result.flux, l_sub_result.flux)
+
+
+def test_arithmetic_after_shift(simulated_spectra):
+    spec = simulated_spectra.s1_um_mJy_e1
+    spec.shift_spectrum_to(redshift = 1)
+
+    # Test that doing arithmetic preserves the shifted spectral axis
+    spec *= 2
+    assert_quantity_allclose(spec.spectral_axis, 2*np.linspace(0.4, 1.05, 100)*u.um)