unexpected behavior when irradiation is off

[gecheline]

The flux out of eclipse shows unexpected variability when irradiation is off. It's especially prominent in eccentric systems and has caused me issues with convergence to physical values when fitting. I'm attaching a figure of a simple system with spherical stars and circular orbit that showcases the issue.

[pmaxted]

Insufficient information provided to recreate the plot shown.

[gecheline]

Here's the script I used for the particular plot I attached, but as I said in the issue, I've seen this type of weird variability for any system, it's just a matter of whether heat_1 and heat_2 are None/0 or at a certain value > 0.

`def lc(times, r1, r2, incl, vsini, sbratio, period, f_c, f_s, sma=1, shape='roche', heat=None):

return ellc.lc(times,t_zero=0,
            radius_1=r1, radius_2=r2,
            incl=incl,
             vsini_1 = vsini,
            sbratio=sbratio, shape_1='roche', shape_2='roche',
            period=period, f_c=f_c, f_s=f_s,
               heat_1 = heat, heat_2 = heat,
             # Limb darkening
            ld_1='lin', ld_2='lin', ldc_1=0.2, ldc_2=0.2,
            a=sma)

r1 = 0.1 r2 = 0.05 vsini = 10. sbratio = 1. period = 100 f_c = 0 f_s = 0 incl = 85. times = np.linspace(0,200,5000)

noheat = lc(times, r1, r2, incl, vsini, sbratio, period, f_c, f_s, sma=100, shape='sphere', heat=None) wheat = lc(times, r1, r2, incl, vsini, sbratio, period, f_c, f_s, sma=100, shape='sphere', heat=0.6)

plt.plot(times, noheat, label='heat_none') plt.plot(times, wheat, label='heat_06') plt.xlabel('Time (d)') plt.ylabel('Flux') plt.legend() `

[pmaxted]

This is not a script that calculates light curves for spherical stars - the variation between the eclipses is the ellipsoidal effect. sbratio=sbratio, shape_1='roche', shape_2='roche', should be sbratio=sbratio, shape_1=shape, shape_2=shape,

[gecheline]

Of course /facepalm/, you're right, I don't pass the shape on to the function so my example really shows roche stars. It still looks off because the flux is higher without irradiation than with irradiation, but I looked into it further and convinced myself it's just the scaling that's causing that.

So let me provide an actual example that showcases my issue. I uploaded two notebooks that reproduce the system I had trouble with, the exact way I modeled it using the ellc backend in Phoebe (I've tried fitting it with pure ellc too and have the same issues with convergence to reasonable values). I also compute the models using the phoebe2 backend for comparison and they don't agree at all both for spherical stars and roche stars with irradiation off (the flux scaling differences are expected because of the pblum-scaling in Phoebe, but it's the differences in the LC shapes that concern me). We've repeatedly tested the parameter conversion from Phoebe to ellc to make sure everything checks out so I don't believe that is the issue here.

Here's the notebook for the system modeled with roche stars: adra roche and spherical stars: adra spheres . I'm attaching the comparison plots too for reference.

I have found a solution for this particular case by implementing workarounds in Phoebe, but I thought it might point to something in ellc worth looking into.

[gecheline]

I spent more time trouble-shooting and simplifying the scripts and the issue seems to boil down to gravity darkening, not irradiation. Irradiation is just making up for the weird trend introduced by gravity darkening, which is why it seems to have "fixed it" in my case. In Phoebe we pass the gravity darkening coefficients directly to ellc under the assumption that the ellc implementation is similar, following surface flux ~ g^beta, where beta is the user-defined coefficient and approx 0.32 for convective and 1. for radiative envelope stars. Is this the case or is the assumption we make in Phoebe incorrect and we need to somehow modify the coefficients before passing them to ellc?

Here's a simple script to reproduce the issue. I'm using spherical stars for which gravity darkening shouldn't have any effect, but when setting gdc to something other than 0/None, there's a flux increase out of eclipse proportional to the gdc values.

import numpy as np
import matplotlib.pyplot as plt
import ellc

times = np.linspace(0,1,200)
lc_kwargs = {
             'radius_1': 0.2, 
             'radius_2': 0.2, 
             'sbratio': 1.0, 
             'incl': 90.,  
             't_zero': 0.0, 
             'period': 1, 
             'f_c': 0.,
             'f_s': 0., 
             'shape_1': 'sphere', 'shape_2': 'sphere', 
                'heat_1': None, 'heat_2': None,

}
lc_q1 = ellc.lc(times, q=1, gdc_1 = 0., gdc_2 = 0., **lc_kwargs)
lc_q2 = ellc.lc(times, q=1, gdc_1 = 0.5, gdc_2 = 0.5, **lc_kwargs)
lc_q3 = ellc.lc(times, q=1, gdc_1 = 1.0, gdc_2 = 1.0, **lc_kwargs)
period = 1.
phases = 0 % period / period
phases = (((times - phases * period) / period) % 1)
phases[phases > 0.5] -= 1
s = np.argsort(phases)

plt.plot(phases[s], lc_q1[s], label='gdc=0')
plt.plot(phases[s], lc_q2[s], label='gdc=0.5')
plt.plot(phases[s], lc_q3[s], label='gdc=1.0')
plt.legend() 

[pmaxted]

From Maxted (2016A&A...591A.111M) ...

For gravity darkening I assume that the specific intensity can be related to the local gravity by a power law with exponent y(l). Note that y(l) is a wavelength dependent quantity, not the bolometric gravity darkening exponent often used in other light curve models. Appropriate values of y(l) and limb-darkening coeffcients for various passbands can be found in Claret & Bloemen (2011).

[gecheline]

but the local gravity for a spherical star is the same everywhere, so in the simplest case I attached the light curve out of eclipse should be flat for any value of gdc, right?

[pmaxted]

The local gravity is not the same everywhere on the surface of a spherical star in a binary because of the gravitational potential from the companion.

[gecheline]

Oh, so ellc is still using the Roche potential for surface gravities even with spherical geometry of the stars (because that's not the case in Phoebe). That's good to know. It still doesn't explain why my system wouldn't converge without irradiation (irradiation really isn't that prominent as to make a huge difference on its own). And there's still a significant difference between the Phoebe-roche light curve and ellc-roche light curve, but I'll try to find the source of the discrepancy on the Phoebe-end. Thank you for your answers and clarifying this!