Extra Galactic Background Light

Generating the EBL

You will need ProSpect v0.8.6 for this vignette to fully run.

Load the libraries we need:

library(ProSpect)
library(celestial)

## Loading required package: RANN

## Loading required package: NISTunits

## Loading required package: pracma

library(foreach)

Encode the Madua and Dickenson 2014 CSFH function:

MD14func = function(z, norm=0.00945){ #norm here is MD14 0.015*0.63 (latter Salp to Chab Driver 2016 correction)
  norm * ((1+z)^2.7) / ((1+((1+z)/2.9)^5.6))
}

Create a function to generate the EBL:

EBL = function(massfunc_z = MD14func, norm=0.00945, Zstart=1e-4, Zfinal=0.02,
               zvec = 10^seq(-3,1, by=0.1), SFRvec=rep(0.01,length(zvec)),
               filters=c('FUV_GALEX', 'NUV_GALEX', 'u_SDSS', 'g_SDSS', 'r_SDSS', 'i_SDSS',
               'Z_VISTA', 'Y_VISTA', 'J_VISTA', 'H_VISTA', 'K_VISTA', 'W1_WISE' , 'W2_WISE',
               'W3_WISE', 'W4_WISE', 'P100_Herschel', 'P160_Herschel', 'S250_Herschel' ,
               'S350_Herschel', 'S500_Herschel'), out='Photom', ...){
  
  data(Dale_NormTot)
  data(BC03lr)

  agevec = cosdistTravelTime(zvec, ref='planck') * 1e9 # Create age vector
  
  if(!is.null(massfunc_z)){ # You can either pass in an SFR function, or an SFRvec of the values at redshifts given by zvec
    SFRvec = massfunc_z(zvec, norm=norm)
  }
  
  covolvec = cosdistCoVol(zvec, OmegaM = 0.3, OmegaL=0.7, H0=70) * 1e9 #  x1e9 to get to Mpc^3 (output is Gpc^3)
  volweights = (c(0,diff(covolvec)) + c(diff(covolvec),0))/2 # This evens out the volume weights around bin edges
  filtout = foreach(i = filters)%do%{approxfun(getfilt(i))} # Generate our filter sets for quick generation
  
  tempSFH = approxfun(agevec, SFRvec, yleft = 0, yright = 0) # Generate our temporary SFH function to create our ZH
  Zvec = Zfunc_massmap_box(agevec, massfunc=tempSFH, Zstart=Zstart, Zfinal=Zfinal) # Create out closed box ZH
  
  fluxz = foreach(i = 1:length(agevec), .combine='rbind')%do%{
    tempSFH = approxfun(agevec - agevec[i], SFRvec, yleft = 0, yright = 0) # Create a SFH from the current time to the end of the Universe
    tempZH = function(x,...){approxfun(agevec - agevec[i], Zvec, yleft = 0, yright = 0)(x)} # Create a ZH from the current time to the end of the Universe
    
    if(out=='Photom'){
      ProSpectSED(SFH = SFHfunc, massfunc=tempSFH, z = zvec[i], agemax = max(agevec) - agevec[i], filtout=filtout, # Run ProSpect SED for our target tempSFH and tempZH
                Dale=Dale_NormTot, speclib=BC03lr, OmegaM = 0.3, OmegaL=0.7, H0=70, Z=tempZH, returnall=FALSE, ...)
    }else if(out=='Spec'){
      temp=ProSpectSED(SFH = SFHfunc, massfunc=tempSFH, z = zvec[i], agemax = max(agevec) - agevec[i], filters=NULL, # Run ProSpect SED for our target tempSFH and tempZH
                Dale=Dale_NormTot, speclib=BC03lr, OmegaM = 0.3, OmegaL=0.7, H0=70, Z=tempZH, returnall=TRUE, ...)$FinalFlux
      temp=approxfun(temp)
      temp(10^seq(2, 7, by=0.01))
    }
  }
  
  fluxz = fluxz * volweights # Weight by volume elements
  flux = colSums(fluxz) # Sum epochs
  flux = flux * 1e-26 # x1e-26 to get to W/m2/Hz (output is Jansky)
  flux = flux * 1e9 # to get to nW/m2/Hz
  flux = flux / (4*pi) # to get to nW/m2/Hz/sr
  
  if(out=='Photom'){
    tempcen = cenwave[cenwave$filter %in% filters,'cenwave']
    tempfreq =  299792458/(tempcen/1e10) # Convert central wavelength from wavelength to freq
    flux = flux * tempfreq # Scale to nW/m2/sr
    output=data.frame(wave=tempcen, flux=flux)
  }else if(out=='Spec'){
    tempwave = 10^seq(2, 7, by=0.01)
    tempfreq =  299792458/(tempwave/1e10) # Convert central wavelength from wavelength to freq
    flux = flux * tempfreq # Scale to nW/m2/sr
    output=data.frame(tempwave, flux=flux)
  }
  return(invisible(output))
}

The wavelength is in Angstroms and the flux is nW/m2/sr (as per Driver et al 2016).

Now we can generate some EBL models (need to up the default zvec resolution to get this looking smooth though):

magplot(EBL(out='Spec', zvec = 10^seq(-3,1, by=0.01)), xlab='Wave / Ang', ylab='nW/m^2/sr',
        ylim=c(0,13), log='x', type='l', grid=TRUE)
lines(EBL(norm=0.008, Zfinal=0.02, out='Spec', zvec = 10^seq(-3,1, by=0.01)), col='red')
lines(EBL(norm=0.008, Zfinal=0.01, out='Spec', zvec = 10^seq(-3,1, by=0.01)), col='blue')
lines(EBL(norm=0.008, Zstart=0.02, Zfinal=0.02, out='Spec', zvec = 10^seq(-3,1, by=0.01)), col='darkgreen')
lines(EBL(norm=0.008, Zstart=0.01, Zfinal=0.01, out='Spec', zvec = 10^seq(-3,1, by=0.01)), col='brown')

Fitting the EBL

Can we fit the EBL? Let’s try with 5% error, where we just use optical and NIR filters:

set.seed(1)
filters=c('u_SDSS', 'g_SDSS', 'r_SDSS', 'i_SDSS', 'Z_VISTA', 'Y_VISTA', 'J_VISTA',
          'H_VISTA', 'K_VISTA')
EBLin = EBL(norm=0.008, Zfinal=0.01, filters=filters)
EBLin$fluxerr = EBLin$flux*0.05
EBLin$flux = rnorm(9, mean=EBLin$flux, sd=EBLin$fluxerr)

Take a look at our mock EBL:

magplot(EBLin[,1:2], ylim=c(0,13), xlab='Wave / Ang', ylab='nW/m^2/sr', grid=TRUE)
magerr(EBLin$wave, EBLin$flux, ylo=EBLin$fluxerr)

Here we will fit for the normalisation and the final metallicity in our closed box model (notice we do the fit in log space):

# Make likelihood function:
EBLlike=function(par, Data){
  -sum(dnorm(x = EBL(norm=10^par[1], Zfinal=10^par[2], filters=filters)$flux,
             mean = Data$flux, sd = Data$fluxerr, log = TRUE))
}
# Fit it (this takes about a minute):
EBLfit = optim(par=c(-2,-1.7), EBLlike, Data=EBLin, method = 'L-BFGS-B',
               lower=c(-3,-3), upper=c(-1,-1), hessian = TRUE)

And now we can check the outputs. Here I basically get back the model we put in:

print(10^EBLfit$par) # model parameters, I see 0.008235409, 0.012355299
print(sqrt(diag(solve(EBLfit$hessian)))) # dex errors, I see 0.01309136, 0.11251226

Compare our fitted model to the observations:

magplot(EBLin[,1:2], ylim=c(0,13), xlab='Wave / Ang', ylab='nW/m^2/sr', grid=TRUE)
magerr(EBLin$wave, EBLin$flux, ylo=EBLin$fluxerr)
lines(EBL(norm=10^EBLfit$par[1], Zfinal=10^EBLfit$par[2]), col='red')

SFH Normalisation

Looking at the EBLs generated above, it is clear that even for very different metallicity evolutions, the K band EBL flux is a very good estimator of the SFH normalisation.

We can check this:

set.seed(666)
norm_ran = runif(20, 0.007, 0.015)
Zstart_ran = runif(20, 1e-4, 0.04)
Zfinal_ran = runif(20, Zstart_ran, 0.04)

Generate the EBLs:

Knorm = foreach(i=1:20, .combine='c')%do%{EBL(norm=norm_ran[i], Zstart = Zstart_ran[i],
                                            Zfinal=Zfinal_ran[i], filters='K_VISTA')[1,2]}

Find the best fit line:

KvEBL = lm(Knorm ~ norm_ran)
print(KvEBL)

## 
## Call:
## lm(formula = Knorm ~ norm_ran)
## 
## Coefficients:
## (Intercept)     norm_ran  
##    -0.06515   1066.69206

Plot the K band EBL flux versus the input normalisation (for very different ZH remember):

magplot(norm_ran, Knorm, xlab='SFH Norm (Msol/Yr/Mpc^3)', ylab='K band EBL flux (nW/m2^2/sr)',
        grid=TRUE)
abline(KvEBL, col='red')