Day/night side temperature of a planet from phase curve

[LucaNap]

Hey there, I've just started using Starry and I want to say it's working great.

I am wondering if there's a way to directly estimate the day (or night) side temperature of a planet in a phase curve.. as far as I understand the planet luminosity (Fp/fs) is the average during the orbit, is it not? Then basically all we would need are the Fp/fs estimates at phase 0.0 (mid transit) and 0.5 (mid eclipse) as if the planet didn't transit/eclipse, and then evaluate the brightness temperature.. but perhaps there's already a way to do this? Thank you

[rodluger]

Hey there, I've just started using Starry and I want to say it's working great.

Cool!

I am wondering if there's a way to directly estimate the day (or night) side temperature of a planet in a phase curve.. as far as I understand the Fpfs is the average during the orbit, is it not?

I'm not sure what you mean by FpFs being the average during the orbit. The quantity you measure from a phase curve or secondary eclipse is the flux of the star Fs plus the flux of the planet Fp in some arbitrary units, which you can usually use to figure out the ratio Fp / Fs at any given point along the orbit. So it's not an average over the orbit.

Then basically all we would need are the Fpfs estimates at phase 0.0 (mid transit) and 0.5 (mid eclipse) as if the planet didn't transit/eclipse (and evaluate the brightness temperature).. but perhaps there's already a way to do this? Thank you

With regards to the temperature, FpFs is the planet emission integrated over the observer-facing hemisphere, so there's no unique solution to things like the "dayside temperature" of the planet: it will depend on how the albedo, emissivity, and heat transport vary over the surface. In the instant re-radiation limit, you could assume that the emission follows the illumination profile, so the intensity map is a dipole centered on the sub-stellar point. This might be a decent approximation in the case of, e.g., a tidally-locked, rocky exoplanet. In this limit, you could model the planet with a ydeg=1 surface map and solve for the coefficients using map.solve() or MCMC.

Whatever model structure you assume, the spherical harmonic coefficients you end up with describe the intensity profile of the planet (which you can access by calling, e.g., map.intensity()). You then have to use Planck's law to convert these intensities to a temperature. Perhaps @ericagol can recommend good references on how this is done in practice?

[LucaNap]

I'm not sure what you mean by FpFs being the average during the orbit. The quantity you measure from a phase curve or secondary eclipse is the flux of the star Fs plus the flux of the planet Fp in some arbitrary units, which you can usually use to figure out the ratio Fp / Fs at any given point along the orbit. So it's not an average over the orbit.

Let's say I have a phase curve of a star-planet system (total flux over orbital phase), as far as I understand in order to derive Fp/fs all I have to do is normalize the flux so that total flux = 1.0 during the secondary eclipse, and then subtract 1. At this point the estimated "flux_model" is Fp/fs over the orbital phase... but how is it related to the planet luminosity defined for example as amp (10*log_amp) in the HotJupiterPhaseCurve notebook? When I run the netbook code with planet[1,0]=0.5, the Fp/fs value is close to zero in the night side (transit) and close to 2amp in the day side (secondary eclipse), while the average of Fp/fs outside of the transit and secondary eclipse is the exact value of the defined planet amp (1 ppt), that's what I was trying to say in the opening post which was poorly written... but most importantly:

With regards to the temperature, FpFs is the planet emission integrated over the observer-facing hemisphere, so there's no unique solution to things like the "dayside temperature" of the planet: it will depend on how the albedo, emissivity, and heat transport vary over the surface. In the instant re-radiation limit, you could assume that the emission follows the illumination profile, so the intensity map is a dipole centered on the sub-stellar point. This might be a decent approximation in the case of, e.g., a tidally-locked, rocky exoplanet. In this limit, you could model the planet with a ydeg=1 surface map and solve for the coefficients using map.solve() or MCMC.

Yes! That's what I am trying do to. I have set ydeg=1 and I am solving for log_amp, offset and planet[1,0].. and I find reasonable parameters (+ mcmc errors) with a very good looking fit. I get that the real temperature depends on lots of things, but a "brightness" temperature can easily be estimated with Planck's law and star/planet radius', as you have said (but any suggestion is very appreciated!). The issue is that in order to estimate such temperatures I need Fp/fs model values during both eclipses, and right now the only way to do it is by taking the average of Fpfs just before ingress and after egress, as of course the actual Fp/fs is zero during secondary eclipse (because the planet can't be seen) and less than zero during the first eclipse (because part of the star is obscured). Basically I'd need the expected Fp/fs at phase 0.0 and 0.5 as if there was no eclipse, like in the following figure... or am I doing something wrong?

Sorry for the (probably) dumb questions but I was not able to find these informations on the docs. Thank you very much!

[rodluger]

I think I understand. No worries, these aren't dumb questions!

You can certainly estimate the temperature from Fp/Fs, but there's a bit of an inconsistency in what you're doing. This might be fine -- it's a very simple model, and is typically what people do for exoplanets -- but it's worth thinking about it. Because you're modeling the planet as a dipole, there's no single dayside temperature: the temperature peaks at the center of the hotspot and falls off monotonically toward the nightside. So it doesn't really make sense to try to estimate a single temperature from the hemispherically-averaged flux Fp/Fs. This might give you an answer in the right ballpark, but I'd recommend a different approach. Once you've fit for the spherical harmonic coefficients (and updated their values in the planet map), you can call

I = map.intensity(lat=lat, lon=lon)

to get the intensity at a certain latitude and longitude on the surface of the planet. This is (up to a constant) equal to the spectral radiance B that goes into Planck's law, so if I define a function

EDIT: This is wrong; see @LucaNap's comment below!

def T(B):
   return (2 * h * c ** 2 / lam ** 5) * 1 / (np.exp(h * c / (lam * k * T)) - 1)

(with appropriate values for the constants h, c, k, and the wavelength lam), I can simply call T(I) to get the temperature anywhere on the surface under my simple dipole model. Do it at latitude=longitude=0 to get the temperature at the center of the hotspot.

If you want to compute the intensity over the entire portion of the planet facing the observer, it's easier to do

I = map.render(res=300, theta=theta)

which will return the value of the intensities on a square 300x300 grid enclosing the unit disk at rotational phase theta. (it's what you see when you call map.show()). You can then call your temperature function on this array to get the temperature everywhere on the dayside. You can then take integrals, averages, etc.

The main thing to be careful about is units. Everything in starry is unitless, since it's defined relative to the star. Just remember that in starry, the quantity Fp/Fs is the integral of the intensity over the projected disk of the planet. So if you compute I using map.render as above, you should find that

Fp/Fs = np.nanmean(I) * np.pi

So figure out the conversion factor for Fp/Fs that gets you into the proper units for flux, and use that to scale the intensity when computing the temperature.

I'm sorry this is all so complicated. I developed starry as a tool for mapping stellar/exoplanet surfaces, so getting temperatures out isn't all that easy. Let me know how this goes. I'll try to put together a tutorial on this soon!

[LucaNap]

Ok, I understand it way better now! By the way, in order to get the temperature from the intensity anywhere on the surface I think you want a function like

def T(I):
   return ( h * c / ( lam *k_B ) / np.log(1 + 2 * h * c**2 / ( I * lam**5)) )

unless I am mistaken, since yours was returning the spectral radiance B(T) over a temperature T. So, for the day/night side intensity I should do something like

I = map.render(res=300, theta=offset + 180.0) # night side or
I = map.render(res=300, theta=offset) # day side

At this point I think I have to set the proper units and scale by multiplying I to the intensity of the star itself, the ratio of the surfaces (star over planet) and pi

def Intensity(T):
   return (2 * h * c ** 2 / lam ** 5) * 1 / (np.exp(h * c / (lam * k_B * T_star)) - 1)

Fp_fs = I * Intensity(T_star) * np.pi * (r_s / r_p)**2

Now np.nanmean(T(Fp_fs)) should be an approximated temperature. I have tried this code and I get the temperature I roughly expected to for the day side. Thank you for all your help! Looking forward to the tutorial.

[rodluger]

Yes, that all makes sense to me! Glad to hear it's working!

[LucaNap]

I still have a small doubt remaining: shouldn't the overall temperature be higher at theta=offset rather than theta=0.0 ?

[rodluger]

Good question.

The way it is modeled in the hot Jupiter example, the offset is telling starry.System how to orient the planet during secondary eclipse. This is entirely separate from the planet map stuff, which is all defined on a fixed latitude-longitude grid in a static frame. So things like b.map.intensity, b.map.render, etc., don't know about offset -- everything is defined relative to the center of the coordinate system. Specifically, this line

b.map[1, 0] = 0.5

orients the dipole such that its maximum intensity is at lon=0 and minimum intensity is at lon=180. When you call render(), the theta parameter is the angular phase of the map relative to the phase at which lon=0 is at the sub-observer point (so theta -- in this context -- is really just the longitude of the sub-observer point).

I'm sorry this was confusing! Hope it makes more sense now. I'll be sure to explain all this in the tutorial -- it might just take me a while to write it up.

[LucaNap]

Thank you @rodluger. I am wondering if you have any advices on how to estimate the uncertainty on these temperatures, since I don't think you can use np.nanmean() in the pymc3 model (but only outside pymc3 with starry.config.lazy = False). Is there a workaround?

Hopefully this conversation will be useful to other people!

[rodluger]

There's probably a way to do this that avoids the numerical integration (i.e., pixelizing the surface and computing the mean), but in the time being you could do

idx = np.where(~np.isnan(starry.Map(ydeg=1).render(res=300).eval()))

which will give you the indices of all the non-NaN values in the projected disk. You only need to do this once at the beginning -- these indices won't change. Then you can just do

tt.mean(Fp_fs[idx])

within the pymc3 model and it should work, where you should make sure to import

import theano.tensor as tt

The reason for the extra step here is that theano doesn't have a nanmean function, so we need to explicitly tell it not to average over elements that are NaN.

[rodluger]

And yes -- please keep posting issues / questions, as this is helpful to others (and to me!)