Could Spectrum1D object also store uncertainty from the spectral_axis?

The wavelength axis in a spectrum might have a variance or uncertainty associated with it, as well as the flux axis. Could we store this within a `Spectrum1d` object?

One use case: when making wavelength solutions from arc lamp lines using a Gaussian Process

The workaround currently would be to use the `meta` OrderedDict that is within `Spectrum1d`, which can be assigned easily enough:
````python
spec = Spectrum1D(spectral_axis = wave_pred * u.angstrom, flux=flux,
uncertainty=StdDevUncertainty(fluxerr),
meta={'spectral_uncertainty':wave_var * u.angstrom} )
````
but `meta` tends to get dropped when manipulating the spectrum (e.g. when doing arithmetic)... this is a bug related to #709, and not my point here except to say `meta` doesn't really work for preserving this extra field yet.

Start by reviewing Spectrum1D's existing spectral_axis and uncertainty handling, then follow the metadata-loss behavior referenced in #709. Use the Gaussian Process wavelength-solution example to clarify how spectral-axis uncertainty should be represented and preserved through arithmetic; the issue names no files or tests.

Written by the indexing model from the issue text.