Covariance between spectral elements

Is there any interest in supporting covariance matrices between different spectral elements? In `sncosmo`, we're looking at how we can fit models of supernovae to spectral data, and are interested in basing this off of `specutils`. One challenge is that we often have covariance information between the different spectral elements that needs to be included in the fit. A discussion of this can be found in sncosmo/sncosmo-eps#3.

At the most basic level, for a Spectrum1D object this involves simply storing a matrix of the covariance between the different spectral elements. It would be nice if the different ways of manipulating spectra could propagate the covariance where appropriate (e.g. resampling). A lot of those operations introduce covariance between spectral elements even if the spectral elements were independent before. I would be happy to help implement this if there is interest.

Start by reviewing the Spectrum1D object and the discussion in sncosmo/sncosmo-eps#3. Determine how covariance between spectral elements should be stored and which spectrum-manipulation operations, especially resampling, must propagate it. Done means the representation and propagation behavior are agreed and implemented with coverage for the relevant operations.

Written by the indexing model from the issue text.