Material prepared for the 52nd Saas-Fee Advanced Course on "The circum-galactic medium across cosmic time: an observational and modelling challenge"

by Alejandro Benitez-Llambay & Michele Fumagalli

Contact for questions or comments: alejandro.benitezllambay@unimib.it

Ly-alpha line profiles from MUSE cubes

In this exercise, you have access to MUSE datacubes of bright Ly-alpha emitters. These have already been reduced and centered around the relevant source. Your task is:

1. Read the files and familiarize yourself with the data.
2. Extract the spectrum of the source and its associated uncertainty. You can do this for the entire spatial range or focus on select regions.
3. Find the line profile and transform the spectrum's wavelength into relative velocity with respect to the peak of the Ly-alpha line.
4. Fit a single gaussian profile to the line, considering the intrinsic uncertainties.
5. Fir multiple (n >= 2) gaussian profiles to the line.

BONUS EXERCISE:

6) Estimate the redshift to the sources:

Note:

Consider $f=\displaystyle\sum_{i=0}^{N-1} x_{i}$, in which $x_{i}$ has an associated uncertainty $\sigma_{i}$.

The uncertainty of $f$ can be calculated via:

$\sigma_{f}^2 = \displaystyle\sum_{i=0}^{N-1} \left | \displaystyle\frac{\partial f}{\partial x_{i}}\right |^2 \sigma_{i}^2 = \displaystyle\sum_{i=0}^{N-1} \sigma_{i}^2 = \displaystyle\sum_{i=0}^{N-1} Var(x_{i})$,

where $Var(x_{i})$ is the variace of $x_{i}$.

The uncertainty of $f$ thus becomes:

$\sigma_{f} = \sqrt {\displaystyle\sum_{i=0}^{N-1} Var(x_{i})}$.

The muse datacubes have stored the intrinsic variance in the cell.

Example solution

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from scipy.optimize import curve_fit
import glob

To carry out (1), we write a helper class that can handle the relevant MUSE datacube data

# We define a helper class to handle the datacube
class muse_file():
    def __init__(self, filename):
        current_cube = fits.open(filename)
        #flux in electron counts
        self.flux   = current_cube[0].data 
        #variance as reported by the intrument
        self.var    = current_cube[1].data 
        #header with useful information
        self.header = current_cube[0].header  
        #dimensions of the cube
        self.nz, self.ny, self.nx = self.flux.shape  

        #zero point for wavelength and space.
        self.l0 = self.header["CRVAL3"]
        self.x0 = self.header["CRVAL1"] 
        self.y0 = self.header["CRVAL3"]

        #Spatial sampling in arcsec
        self.dx = self.header['CD2_2'] * 60 * 60
        self.dy = self.header['CD2_2'] * 60 * 60

        #Wavelength sampling is done every 0.125 A
        self.dl = self.header["CD3_3"]

        self.l = self.l0 + np.arange(0,self.nz,1)*self.dl
        self.x = self.x0 + np.arange(0,self.nx,1)*self.dx
        self.y = self.y0 + np.arange(0,self.ny,1)*self.dy

Next, we become familiar with the data (Task 1) and extract the spectrum considering the uncertaintys (Task 2), for the whole sample.

input_files = glob.glob('../MUSECubes/Lya/*.fits')

fig = plt.figure(1,figsize=(10,10), tight_layout=True)
axes = [fig.add_subplot(int(len(input_files)/2)+1,2,i+1) for i in range(len(input_files))]

for input_file, ax in zip(input_files, axes):
    cube = muse_file(input_file)
    
    #Get the mean spectrum within the sampled spatial range:
    spectrum_sum = np.nansum(cube.flux, axis=(1,2))

    #uncertainty associated with the stack
    spectrum_var = np.nansum(cube.var, axis=(1,2))
    spectrum_sigma = np.sqrt(spectrum_var)

    #Plot the extracted spectrum with 1-sigma error bars
    ax.errorbar(cube.l, spectrum_sum, yerr=spectrum_sigma, capsize=3, color='k', drawstyle='steps-mid')    
    ax.grid()
    #ax.set_ylim(-0.5,1.0)
    ax.set_title(input_file)
    
    ax.set_xlabel('Wavelength / A')
    ax.set_ylabel('Flux')
    
plt.show()

in which the units of the y-axis is (10$^{-20}$ erg s$^{-1}$ cm$^{-2}$ $A^{-1}$)

For task 3 we need to find the line profile and transform the spectrum's wavelength into relative velocity with respect to the peak of the Ly-alpha line. To show how this can be done, we consider a simple worked example: ("lae_cube9.fits")

input_file = "../MUSECubes/Lya/lae_cube8.fits"
cube = muse_file(input_file)

#Get the mean spectrum within the sampled spatial range:
spectrum_sum = np.nansum(cube.flux, axis=(1,2))

#uncertainty associated with stack
spectrum_var = np.sum(cube.var, axis=(1,2))
spectrum_sigma = np.sqrt(spectrum_var)

#Find the Lya line:
kmax = np.argmax(spectrum_sum)
lmax = cube.l[kmax]
fmax = spectrum_sum[kmax]

plt.errorbar(cube.l, spectrum_sum, yerr=spectrum_sigma, capsize=3, color='k', drawstyle='steps-mid')    
plt.axvline(cube.l[kmax], ls='--')
plt.xlabel('Wavelength / A')
plt.ylabel('Flux / 10$^{-20}$ erg s$^{-1}$ cm$^{-2}$ $A^{-1}$')
plt.grid()

In this example, finding the flux's maximum does an excellent job of identifying the line. With the position of the line, we can now convert wavelength into velocity relative to the line and try to fit a single Gaussian profile

With the location of the line identified, we now transform the x-axis to velocity, and complete Tasks (4) and (5) by fitting one and two Gaussians to the data:

def gaussian(x, A, x0, sigma):
    return A * np.exp(-0.5 * (x-x0)**2 / sigma**2)

c = 3e5 #km /s
delta_v = c * (cube.l-lmax)/lmax

params, pcov = curve_fit(
    gaussian, delta_v, spectrum_sum, sigma=spectrum_sigma, p0=[1.0, 0.0, 100.0]
)

plt.errorbar(delta_v, spectrum_sum, yerr=spectrum_sigma, capsize=3,  color='k', drawstyle='steps-mid')    
plt.plot(delta_v, gaussian(delta_v, params[0], params[1], params[2]))
plt.xlabel(r'$\Delta V [km/s]$')
plt.ylabel('Flux / 10$^{-20}$ erg s$^{-1}$ cm$^{-2}$ $A^{-1}$')
plt.grid()

print("Single gaussian parameters:")

print(r" A = %.3f +- %.3f"%(params[0], np.sqrt(pcov[0][0])))
print(r" v0 = %.3f +- %.3f"%(params[1], np.sqrt(pcov[1][1])))
print(r" sigma= %.3f +- %.3f"%(params[2], np.sqrt(pcov[2][2])))

Single gaussian parameters:
 A = 703.180 +- 55.022
 v0 = 19.577 +- 19.746
 sigma= 217.135 +- 19.607

A single Gaussian profile does a god job of fitting the line. Still, there is a small excess in the tails

To fit two Gaussians, we shall consider a different example (lae_cube15.fits)

# We repeat the same excersice but now for two gaussians:

#This functions reuses the Gaussian defined above
def bigaussian(x, A1, x1, sigma1, A2, x2, sigma2):
    return gaussian(x, A1, x1, sigma1) + gaussian(x, A2, x2, sigma2)

input_file = "../MUSECubes/Lya/lae_cube8.fits"
cube = muse_file(input_file)

#Get the stack:
spectrum_sum = np.nansum(cube.flux, axis=(1,2))

#uncertainty
spectrum_var = np.nansum(cube.var, axis=(1,2))
spectrum_sigma = np.sqrt(spectrum_var)

#Find the Lya line:
kmax = np.argmax(spectrum_sum)
lmax = cube.l[kmax]
fmax = spectrum_sum[kmax]

c = 3e5 #km /s
delta_v = c * (cube.l-lmax)/lmax

params, pcov = curve_fit(
    bigaussian, delta_v, spectrum_sum, sigma=spectrum_sigma, p0=[800.0, 0.0, 100.0,800,-100.0,100.0]
)


plt.errorbar(delta_v, spectrum_sum, yerr=spectrum_sigma, capsize=3, color='k', drawstyle='steps-mid')    

plt.plot(delta_v, gaussian(delta_v, params[0], params[1], params[2]), ls='--')
plt.plot(delta_v, gaussian(delta_v, params[3], params[4], params[5]), ls='--')
plt.plot(delta_v, bigaussian(delta_v, params[0], params[1], params[2], params[3], params[4], params[5]), color='k')

plt.xlabel(r'$\Delta V [km/s]$')
plt.ylabel('Flux / 10$^{-20}$ erg s$^{-1}$ cm$^{-2}$ $A^{-1}$')
plt.grid()


print("Double gaussian parameters:")
print(r" (A1,A2) = (%.3f +- %.3f,%.3f +- %.3f)"%(params[0], np.sqrt(pcov[0][0]), params[3], np.sqrt(pcov[3][3])))
print(r" (v01,v02) = (%.3f +- %.3f,%.3f +- %.3f)"%(params[1], np.sqrt(pcov[1][1]), params[3], np.sqrt(pcov[4][4])))
print(r" (sigma1,sigma2) = (%.3f +- %.3f,%.3f +- %.3f)"%(params[2], np.sqrt(pcov[2][2]), params[5], np.sqrt(pcov[5][5])))

Double gaussian parameters:
 (A1,A2) = (539.771 +- 127.395,308.368 +- 110.311)
 (v01,v02) = (67.263 +- 22.645,308.368 +- 100.849)
 (sigma1,sigma2) = (119.545 +- 29.850,327.412 +- 63.405)

To solve the bonus excercise, we need to bring the observed spectrum to vacuum, for which we know the rest-fram Lya wavelength, 1215.67 A. The air to vaccum conversion is done using the following function:

def air2vac(wave):
    sigma2 = (1e4/wave)**2.    
    fact = 1.+5.792105e-2/(238.0185-sigma2)+1.67917e-3/(57.362-sigma2)
    return wave*fact

input_file = "../MUSECubes/Lya/lae_cube8.fits"
cube = muse_file(input_file)

#Get the mean spectrum within the sampled spatial range:
spectrum_sum = np.nansum(cube.flux, axis=(1,2))
l_vac = air2vac(cube.l)

#Find the Lya line:
kmax = np.argmax(spectrum_sum)
lmax = l_vac[kmax]

lya = 1215.67

c = 3e5 #km/s
redshift = (lmax-lya)/lya

print("Redshift of %s is z=%.2f"%(input_file, redshift))

Redshift of ../MUSECubes/Lya/lae_cube8.fits is z=3.13