The absolutely Luminosity about stars.

[d11122]

hello,I'm a graduate students to use Phoebe2.3, I have a question in my work.
I have fitted the light curve about ztf data,I‘m use the normalized data.When I analyse the result of mcmc ,I want to compute the absolutely Luminosity of the primary and secondary star (not the passband absolutely Luminosity !). I read all docs about Phoebe 2.3,but I don't find the answer.Dose Phoebe provide the function？
Thank you!
Daimin

[kecnry]

Hi Daimin,

Great question.... not directly, but you can get a decent approximation by having phoebe calculate the absolute passband luminosity in a bolometric passband (which doesn't quite capture the entire spectrum, but should cover most of it). You would need to do something like the following:

b.add_dataset('lc', passband='Bolometric:900-40000', pblum_mode='absolute', dataset='lcbol')
pblums = b.compute_pblums(dataset='lcbol')

which for the default binary gives the following output:

{'pblum@primary@lcbol': <Quantity 4.49281145e+26 W>, 'pblum@secondary@lcbol': <Quantity 4.49281145e+26 W>}

the values of this dictionary are astropy quantities, so can be converted to solar units, if you'd like:

pblums['pblum@primary@lcbol'].to('solLum')

which gives:

<Quantity 1.1736707 solLum>

Hope this helps!
Kyle

[d11122]

thank you very much ,I will do with your responses!

[d11122]

hello ，
dear Kyle!
I have gotten the absolute luminosity which I wanted! but I have one more small question, how to calculate the error? I also try to look up the docs,I got nothing.
Thank you again from my heart!
Daimin

[kecnry]

Unfortunately there is no easy way to do that, although I will note this as something to consider implementing for one of the next releases.

Below are a few options you have to do this manually... but I realize neither of them are simple or pretty. I'll try to adapt much of this into a tutorial on the website so its easier for anyone in the future to find until something more convenient is implemented.

Your two options are to:

1. Use Stefan-Boltzmann approximation which assumes blackbodies and spherical stars, but is algebraic so allows to use constraints to easily and quickly propagate distributions.
2. Use the bolometric passband (which doesn't cover the whole spectrum, but a pretty wide range) and compute the luminosity from the meshes and atmosphere models themselves for a large number of samples in requiv and teff manually.

In either case, you'll need to have uncertainties or posteriors (from MCMC for example) for at least the radii and temperatures of your stars.

1: Stefan Boltzmann approximation via constraints

You can use PHOEBE to propagate these distributions to other constrained parameters, which means that we can do the same for luminosity by creating a custom parameter/constraint.

Unfortunately these must be algebraic, so cannot make calls to compute_pblums (which requires generating the mesh under-the-hood). However, if your stars are sufficiently detached that Stefan-Boltzmann provides a good approximation, then this is your easiest option. You could do something like the following:

import phoebe
from phoebe import u,c
import numpy as np

b = phoebe.default_binary()

lum_sb_def = phoebe.parameters.FloatParameter(qualifier='lum_sb',
                                              default_unit=u.solLum,
                                              value=1, 
                                              description='Stefan-Boltzmann approximation of luminosity',
                                              latexfmt=r'L_\mathrm{{ sb, {component} }}')
                                                                                 
for component in b.hierarchy.get_stars():
    lum_sb_param, created = b.get_or_create('lum_sb', lum_sb_def.copy(), context='component', component=component)
    
    lhs = b.get_parameter(qualifier='lum_sb', component=component)
    R = b.get_parameter(qualifier='requiv', component=component)
    T = b.get_parameter(qualifier='teff', component=component)
    rhs = 4*np.pi*R**2*c.sigma_sb.to('solLum / (solTeff4 solRad2)').value*T**4
    
    b.add_constraint(lhs, rhs)

print(b.filter(qualifier='lum_sb', context='component'))

which gives

ParameterSet: 2 parameters
C        lum_sb@primary@component: 1.1676082413206266 solLum
C      lum_sb@secondary@component: 1.1676082413206266 solLum

which we can then compare to the "bolometric" value computed from the meshes

b.add_dataset('lc', passband='Bolometric:900-40000', pblum_mode='absolute', dataset='lcbol')
print({k: v.to('solLum') for k,v in b.compute_pblums(dataset='lcbol').items()})

which gave

{'pblum@primary@lcbol': <Quantity 1.1736707 solLum>, 
'pblum@secondary@lcbol': <Quantity 1.1736707 solLum>}

So in this case, these are not identical, but are pretty close. Its up to you which approximation you think is most appropriate for your system.

For the sake of this example, I'll just create mock gaussian distributions on requiv and teff of the primary star. In reality, you'll either need to use posteriors from MCMC in their place or set the distributions based on known uncertainties.

b.add_distribution({'requiv@primary': phoebe.gaussian_around(0.1),
                    'teff@primary': phoebe.gaussian_around(100)},
                   distribution='mocked_posteriors')

Now you can easily propagate these distributions through to the luminosity parameter you created (and can also call uncertainties_from_distribution_collection to access the uncertainties directly):

_ = b.plot_distribution_collection(parameters=['requiv@primary', 'teff@primary', 'lum_sb@primary'], show=True)

2: "Bolometric" Luminosity using atmospheres by manually sampling and calling compute_pblums

Alternatively, if the Stefan-Boltzmann approximation is insufficient, you can manually sample from your posteriors, set those values, call compute_pblums, store those sampled computed luminosities, and create and plot the distribution from those. Note I am using 1000 samples below and it still takes a few minutes to compute, but you may want to use even more samples to get cleaner distributions out. This is unfortunately quite expensive as the system needs to be re-initialized and re-meshed N times.

N = 1000
samples = b.sample_distribution_collection(sample_size=N)
requivs = samples['requiv@primary@star@component']
teffs = samples['teff@primary@star@component']

lums = np.empty(N)

# loop over each of the sampled values of requiv and teffs from the (posterior) distributions
for i, (requiv, teff) in enumerate(zip(requivs, teffs)):
    b.set_value(qualifier='requiv', component='primary', context='component', value=requiv)
    b.set_value(qualifier='teff', component='primary', context='component', value=teff)
    
    # compute and store the bolometric passband luminosity
    pblums = b.compute_pblums(dataset='lcbol')   
    lums[i] = pblums['pblum@primary@lcbol'].to('solLum').value

# create a multivariate samples distribution object
mvs = phoebe.dependencies.distl.mvsamples(np.array([requivs, teffs, lums]).T,
                                          labels=['requiv', 'teff', 'lum'],
                                          labels_latex=[r'R_\mathrm{{ equiv, primary }}',
                                                        r'T_\mathrm{{ eff, primary }}',
                                                        r'L_\mathrm{{ bol, primary }}'])

_ = mvs.plot(show=True)

You'll notice this isn't the exact same result as above. Increasing the number of samples should help get rid of random noise and bring the sample distributions of requiv and teff closer to their defined gaussian distribution, but there will always be a slight difference between the two different "methods" of estimating the luminosity.

[d11122]

Thank you for your careful and patient answers！