Handle reading/writing of covariance matrices

[kelle]

Some spectra have relevant covariance matrices which should be stored in the FITS file.

• First, figure out how to get the covariance matrices into a FITS HDU. (Read/Write covariance matrixes as FITS files astropy#16648)
• Then, make a specutils reader/writer which reads/writes a spectrum and a covaraince matrix.

from @kbwestfall:
I would envision it as two different tables in 1 file. The current I/O basically works directly with HDUs. So you could have the main code in core (maybe in NDData actually) that constructs the HDUs that then get added to the spectrum HDU(s).