Evaluating the source and nature of ionizing photons.

 Objective: we want to understand which populations produce large counts of ionizing photons and which stars dominate that production

This code is in IDL (tested on IDL v8.5.1 for Mac). "readcol" is from the IDL astronomy library. You will need to have downloaded or copied the files from imf135_300 of the BPASS v2.2.1 release.

You can download this tutorial as an IDL procedure here.

;; ********************************************************************************
;; First setup BPASS variables - modify "root" to match your path to models.

ext='imf135_300'
root='BPASSv2.2.1_bin-imf135_300/' 

;; Metallicities

Zlab=['zem5','zem4','z001','z002','z003','z004','z006','z008','z010','z014','z020','z030','z040']
Z=[0.0001,0.0001,0.001,0.002,0.003,0.004,0.006,0.008,0.010,0.014,0.020,0.030,0.040]
numz=n_elements(zlab)

;; Age and time steps (the width of each time bin varies due to
;; logarithmic spacing)

n_age=51
logage=findgen(n_age)*0.1 + 6.0
dage=fltarr(n_age)
dage[0]=10.^6.05
for ll=1,n_age-1 do dage[ll]=10^(logage[ll]+0.05)-10^(logage[ll]-0.05)

;; Some convenient constants (which we may or may not need)

lsun=3.846d26 ; W, J/s
c=3.d8 ; m/s
h=6.62606957d-34 ; m2 kg / s
a= 1d-10 ; m/AA

;; ********************************************************************************
;; Part I: Plotting the ionizing photon production rates

;; We want to know about ionizing photons so first check the manual

;; We have two choices:
;; (i) read in the full SED files (51 columns x 10,000 rows per metallicity)
;; (ii) use the precalculated ionizing photon flux files (51 cols x 5 rows per Z)

;; For option i (which we're not going to use)
;;
;; model=dblarr(n_age,30000,numz)
;; for zn=0,numz-1 do begin
;; openr, 1, root+'spectra-bin-'+ext+'.'+zlab[zn]+'.dat'
;; readf, 2, lama, a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,$
;; a18,a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,a35,a36,$
;; a37,a38,a39,a40,a41, a42,a43,a44,a45,a46,a47,a48,a49,a50,a51, format='(52(e16))'
;; lam[i]=lama
;; model[0:40,i]=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,$
;; a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,a35,a36,$
;; a37,a38,a39,a40,a41]
;; close, 1
;; endfor
;;
;; sub2=where(lam lt 912. and lam gt 10.)
;; uvphot=dblarr(numz,n_age)
;; For zn=0,numz-1 do for ll=0,n_age-1 do uvphot[zn,ll]=alog10(total(model(ll,sub2)*lsun/(h*c/(a*lam[sub2]))))

;; We're going to use option ii for now. [Note: if copy/pasting from command line you'll have to preceed for loops with ".run" and follow them with "end".]

uvphot=dblarr(numz,n_age)
for zn=0,numz-1 do begin
 readcol, root+'ionizing-bin-'+ext+'.'+zlab[zn]+'.dat', lage, nion
 uvphot[zn,*]=nion
endfor

;; Okay, let's take a first look at these

plot, logage, uvphot[0,*],xtit='Log (Age/yrs)', ytit='log (N!Dion!N per s)',xr=[6,10],yr=[46,54],xs=1,ys=1
for zn=1,numz-1 do oplot, logage, uvphot[zn,*],line=zn mod 6

;; First thing to note: the ionizing photon flux drops with metallicity so is dominated by the lowest metallicity you have

;; Second thing to note: The two extreme metallicities drop off faster than the others (which are similar).

;; CAUTION: This assumes that the entire population is the same age and has an initial total mass of 10^6 Msun.

;; Let's think about a slightly different star formation history: Constant star formation at a rate of 1 Msun/yr

;; Consider the stars which formed a time log(years)=6.0+i*0.1 before present. They will have logage[i] but how many of them are there? Well, that will depend on the width of the age bin, dage[i] since 1 Msun of stars forms per year of bin width. So the wide bins at old ages still contribute a large amount to the number of stars and photon flux, even though each star is emitting relatively little UV.

plot, logage, uvphot[0,*]+alog10(dage),xtit='Log (Age/yrs)', ytit='total log (N!Dion!N/s) from bin',xr=[6,10],yr=[55,60],xs=1,ys=1
for zn=1,numz-1 do oplot, logage, uvphot[zn,*]+alog10(dage),line=zn mod 6

;; Okay, but that's not quite what we want. We're interested in continuous star formation which started a time logage[i] ago, not just stars which formed logage[i] before present, so we will see the sum of all ages *up to* logage[i]
;; i.e.

cuvphot=dblarr(numz,n_age)
for zn=0,numz-1 do for ll=0,n_age-1 do for kk=0,ll do cuvphot[zn,ll]=cuvphot[zn,ll]+dage[kk]*(1d1^uvphot[zn,kk])
cuvphot=alog10(cuvphot)

;; Now let's take a look at this.

plot, logage, cuvphot[0,*],xtit='Log (Age/yrs)', ytit='total log (N!Dion!N/s) from bin',xr=[6,10],yr=[58.5,60],xs=1,ys=1
for zn=1,numz-1 do oplot, logage, cuvphot[zn,*],line=zn mod 6

It should look like this.

;; Things to note:
;; - There is still a clear trend with metallicity
;; - Changing the star formation history has had a far larger effect than changing the metallicity for single bursts
;; - The ionizing photon flux stabilises somewhere around logage=8.

;; So what stars are driving this evolution? The stars which dominate at logage=6 and the obvious massive OB stars. But what's still contributing at logage=7-8? To find out, we need to look at another output files.

;; ********************************************************************************
;; Part II: Figuring out the stellar populations

;; We want to know about the contributions of different types of stars. A good place to look for this is in the "numbers" files

;; Check the manual: These have 51 rows and 21 columns - an age column and numbers of ten different stellar types divided into two luminosity bins.

types=['O','Of','B','A','YSG','K','M','WNH', 'WN','WC']
ntype=n_elements(types)
starcount=fltarr(numz, n_age,n_elements(types),2)

for zn=0,numz-1 do begin
  openr, 1, root+'numbers-bin-'+ext+'.'+zlab[zn]+'.dat'
  temp=dblarr(21)
  for ll=0,n_age-1 do begin
    readf,1, temp
    starcount[zn,ll,*,0]=temp[1:10]
    starcount[zn,ll,*,1]=temp[11:20]
  endfor 
  close, 1
endfor

;; okay, let's take a look at the very luminous stars at Z=0.1 Zsun metallicity:

zn=3 ;; z002
plot, logage, starcount[zn,*,0,0],xr=[6,10],yr=[1,2000],xtit='Log (Age/yrs)', ytit='Number of Stars',/ylog
for i=1,ntype-1 do oplot, logage, starcount[zn,*,i,0],line=i mod 6

;; This is fine but hard to tell which line is which, so let's look at particular things

;; Let's put a key on:

for i=0,ntype-1 do plots, [0.85,0.9],0.9-0.05*i*[1,1],/norm,line=i mod 6
for i=0,ntype-1 do xyouts, 0.9,0.9-0.05*i,types[i],/norm

;; Still some confusion so where are the Wolf-Rayets? Let's mark them in colour.

for i=0,ntype-1 do if strmatch(types[i],'*W*') eq 1 then oplot, logage, starcount[zn,*,i,0],line=i mod 6,col=250
for i=0,ntype-1 do if strmatch(types[i],'*W*') eq 1 then plots, [0.85,0.9],0.9-0.05*i*[1,1],/norm,line=i mod 6,col=250

Your plot should look like this.

;; Hmm, so there's a tail of these extending to logage=7.3, but most are short lived peaking at logage=6.4. The dotted line might be giving us a late kick between logage 7 and 8. Let's check which that is.

;; Right so these are the WNHs - Stripped envelope massive stars, which still have some hydrogen in their atmospheres.

;; But we've only looked at the brightest stars. What about the less luminous but more numerous ones? How numerous?

print, max(starcount[zn,*,*,1])

;; We're going to have to adjust the y axis to accomodate them.

plot, logage, starcount[zn,*,0,1],xr=[6,10],yr=[1,1e6],xtit='Log (Age/yrs)', ytit='Number of Stars',/ylog
for i=1,ntype-1 do oplot, logage, starcount[zn,*,i,1],line=i mod 6

;; colour the stripped stars

for i=0,ntype-1 do if strmatch(types[i],'*W*') eq 1 then oplot, logage, starcount[zn,*,i,1],line=i mod 6,col=250

;; and the key (moved down to the lower left)

for i=0,ntype-1 do plots, [0.1,0.15],0.55-0.045*i*[1,1],/norm,line=i mod 6
for i=0,ntype-1 do xyouts, 0.15,0.55-0.045*i,types[i],/norm
for i=0,ntype-1 do if strmatch(types[i],'*W*') eq 1 then plots, [0.1,0.15],0.55-0.045*i*[1,1],/norm,line=i mod 6,col=250

Your plot should look like this.

;; Okay, now we've found something interesting.

;; Most of the low mass stars (BAFGKM) are pretty much unchanged over a Hubble time. But there's a population of dwarf O stars and stripped envelope stars which kicks in about a logage of 7-8 and are likely to be hot and numerous, even if not individually luminous. These are binary products

;; By sheer numbers, the dwarf WNH stars likely dominate the ionizing flux. The dwarf Os cant be dominating because they kick in too late to explain the evolution we saw earlier.

;; So let's focus on these stripped stars and now just plot these.

plot, logage, starcount[zn,*,0,0],xr=[6,10],yr=[1,2000],xtit='Log (Age/yrs)', ytit='Number of Stars',/ylog,/nodata

;; L>4.9

for i=0,ntype-1 do if strmatch(types[i],'*W*') eq 1 then oplot, logage, starcount[zn,*,i,0],line=i mod 6,col=250

;; L<4.9

for i=0,ntype-1 do if strmatch(types[i],'*W*') eq 1 then oplot, logage, starcount[zn,*,i,1],line=i mod 6

Your plot should look like this.

;; So the low luminosity WNHs might well be the dominant population, even above the high luminosity ones, given both their number and the bin-width effect which will boost them in a continuous SF scenario... as long as they're not too faint.

;; What about their spectra? Well, in BPASS these will be assigned stripped-envelope PoWR spectra, despite being lower luminosity, so will be far-UV bright. See also Gotberg et al (2018A&A...615A..78G) for spectra we don't include yet.

;; ********************************************************************************
;; Part III: Plotting a HR diagram.

;; Getting luminosities for individual stars out of BPASS, together with their weighting, is difficult, but there is one way of checking the population properties: the HR diagram.

;; See manual on HR diagram files

;; Note, reading in the HR diagrams is very slow and the format is painful. We're only going to read them at one metallicity

zn=3                            ; z=0.002

file=root+'hrs-bin-'+ext+'.'+zlab[zn]+'.dat'

hr=fltarr(100,100,51,9)
openr,5,file
dummy=fltarr(100)

for iX=0,8 do begin
   for iT=0,50 do begin
      for i=0,99 do begin
         readf,5,dummy
         hr(i,0:99,iT,iX)=dummy
      endfor
   endfor
endfor
close,5

;; We now have 51x9 grids, each with 100x100 bins
;; the 51 are the age bins
;; the 9 breaks down into:
;; (i) 3 grids of log(T) vs log (L)
;; (ii) 3 grids of log(T) vs log(G)
;;(iii) 3 grids of log(T) vs log(TG)

;; We only want the first three. These contain info on hydrogen rich stars, partly-stripped stars and hydrogen-free stars respectively.

;; We've identified WNH stars as of interest. These are partly-stripped so we want

iX=1

;; We are most interested in ages around 10^8 years so

iT=21

logT=findgen(100)*0.1+0.1
logL=findgen(100)*0.1-2.9

;So let's plot this (reversing the temperature axis by convention)

set_plot, 'ps'
device, filename='tutorial_hr.ps',xs=20,ys=20,xoff=0,yoff=0,/times,/encap,/col

plot, logT, logL,xs=1,ys=1,position=[0.1,0.1,0.95,0.95],xr=[max(logT),min(logT)],yr=[min(logL),max(logL)],xtit='Log (T)', ytit='Log (L)'
tv, 255-bytscl(congrid(reverse(alog10(hr[*,*,iT,iX])),500,500),min=-5.,max=9),0.1,0.1,xsize=0.85,ysize=0.85,/norm

device, /close
set_plot, 'x'
spawn, 'pstopdf tutorial_hr.ps; open tutorial_hr.pdf'

You should have this plot.

;; So that's okay, but frankly fairly unenlightning. The long streak heading down and to the right is the White Dwarf cooling sequence which is just starting to appear at this age. Let's zoom in a bit and play around and look at all three compositions in a 3 colour diagram

lhr=dblarr(100,100)
collhr=dblarr(100,100,3)

lhr[*,*]=-99
iX=0 & iT=21
for i=0,99 do for j=0,99 do if hr[i,j,iT,iX] ne 0 then lhr[100-i,j]=alog10(hr[i,j,iT,iX])
collhr[*,*,0]=lhr ;; Red = Hydrogen

lhr[*,*]=-99
iX=1 & iT=21
for i=0,99 do for j=0,99 do if hr[i,j,iT,iX] ne 0 then lhr[100-i,j]=alog10(hr[i,j,iT,iX])
collhr[*,*,1]=lhr ;; Green = partly stripped

lhr[*,*]=-99
iX=2 & iT=21
for i=0,99 do for j=0,99 do if hr[i,j,iT,iX] ne 0 then lhr[100-i,j]=alog10(hr[i,j,iT,iX])
collhr[*,*,2]=lhr ;; Blue = stripped

tmin=50 & tmax=70 & lmin=20 & lmax=80

set_plot, 'ps'
device, filename='tutorial_hr2.ps',xs=22,ys=22,xoff=0,yoff=0,/times,/col,encap=1
!p.font=0 & !p.charsize=2 & !p.thick=3
 sun=sunsymbol()
 plot, logT, logL,xs=1,ys=1,position=[0.15,0.15,0.97,0.99],xr=[logt[tmax],logt[tmin]],yr=[logl[lmin],logl[lmax]],xtit='log T/K', ytit='Log L/L'+sun
 tv, bytscl((collhr[tmin:tmax,lmin:lmax,*]),min=-5,max=0),3.3,3.3,xsize=0.82*22,ysize=0.84*22,/centimeters,true=3
 plot, logT, logL,xs=1,ys=1,position=[0.14,0.15,0.97,0.99],xr=[logt[tmax],logt[tmin]],yr=[logl[lmin],logl[lmax]],/noerase,col=255,xtickn=[' ',' ',' ',' ',' ',' ',' '],ytickn=[' ',' ',' ',' ',' ',' ',' '],/nodata

 xyouts, 0.25, 0.25, 'log(age)='+strtrim(logage(iT),2),/norm,col=255

device, /close
spawn, 'pstopdf tutorial_hr2.ps; open tutorial_hr2.pdf'
set_plot, 'x'

You should have this plot.

;; Note that the scalings here are logarithmic and so the colour scale overemphasises the low number counts relative to the high counts, but they give you the idea! The narrow hydrogen main sequence is clear, together with the large number of stars partly stripped by neighbours and hence burning hotter than might otherwise be expected at their luminosity. Towards the bottom left, the white dwarf sequence shows a range of compositions.

;; Exercise: contrast this with the same plot for logage=6.4 and 9.5

End