Forward modelling observed temperatures

[jsimpson-astro]

Hi PHOEBE team,

I am running a relatively simple model grid on some irradiated MS + NS binaries with PHOEBE, with the goal of simulating the temperature that would be observed (i.e. in optical spectra) at different orbital phases and in phase bins of varying width.
Currently I have something working that gives reasonable temperatures, but I wanted to check the logic is correct.

What I am doing right now:

• Computing meshes over a number of phases within each phase bin: I use 100 phase points to approximate an observation in a phase bin, around a given central phase with a phase bin width of 0.25. So e.g. at phase 0.5, I would compute 100 meshes between phase 0.375 and phase 0.625.
• At each point, for each mesh, calculating the flux-weighted temperature over the MS star: I take the teffs, visibilities, areas, and mus/nws (as I understand, they are the same) and use the visibilities to mask each of them where visibilities != 1 (eclipses are off). Then I do np.average(teffs**4, weights=mus*areas)**0.25 or ((mus * areas * teffs**4).sum() / (mus * areas).sum())**0.25.
• Then, I average the temperature calculated at each point together, again flux-weighted, throughout the phase bin by doing (obs_teffs**4).mean()**0.25.

So then I have (I hope) temperatures at a given phase bin that should approximate what would be observed from a spectroscopic measurement over the same phases.
Is there anything crucial I've missed here within PHOEBE? I had a look at how flux-weighted RVs are done, with abs_intensities weighting, but I do not have this as I am only working with the mesh dataset.

Thank you for your help, and your hard work as always!
Jordan

[Raindogjones]

Hi Jordan,
On the face of it, that looks reasonable to me but with the pretty big caveat that irradiation in real life is not as simple as absorbing a fraction (the albedo, irrad_frac_refl_bol in phoebe) of the incoming flux and that increasing the temperature of the surface element upon which it is incident. Put simply, an irradiated atmosphere is not the same as the atmosphere of the same star with an increased effective temperature, and this alters the emergent spectrum. See, for example, the discussions of this in Wilson (1990) and Horvat et al. (2019).
Cheers,
Dave

[jsimpson-astro]

Hi Dave,
Thanks for the quick response and the references. I'm treating this grid as more of a toy system (with irrad_frac_refl_bol = 1) just to explore the relative effects with observing different phases, temperatures, inclinations etc..
I'll keep the caveats you mentioned in mind though, and will definitely have to look into it further for future modelling.
Thanks for your help!
Jordan

[kecnry]

This is a cool idea! We are hoping to eventually have a much easier way to access this sort of information (see #407 (comment)), so please keep us updated how it goes or if you need access to anything else.

If you want to account for any type of irradiation or passband effects by weighing by intensities instead of temperatures, you should be able to create a light curve dataset (with a bolometric passband if you want irradiation but not an actual passband), with no times, and then request the abs_intensities columns to be populated. I'd be curious to see the difference it makes in the result if you decide to try it!

[Raindogjones]

Nice idea, Kyle. Using a bolometric light curve would be a good, maybe simpler, proxy for the temperature as it is just proportional to T^4, so no need to worry about weighting temperatures from the mesh.

[jsimpson-astro]

Hi Kyle, thanks for the suggestion - I gave it a try with a couple of models, both heavily irradiated, to compare the two calculations. I changed the calculation for summing the temperature of each mesh to np.average(teffs, weights=abs_intensities*areas*mus) instead of np.average(teffs**4, weights=areas*mus)**0.25. The averaging of each point is still the same however, the mean of T^4 then to the 1/4. The results are below - you can see it biases the 'observations' more to higher temperatures which does look a bit more in line with observations.

Could I ask what the abs_intensities mean vs. the intensities? And also, if I wanted to weight the averaging of each point by the luminosity of the secondary at each point, should I just weight that average by the sum of abs_intensities*areas*mus?
Thanks!
Jordan

[kecnry]

Great, thanks for the update!

Could I ask what the abs_intensities mean vs. the intensities?

intensities is just abs_intensities multiplied by a single scale factor such that the star's luminosity matches the prescribed passband luminosity (assuming pblum_mode is left at the default of 'component-coupled'). Since you're using it as a weight, I don't expect any difference between the two.

And also, if I wanted to weight the averaging of each point by the luminosity of the secondary at each point, should I just weight that average by the sum of abs_intensities*areas*mus?

Technically the luminosity would be np.sum(abs_normal_intensities*areas*ldint)*ptfarea*np.pi (or again, normal_intensities to get the scaled luminosity), but dropping the constant scalars is fine for the weight. You definitely do want the normal instead of projected intensities though. The mesh does also optionally expose pblum_ext or abs_pblum_ext (the ext meaning extrinsic - since these will be observed luminosities as a function of time including spots, phase-dependent distortion, irradiation, etc, whereas the input pblum parameter is intrinsic with these effects excluded), which you could then use directly.