ProGeny: The Children of Stars

Aaron Robotham

2025-02-12

Load what we will need:

library(ProGeny)
## Isochrones: Download with progenyIsoDownload. Load with read.fst
## Atmospheres: Download with progenyAtmosDownload. Load with progenyAtmosLoad
## Interpolate grids with progenyInterpGrid_All
## Update the best match between Isochrones and Atmospheres with progenyInterpBest
## Generate your SSP using the input above with progenyMakeSSP, e.g. minimally:
## 
## Iso  = read.fst('path/to/Iso.fst', as.data.table=TRUE)
## Spec_combine = progenyAtmosLoad('path/to/atmos/')
## Interp_combine = progenyInterpGrid_All(Iso=Iso, Spec_combine=Spec_combine)
## Iso = progenyInterpBest(Iso=Iso, Interp_combine=Interp_combine)
## progenyMakeSSP(Iso=Iso, Spec_combine=Spec_combine, Interp_combine=Interp_combine)
library(ProSpect)
## Loading required package: ProSpectData
## Loading required package: Rcpp
## 
## Attaching package: 'ProSpect'
## The following objects are masked from 'package:ProGeny':
## 
##     binlims, IMF_Chabrier, IMF_Kroupa, IMF_Salpeter
library(data.table)
library(Rfits)
library(fst)
library(magicaxis)

The main thing ProGeny does (and the reason it was designed at all) is to create Simple Stellar Population (SSP) models that can directly load into ProSpect. This is useful because it gives us full control over what isochrones and stellar atmospheres we use, and how we map them together. We can also control our initial mass function (IMF).

The basic path to doing this is to firstly load our desired isochrones (MIST, in this case). On my machine this is at:

MistIso  = read.fst('~/Dropbox (Personal)/ProGeny_isochrone/MistIso.fst', as.data.table=TRUE)
#ParsecIso  = read.fst('~/Dropbox (Personal)/ProGeny_isochrone/ParsecIso.fst') #example if using PARSEC
#BastiIso  = read.fst('~/Dropbox (Personal)/ProGeny_isochrone/BastiIso_New.fst') #example if using BaSTI

If you do not have one already, you can get them using the progenyIsoDownload function. Currently we provide MIST (preferred, for reasons we shall get to), PARSEC and BaSTI as base isochrones.

Here we load four different base spectra options: C3K, PHOENIX (Husser) MILES and BaSeL (mostly styled as BASEL below). To get these locally you will need to run progenyAtmosDownload. The below might take about a minute to run (depending on the speed of your HDD or SSD).

Spec_C3K = progenyAtmosLoad('~/Dropbox (Personal)/ProGeny_atmos/', base='combine_C3K_conroy', wavegrid=BC03lr$Wave)
Spec_PHOENIX = progenyAtmosLoad('~/Dropbox (Personal)/ProGeny_atmos/', base='combine_PHOENIX_husser', wavegrid=BC03lr$Wave)
Spec_MILES = progenyAtmosLoad('~/Dropbox (Personal)/ProGeny_atmos/', base='combine_MILES_vazdekis', wavegrid=BC03lr$Wave)
Spec_BASEL = progenyAtmosLoad('~/Dropbox (Personal)/ProGeny_atmos/', base='combine_BASEL_wlbc', wavegrid=BC03lr$Wave)

Below we interpolate our spectra onto the isochrone grid of Teff, logG and logZ. The below might take a minute or two to run (depending on your CPU speed). For this example we will use MIST, but keen users should try switching to alternative isochrones.

Interp_PHOENIX = progenyInterpGrid_All(Iso = MistIso, Spec_combine = Spec_PHOENIX)
Interp_C3K = progenyInterpGrid_All(Iso = MistIso, Spec_combine = Spec_C3K)
Interp_MILES = progenyInterpGrid_All(Iso = MistIso, Spec_combine = Spec_MILES)
Interp_BASEL = progenyInterpGrid_All(Iso = MistIso, Spec_combine = Spec_BASEL)

We need to update the isochrones with the origin of the ‘best’ spectra (since we can only have one per isochrone element). The below should be very fast. Make sure you pass back in the appropriate isochrone (MIST in this case). Note also you will need to change the arguments for label_AGB and label_white if switching isochrones (since different isochrones use different labelling systems, as suggested by the comments).

Iso_C3K = progenyInterpBest(Iso = MistIso,
                        Interp_combine = Interp_C3K,
                        label_AGB = 4:5, #7:8 [PC] / 4:5 [Mist]
                        label_white = 6:8 # 9 [PC] / 6:8 [Mist]
                        )

Iso_PHOENIX = progenyInterpBest(Iso = MistIso,
                        Interp_combine = Interp_PHOENIX,
                        label_AGB = 4:5, #7:8 [PC] / 4:5 [Mist]
                        label_white = 6:8 # 9 [PC] / 6:8 [Mist]
                        )

Iso_MILES = progenyInterpBest(Iso = MistIso,
                        Interp_combine = Interp_MILES,
                        label_AGB = 4:5, #7:8 [PC] / 4:5 [Mist]
                        label_white = 6:8 # 9 [PC] / 6:8 [Mist]
                        )

Iso_BASEL = progenyInterpBest(Iso = MistIso,
                        Interp_combine = Interp_BASEL,
                        label_AGB = 4:5, #7:8 [PC] / 4:5 [Mist]
                        label_white = 6:8 # 9 [PC] / 6:8 [Mist]
                        )

And finally we can put these steps together to generate our SSPs with our desired target IMF (Chabrier in this case). The below will take a good few minutes to run (depending on your CPU speed and number of cores).

PG_Ch_C3K = progenyMakeSSP(Iso = Iso_C3K,
                               IMFfunc = IMF_Chabrier,
                               masslow = 0.1,
                               massmax = 100,
                               Spec_combine = Spec_C3K,
                               Interp_combine = Interp_C3K
               )
## Generating spec grids:
## Generating evo grids:
PG_Ch_PHOENIX = progenyMakeSSP(Iso = Iso_PHOENIX,
                               IMFfunc = IMF_Chabrier, #but perhaps try IMF_Kroupa or IMF_kroupa
                               masslow = 0.1,
                               massmax = 100,
                               Spec_combine = Spec_PHOENIX,
                               Interp_combine = Interp_PHOENIX
               )
## Generating spec grids:
## Generating evo grids:
PG_Ch_MILES = progenyMakeSSP(Iso = Iso_MILES,
                               IMFfunc = IMF_Chabrier,
                               masslow = 0.1,
                               massmax = 100,
                               Spec_combine = Spec_MILES,
                               Interp_combine = Interp_MILES
               )
## Generating spec grids:
## Generating evo grids:
PG_Ch_BASEL = progenyMakeSSP(Iso = Iso_BASEL,
                               IMFfunc = IMF_Chabrier,
                               masslow = 0.1,
                               massmax = 100,
                               Spec_combine = Spec_BASEL,
                               Interp_combine = Interp_BASEL
               )
## Generating spec grids:
## Generating evo grids:

We can check this with our ProSpect speclib_check function:

speclib_check(PG_Ch_C3K) #just checking this one since it our fiducial SSP

Hopefully it all passes! I won’t do this for all of the SSPs here though (too much clutter).

Spectral Outputs

Now for some plots. Below we can compare a specific target metallicity (Z) and age (in log(Age/Years) units) for fixed isochrone (MIST) and IMF (Chabrier), so we only vary the stellar spectral libraries used (PHOENIX, C3K, MILES and BaSeL). Here we loop around the metallicities available in BC03, and ages 105 – 1010 Years in 1 dex spacing:

data(BC03lr)

for(Z in BC03lr$Z){
  for(logAge in 7:10){
    temp_PHOENIX = cbind(PG_Ch_PHOENIX$Wave, PG_Ch_PHOENIX$Zspec[[min(which(PG_Ch_PHOENIX$Z >= Z))]][min(which(PG_Ch_PHOENIX$Age >= 10^logAge)),])
    temp_C3K = cbind(PG_Ch_C3K$Wave, PG_Ch_C3K$Zspec[[min(which(PG_Ch_C3K$Z >= Z))]][min(which(PG_Ch_C3K$Age >= 10^logAge)),])
    temp_MILES = cbind(PG_Ch_MILES$Wave, PG_Ch_MILES$Zspec[[min(which(PG_Ch_MILES$Z >= Z))]][min(which(PG_Ch_MILES$Age >= 10^logAge)),])
    temp_BASEL = cbind(PG_Ch_BASEL$Wave, PG_Ch_BASEL$Zspec[[min(which(PG_Ch_BASEL$Z >= Z))]][min(which(PG_Ch_BASEL$Age >= 10^logAge)),])
    
    magplot(NA, NA, log='xy', xlim=c(1e2,1e5), ylim=c(1e-8,3), xlab='Wavelength / Ang', ylab='Lsol / Ang (for 1 Msol)')
    lines(temp_PHOENIX, col='blue')
    lines(temp_C3K, col='darkgreen')
    lines(temp_MILES, col='red')
    lines(temp_BASEL, col='orange')
    lines(BC03lr$Wave, BC03lr$Zspec[[min(which(BC03lr$Z >= Z))]][min(which(BC03lr$Age >= 10^logAge)),], col='black')
    legend('topleft', legend=c('PG_Ch_PHOENIX', 'PG_Ch_C3K', 'PG_Ch_MILES', 'PG_Ch_BASEL', 'BC03'), col=c('blue', 'darkgreen', 'red', 'orange', 'black'), lty=1)
    legend('topright', legend=c(paste('Z:',Z), paste('log(Age/Yr):',logAge)))
  }
}

The four different ‘base’ stellar atmosphere sets included with ProGeny show remarkable consistency for the most part. The main differences are caused by their differing native resolution (where PHOENIX is the highest).

Using modifications to the above general approach, in the ProGeny I paper we also look at the impact of switching isochrone (keeping the stellar spectra and IMF fixed), and of switching IMF (keeping the isochrones and stellar spectra fixed). In general we remark that switching isochrones has the most impact, followed by stellar spectra. IMFs can have a large impact in mass-to-light if switching from Chabrier or Kroupa to Salpeter, but the spectral shape is fairly consistent.

Evolution of Stellar Mass Remaining

It is also interesting to compare the stellar evolutionary tracks, in particular the mass remaining in stars with time, which is relevant for calculating standard stellar masses (which needs mass in stars, not total mass formed). Here we use PG_Ch_PHOENIX, but all four will be the same since we are using the same MIST isochrones in all cases.

for(Z in BC03lr$Z){
  magplot(PG_Ch_PHOENIX$Age, PG_Ch_PHOENIX$Zevo[[min(which(PG_Ch_PHOENIX$Z >= Z))]]$SMstar, type='l', ylim=c(0,1),
          log='x', xlab='Age / Yr', ylab = 'SMfrac', col='blue')
  lines(BC03lr$Age, BC03lr$Zevo[[min(which(BC03lr$Z >= Z))]]$SMstar, col='black')
  legend('bottomleft', legend=c('ProGeny [MIST]', 'BC03'), col=c('blue', 'black'), lty=1)
  legend('topright', legend=paste('Z:',Z))
}

Colour Outputs

In the ProGeny I paper we also look at the impact of colour differences. Below we look at the impact of switching stellar spectral libraries, where using broadband colours make it easier to see small differences because the spectra covers many dex in luminosity.

NUV_filt = getfilt('NUV')
r_filt = getfilt('r_VST')
Ks_filt = getfilt('Ks_VISTA')
    
for(Z in BC03lr$Z){
  mag_NUV_r_C3K = {}
  mag_r_Ks_C3K = {}
  mag_NUV_r_PHOENIX = {}
  mag_r_Ks_PHOENIX = {}
  mag_NUV_r_MILES = {}
  mag_r_Ks_MILES = {}
  mag_NUV_r_BASEL = {}
  mag_r_Ks_BASEL = {}
  mag_NUV_r_BC03 = {}
  mag_r_Ks_BC03 = {}
  
  for(logAge in 6:10){
    temp_C3K = specReBin(PG_Ch_C3K$Wave, PG_Ch_C3K$Zspec[[min(which(PG_Ch_C3K$Z >= Z))]][min(which(PG_Ch_C3K$Age >= 10^logAge)),], BC03lr$Wave)
    temp_PHOENIX = specReBin(PG_Ch_PHOENIX$Wave, PG_Ch_PHOENIX$Zspec[[min(which(PG_Ch_PHOENIX$Z >= Z))]][min(which(PG_Ch_PHOENIX$Age >= 10^logAge)),], BC03lr$Wave)
    temp_MILES = specReBin(PG_Ch_MILES$Wave, PG_Ch_MILES$Zspec[[min(which(PG_Ch_MILES$Z >= Z))]][min(which(PG_Ch_MILES$Age >= 10^logAge)),], BC03lr$Wave)
    temp_BASEL = specReBin(PG_Ch_BASEL$Wave, PG_Ch_BASEL$Zspec[[min(which(PG_Ch_BASEL$Z >= Z))]][min(which(PG_Ch_BASEL$Age >= 10^logAge)),], BC03lr$Wave)
    temp_BC03 = cbind(BC03lr$Wave, BC03lr$Zspec[[min(which(BC03lr$Z >= Z))]][min(which(BC03lr$Age >= 10^logAge)),])
    
    mag_NUV_r_C3K = c(mag_NUV_r_C3K, magABcalc(temp_C3K, filter=NUV_filt) - magABcalc(temp_C3K, filter=r_filt))
    mag_r_Ks_C3K = c(mag_r_Ks_C3K, magABcalc(temp_C3K, filter=r_filt) - magABcalc(temp_C3K, filter=Ks_filt))
    mag_NUV_r_PHOENIX = c(mag_NUV_r_PHOENIX, magABcalc(temp_PHOENIX, filter=NUV_filt) - magABcalc(temp_PHOENIX, filter=r_filt))
    mag_r_Ks_PHOENIX = c(mag_r_Ks_PHOENIX, magABcalc(temp_PHOENIX, filter=r_filt) - magABcalc(temp_PHOENIX, filter=Ks_filt))
    mag_NUV_r_MILES = c(mag_NUV_r_MILES, magABcalc(temp_MILES, filter=NUV_filt) - magABcalc(temp_MILES, filter=r_filt))
    mag_r_Ks_MILES = c(mag_r_Ks_MILES, magABcalc(temp_MILES, filter=r_filt) - magABcalc(temp_MILES, filter=Ks_filt))
    mag_NUV_r_BASEL = c(mag_NUV_r_BASEL, magABcalc(temp_BASEL, filter=NUV_filt) - magABcalc(temp_BASEL, filter=r_filt))
    mag_r_Ks_BASEL = c(mag_r_Ks_BASEL, magABcalc(temp_BASEL, filter=r_filt) - magABcalc(temp_BASEL, filter=Ks_filt))
    mag_NUV_r_BC03 = c(mag_NUV_r_BC03, magABcalc(temp_BC03, filter=NUV_filt) - magABcalc(temp_BC03, filter=r_filt))
    mag_r_Ks_BC03 = c(mag_r_Ks_BC03, magABcalc(temp_BC03, filter=r_filt) - magABcalc(temp_BC03, filter=Ks_filt))
  }
  
  magplot(mag_NUV_r_C3K, mag_r_Ks_C3K, type='l', col='blue', xlim=c(-2,8), ylim=c(-2.5,2), xlab='NUV - r', ylab='r - Ks')
  points(mag_NUV_r_C3K, mag_r_Ks_C3K, cex=seq(0.5,2,len=5), col='blue')
  lines(mag_NUV_r_PHOENIX, mag_r_Ks_PHOENIX, col='darkgreen')
  points(mag_NUV_r_PHOENIX, mag_r_Ks_PHOENIX, cex=seq(0.5,2,len=5), col='darkgreen')
  lines(mag_NUV_r_MILES, mag_r_Ks_MILES, col='red')
  points(mag_NUV_r_MILES, mag_r_Ks_MILES, cex=seq(0.5,2,len=5), col='red')
  lines(mag_NUV_r_BASEL, mag_r_Ks_BASEL, col='orange')
  points(mag_NUV_r_BASEL, mag_r_Ks_BASEL, cex=seq(0.5,2,len=5), col='orange')
  lines(mag_NUV_r_BC03, mag_r_Ks_BC03, col='black')
  points(mag_NUV_r_BC03, mag_r_Ks_BC03, cex=seq(0.5,2,len=5), col='black')
  legend('bottomright', legend=c('PG_Ch_C3K', 'PG_Ch_PHOENIX', 'PG_Ch_MILES', 'PG_Ch_BASEL', 'BC03'), col=c('blue', 'darkgreen', 'red', 'orange', 'black'), lty=1)
  legend('topright', legend=paste('Z:',Z))
  legend('bottom', legend=c(expression(10^6*' [Yr]'), expression(10^7*' [Yr]'), expression(10^8*' [Yr]'), expression(10^9*' [Yr]'), expression(10^10*' [Yr]')), pt.cex=seq(0.5,2,len=5), pch=1)
}