Directly compare subpixel shifts between two spectra — and get believable errors
I think using cross-correlation and interpolating the peak would work fine.  As described in Is up-sampling prior to cross-correlation useless?, interpolating or upsampling before the cross-correlation doesn't actually get you any more information.  The information about the sub-sample peak is contained in the samples around it. You just need to extract it with minimal error.  I gathered some notes here.

The simplest method is quadratic/parabolic interpolation, which I have a Python example of here.  It's supposedly exact if your spectrum is based on a Gaussian window, or if the peak happens to fall exactly on the midpoint between samples, but otherwise has some error.  So in your case you probably want to use something better.

Here is a list of more complicated, but more accurate estimators.  "Of the above methods, Quinn's second estimator has the least RMS error."

I don't know the math, but this paper says that their parabolic interpolation has theoretical accuracy of 5% of the width of an FFT bin.

Using FFT interpolation on the cross-correlation output does not have any bias error, so that's the best if you want really good accuracy.  If you need to balance accuracy and computation speed, it's recommended to do some FFT interpolation and then follow it with one of the other estimators to get a "good enough" result.

This just uses the parabolic fit, but it outputs the right value for the offset if the noise is low:

def parabolic_polyfit(f, x, n):
    a, b, c = polyfit(arange(x-n//2, x+n//2+1), f[x-n//2:x+n//2+1], 2)
    xv = -0.5 * b/a
    yv = a * xv**2 + b * xv + c

    return (xv, yv)

arraysize = 1001
fluxlevel = 100.0
noise = 0.3 # 2.0 is too noisy for sub-sample accuracy
signal_std = 15.0
signal_depth = 40.0
gaussian = lambda x: np.exp(-(mu-x)**2/ (2 * signal_std))
mu = 500.1
a_flux = np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth)
mu = 500.8
b_flux = np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth)

a_flux -= np.mean(a_flux)
b_flux -= np.mean(b_flux)

corr = ss.fftconvolve(b_flux, a_flux[::-1])

peak = np.argmax(corr)
px, py = parabolic_polyfit(corr, peak, 13)

px = px - (len(a_flux) - 1)
print px





The noise in your sample produces results that vary by more than a whole sample, so I reduced it.  Fitting the curve using more of the points of the peak helps tighten the estimate somewhat, but I'm not sure if that's statistically valid, and it actually makes the estimate worse for the lower-noise situation.

With noise = 0.2 and 3-point fit, it gives values like 0.398 and 0.402 for offset = 0.4.

With noise = 2.0 and 13-point fit, it gives values like 0.156 and 0.595 for offset = 0.4.