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

MgII Surface Brightness profile from MUSE datacubes

In this exercise, you have access to MUSE datacubes of 130 MgII absorbers. These have already been reduced and centered around the relevant source. Your task is:

1. Stack all the datacubes to boost the S/N ratio.
2. Find the optimal wavelength window that maximizes the spatial S/N.
3. Calculate the MgII surface brightness profile.

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'] # Now this is in pkcp
        self.dy = self.header['CD2_2']  # Now this is in pkcp

        #Wavelength sampling is done every 1 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 stack the MUSE cubes centred in the MgII line to produce surface brigthness profiles. The sources are already centred, so we proceed with the stack safely.

#read one file to understand the structure
cube = muse_file('../MUSECubes/MgII_sample2/mgii_cube1.fits')
nx = cube.nx
ny = cube.ny
nz = cube.nz

#Now get all the files:
input_files = glob.glob('../MUSECubes/MgII_sample2/*.fits')

fluxes_stack = np.zeros([nz, nx, ny])
fluxes_var = np.zeros([nz, nx, ny])

for i, input_file in enumerate(input_files):
    cube = muse_file(input_file)       
    
    fluxes_stack += np.nan_to_num(cube.flux)
    fluxes_var += np.nan_to_num(cube.var)
    
extent = [
    -cube.dx*cube.nx/2, cube.dx*cube.nx/2,
    -cube.dy*cube.ny/2, cube.dy*cube.ny/2,
]

We now examine the spectrum of the stack:

#uncertainty associated with stack
spectrum_sum = np.sum(fluxes_stack, axis=(1,2))
spectrum_var = np.sum(fluxes_var, axis=(1,2))
spectrum_sigma = np.sqrt(spectrum_var)

#Find the MgII 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='--')


#A window of 8 A selects the emmision and aboroption well.
delta_l = 8
plt.errorbar(cube.l[kmax-delta_l:kmax+delta_l], spectrum_sum[kmax-delta_l:kmax+delta_l], yerr=spectrum_sigma[kmax-delta_l:kmax+delta_l], capsize=3, color='r', drawstyle='steps-mid')    

plt.xlabel('Wavelength / A')
plt.ylabel('Flux / 10$^{-20}$ erg s$^{-1}$ cm$^{-2}$ $A^{-1}$')
plt.grid()

We now plot the spatial distribution of the stack within the selected wavelength range.

#Stack over all sources
flux_stack = np.sum(fluxes_stack[kmax-delta_l:kmax+delta_l,:,:], axis=0) * cube.dl
flux_var = np.sum(fluxes_var[kmax-delta_l:kmax+delta_l,:,:], axis=0) * cube.dl**2
flux_sigma = np.sqrt(flux_var) 

plt.imshow(flux_stack, extent=extent)
plt.xlabel('Spatial distance [pkpc]')
plt.ylabel('Spatial distance [pkpc]')

Text(0, 0.5, 'Spatial distance [pkpc]')

Next, we measure the profile. There are several way to characterize the profile. We can use, for example, the mean profile, the median, etc.

We first consider the mean profile:

x = -cube.dx*cube.nx/2 + np.arange(0,cube.nx) * cube.dx
y = -cube.dy*cube.ny/2 + np.arange(0,cube.ny) * cube.dy

xv, yv = np.meshgrid(x,y)
r = np.sqrt(xv**2 + yv**2)

#consider logarithmic radial bins
bins = 10**np.linspace(0.3,1.4,15)

#Calculate the mean profile in radial annuli
yhist_counts, xhist = np.histogram(r.flatten(), bins=bins)

#Calculate the weighted sum per radial bin
yhist_weight, xhist = np.histogram(r.flatten(), bins=bins, weights=flux_stack.flatten())

#Calculate the weighted uncertainty per radial bien
yhist_err_weight, xhist = np.histogram(r.flatten(), bins=bins, weights=flux_var.flatten())

#The actual profile
profile = yhist_weight / (np.pi * (xhist[1:]**2-xhist[:-1]**2))
profile[yhist_counts > 0] /= (yhist_counts[yhist_counts >0])

#The uncertainty
profile_sigma = np.sqrt(yhist_err_weight / (np.pi * (xhist[1:]**2-xhist[:-1]**2)))
profile_sigma[yhist_counts>0] /= np.sqrt(yhist_counts[yhist_counts>0])

#Geometric mean
xmed = np.sqrt(xhist[1:]*xhist[:-1])
plt.errorbar(np.log10(xmed), profile, yerr=profile_sigma, capsize=3, color='k', marker='o')

plt.xlabel(r'$\rm log_{10} \left ( R / pkpc \right )$')
plt.ylabel(r'Flux / 10$^{-20}$ erg s$^{-1}$ cm$^{-2}$ kpc$^{-2}$')
plt.grid()

#Note: All the previous lines can be called all at once using scipy.stats.binned statistics. However, it is alway useful to know what we are doing.

Note that the percentiles above do not take into account the intrinsic uncertainty in the data, but the typical dispersion of the data.

As a final step to characterize the uncertainty in our measured profile consist in Bootstraping with replacement from the stack and consider the mean profile and standard deviation of the resulting average. We cab compare the uncertainty of this excersice with the uncertainty obatined previously.

from scipy.stats import binned_statistic

#Now get all the files:
input_files = glob.glob('../MUSECubes/MgII_sample2/*.fits')

fluxes = np.zeros([len(input_files),nz, nx, ny])

for i, input_file in enumerate(input_files):
    cube = muse_file(input_file)   
    fluxes[i,:,:,:] = cube.flux    

N = 1000
mean_profiles = []

#Get N random samples with replacement.
for i in range(N):
    random_int = np.random.randint(0, len(input_files), len(input_files))
    
    flux_stack = np.nansum(fluxes[random_int, kmax-delta_l:kmax+delta_l,:,:], axis=(0,1)) * cube.dl  
    mean_flux, bin_edges, bin_counts = binned_statistic(r.flatten(), flux_stack.flatten(), bins=bins, statistic='mean')
    
    #The actual profile
    profile = mean_flux / (np.pi * (bin_edges[1:]**2-bin_edges[:-1]**2))
    mean_profiles.append(profile)

profile_bootstrap = np.mean(mean_profiles, axis=0)
profile_sigma_bootstrap = np.std(mean_profiles, axis=0)

plt.errorbar(np.log10(xmed), profile_bootstrap, yerr=profile_sigma_bootstrap, capsize=3, color='k', marker='o')
plt.xlabel(r'$\rm log_{10} \left ( R / pkpc \right )$')
plt.ylabel(r'Flux / 10$^{-20}$ erg s$^{-1}$ cm$^{-2}$ kpc$^{-2}$')
plt.grid()

These uncertainties can be added in quadrature to the intrinsic uncertainties so that we report a more realistic uncertainty that considers the inherent uncertainty of the measurements and the typical dispersion of the mean, estimated via Bootstrap with replacement.

However, as the bootstrap error is smaller than the propagated instrumental uncertainties, the instrument's uncertainties are a safer lower bound to the errors.