SpectralCoord can't handle pixel units

When I have 2 Spectrum1D objects, which are not wavelength calibrated and have pixels for the spectral unit, I cannot do e.g. addition or subtraction of my spectra. This is an issue for doing things like co-adding spectra, doing sky subtraction, etc.

Here is a MRE:

````python
from specutils import Spectrum1D
import numpy as np
import astropy.units as u
t1 = Spectrum1D(spectral_axis=np.arange(100)*u.pixel, flux=np.random.normal(size=100)*u.adu)
t2 = Spectrum1D(spectral_axis=np.arange(100)*u.pixel, flux=np.random.normal(size=100)*u.adu)
t1-t2
````

which throws this error:
`UnitTypeError: SpectralCoord instances require units equivalent to '(Unit("Hz"), Unit("m"), Unit("J"), Unit("1 / m"), Unit("km / s"))', so cannot set it to 'pix'.`

I believe the workaround here is to do:
````python
t1.subtract(t2, compare_wcs=None)
````

Could we expose this easier to users?

Start with the Spectrum1D subtraction path shown in the minimal reproducible example, especially the compare_wcs=None workaround. Trace how pixel-unit spectral axes are handled during subtraction and determine how that behavior should be exposed without requiring the workaround. Done means two Spectrum1D objects with pixel spectral units can be subtracted or added as described.

Written by the indexing model from the issue text.