This Jupyter notebook outlines a model to simulate the evolution of a brown dwarf when influenced by a close M-dwarf companion. The model includes importing necessary libraries such as numpy, scipy, and matplotlib for numerical computations and plotting, along with a specific library for astronomical simulations. It details the initialization of constants like solar luminosity, radius, and mass, alongside brown dwarf parameters. It further describes the loading of an MDwarf model, the creation and normalization of various functions for the model, and the establishment of classes for modeling the interaction between a brown dwarf and its M-dwarf companion, specifically focusing on the aspects of flaring events, luminosity calculations, and temperature variations over time due to these stellar interactions. The model attempts to quantify how such interactions, particularly flaring events from the M-dwarf, could affect the temperature evolution of a brown dwarf, including detailed methods for calculating incident luminosity, orbital parameters, and resulting surface parameters and temperatures over significant time scales.

snippet.python

import numpy as np
import scipy
import matplotlib.pyplot as plt
from pysep.dsep import load_model
from scipy.interpolate import interp1d
from tqdm.notebook import tqdm
from scipy.optimize import curve_fit
 
try:
    plt.style.use(pubStyle)
except NameError:
    pass

snippet.python

SLUM = 3.846e33 # ergs / s
SRAD = 6.957e+10 # cm
SMASS = 1.989e+33 # g
BDLUM = 1/100_000 # L/Lsolar
BDTEMP = 1000 # K
BDRAD = 1.001e+10 # cm
sb = 5.670374419e-5 # cgs
kb =  1.3807e-16#  cm2 g s-2 K-1 (cgs
BDATMMFRAC = 0.00001  #0.000001
G = 6.67430e-8 # cgs

snippet.python

basedMDwarfModel = load_model("model.dsep")

snippet.python

def PosNormal(mean, sigma):
    x = np.random.normal(loc = mean, scale = sigma, size=1)[0]
    return(x if x>=0 else PosNormal(mean,sigma))

snippet.python

def normalize(a):
    if a.max() == a.min():
        return a
    return a - a.min()/(a.max()-a.min())

snippet.python

def interpTRK(targetAge, df):
    ages = df.AGE
    closestGreater = ages[ages > targetAge].min()
    closestLess = ages[ages < targetAge].max()
    row_closest_less = df[df['AGE'] == closestLess]
    row_closest_greater = df[df['AGE'] == closestGreater]
    fraction = (targetAge - closestLess) / (closestGreater - closestLess)
    interpolated_row = row_closest_less.drop(columns=['AGE']).reset_index(drop=True) + \
                   fraction * (row_closest_greater.drop(columns=['AGE']).reset_index(drop=True) - row_closest_less.drop(columns=['AGE']).reset_index(drop=True))
 
    return interpolated_row

snippet.python

def generate_flare_events_array(timeArraySize, flareProp):
    randomNumbers = np.random.rand(timeArraySize)
    flareEventsArray = randomNumbers < flareProp
 
    return flareEventsArray

snippet.python

def build_flare(t12):
    trise = t12[t12 <= 0]
    tdecay = t12[t12 > 0]
    frise = 1 + 1.941*trise - 0.175*trise**2 - 2.246*trise**3 - 1.1225 * trise**4
    fdecay = 0.6890 * np.exp(-1.6 * tdecay) + 0.3030*np.exp(-0.2783 * tdecay)
    return np.hstack([frise, fdecay])

snippet.python

class TriangulatedSphere:
    _SLUM = 3.846e33
    _SRAD = 6.957e+10
 
    def __init__(self, nTheta, nPhi, times, model, modelAge=5, noiseScale=0.1):
        trk = interpTRK(modelAge, model.trk)
        self.radius = 10**trk.log_R.values[0] * self._SRAD
        self.lum = 10**trk.log_L.values[0] * self._SLUM
 
        self.nTheta = nTheta
        self.nPhi = nPhi
        self.times = times
        self.fluxMap = np.zeros((times.shape[0], nTheta, nPhi))
 
 
        self.thetaValues = np.linspace(0, np.pi, nTheta )
        self.phiValues = np.linspace(0, 2 * np.pi, nPhi)
        self.surfaceArea = 4 * np.pi * self.radius**2
        self.triangleArea = self.surfaceArea / (nTheta * nPhi)
        self.meanFlux = (self.lum)/(4*np.pi * (self.radius)**2)
        self.set_overall_flux(self.meanFlux)
        self.noiseLevel = noiseScale*self.meanFlux
 
    def set_overall_flux(self, flux):
        """Set the same flux value for the entire surface."""
        self.fluxMap[:, :, :] = flux
 
    def set_overall_flux_at_time(self, flux, time):
        tID = self.time_to_index(time)
        self.fluxMap[tID, :, :] = flux
 
    def incriment_overall_flux_at_time(self, df, time):
        tID = self.time_to_index(time)
        self.fluxMap[tID, :, :] += df
 
    def set_flux_at(self, theta, phi, flux, time):
        iTheta, iPhi = self.theta_phi_to_index(theta, phi)
        tID = self.time_to_index(time)
        self.fluxMap[tID, iTheta, iPhi] = flux
 
    def get_flux_at(self, theta, phi, time):
        iTheta, iPhi = self.theta_phi_to_index(theta, phi)
        tID = self.time_to_index(time)
        return self.fluxMap[tID, iTheta, iPhi]
 
    def get_luminosity_over_solid_angle(self, thetaCenter, phiCenter, omega, time, solar=False):
        # Convert time to indices
        if isinstance(time, float):
            tIDs = np.array([self.time_to_index(time)])
        elif isinstance(time, slice):
            tIDs = np.array([self.time_to_index(t) for t in self.times[time]])
        elif isinstance(time, str) and time == '*':
            tIDs = np.arange(len(self.times))
 
        lums = np.zeros(len(tIDs))
 
        thetaGrid, phiGrid = np.meshgrid(self.thetaValues, self.phiValues, indexing='ij')
        distance = self.angular_distance_vectorized(thetaGrid, phiGrid, thetaCenter, phiCenter)
        withinSolidAngle = distance < np.sqrt(omega / np.pi)
 
        for i, tID in enumerate(tIDs):
            if solar:
                lums[i] = np.sum(self.fluxMap[tID, withinSolidAngle] * self.triangleArea / self._SLUM)
            else:
                lums[i] = np.sum(self.fluxMap[tID, withinSolidAngle] * self.triangleArea)
 
        return lums
 
    @staticmethod
    def angular_distance_vectorized(theta1, phi1, theta2, phi2):
        """Vectorized calculation of angular distance."""
        return np.arccos(np.sin(theta1) * np.sin(theta2) + np.cos(theta1) * np.cos(theta2) * np.cos(phi1 - phi2))
 
    def index_to_theta_phi(self, iTheta, iPhi):
        """Convert grid indices to theta and phi values."""
        return self.thetaValues[iTheta], self.phiValues[iPhi]
 
    def theta_phi_to_index(self, theta, phi):
        iTheta = np.argmin(np.abs(self.thetaValues - theta))
        iPhi = np.argmin(np.abs(self.phiValues - phi))
        return iTheta, iPhi
 
    def index_to_time(self, tID):
        return self.times[tID]
 
    def time_to_index(self, time):
        return np.argmin(np.abs(self.times - time))
 
    def evolve(self, flareProb, flareMeanDuration, flareMeanAmplitude, flareMeanLat, flareDurationScale, flareAmplitudeScale, flareLatScale):
        flareDecisions = np.random.uniform(0, 1, size=len(self.times)) < flareProb
        flareDurations = np.abs(np.random.normal(loc=flareMeanDuration, scale=flareDurationScale, size=flareDecisions.sum()))
        flareAmplitudes = np.abs(np.random.normal(loc=flareMeanAmplitude, scale=flareAmplitudeScale, size=flareDecisions.sum()))
 
        flarePhis = np.radians(np.random.choice([-1, 1], size=flareDecisions.sum())) * np.abs(np.random.normal(loc=flareMeanLat, scale=flareLatScale, size=flareDecisions.sum()))
        flareThetas = np.random.uniform(0, 2*np.pi, size=flareDecisions.sum())
 
        noise = np.random.normal(loc=0, scale=self.noiseLevel, size=len(self.times))
 
        fID = 0
        for tid, (time, doFlare, s) in tqdm(enumerate(zip(self.times, flareDecisions, noise)), total=len(self.times), desc="Evolving"):
            if doFlare:
                timeCond = (self.times >= time) & (self.times <= time + flareDurations[fID])
                flareTime = np.linspace(-1, 7, timeCond[timeCond==True].shape[0])
                flareData = build_flare(flareTime)              
                timeCond = (self.times >= time) & (self.times <= time + flareDurations[fID])
                flare = self.meanFlux + (self.meanFlux * flareAmplitudes[fID] * normalize(flareData))
                for flareTimeId, (flareFlux, dft) in enumerate(zip(flare, self.times[timeCond])):
                    self.set_flux_at(flareThetas[fID], flarePhis[fID], flareFlux, time + dft)
                fID += 1
 
            self.incriment_overall_flux_at_time(s, time)

snippet.python

estimateParticlesBBN = lambda mass: mass * (0.75 / 1.67e-27 + 0.25 / 6.64e-27 + 1e-9 / 1.15e-26)

snippet.python

def incident_luminosity(luminosityA, aR, bR, d):
    crossSectionalAreaB = np.pi * bR**2
    surfaceAreaD = 4 * np.pi * d**2
    fractionIncident = crossSectionalAreaB / surfaceAreaD
    incidentLuminosityB = luminosityA * fractionIncident
 
    return incidentLuminosityB

snippet.python

def orbital_params(m1, m2, p):
    Msolar = 1.989e+33 # g
    mt = m1 + m2
    mt = Msolar * mt
    G = 6.67430e-8 # cgs
    a = ((p**2 * G * mt)/(4*np.pi**2))**(1/3) # cm
    return a

$$ L = 4\pi R^{2}\sigma T^{4} $$ $$ U = \frac{3}{2} N k_{B} T $$

snippet.python

def calculate_new_surface_params(iL, temp, mass, dt):
    N = estimateParticlesBBN(mass*BDATMMFRAC) # assume all energy is distributed through the atmosphere
    Q = iL * dt
    dT = Q / ((3/2) * N * kb)
    return temp + dT

snippet.python

def cool_blackbody(radius, initialTemperature, mass, delta_time, internal_energy):
    N = estimateParticlesBBN(mass*BDATMMFRAC) # assume all energy is distributed through the atmosphere
    Ui = 3/2 * N * kb * initialTemperature
 
    initial_luminosity = 4 * np.pi * radius**2 * sb * (initialTemperature)**4
    energy_lost = initial_luminosity * delta_time
 
    Uf = (Ui - energy_lost) + internal_energy
    finalTemperature = (2/3 * Uf)/(N*kb)
    return finalTemperature

All units are kept in cgs (roughly, using pi 10 ^ 7 for sec - > year conversion)

snippet.python

# Define how many time steps are used. The cadence right now has a probelmatically large effect on the outcome. This needs
#  to be smoothed out.
period = 6 * 3600
orbits = 50
tmax = period * orbits
cadence = 30
times = np.arange(0, tmax, cadence)

One major issue right now is that the flaring model predicts flares based on the cadence of the data (i.e. a 5% flare rate means that every time step there is a 5 percent change of a flare event. This means that if you increase the observational frequency the number of flares will go up. This needs to be chaged to a flare rate / time)

snippet.python

sphere = TriangulatedSphere(model=basedMDwarfModel, times=times, nTheta=5, nPhi=5, noiseScale=0.001)
sphere.evolve(0.05, 1000*60, 0.25, 30, 5*60, 1, 10)
#sphere.evolve(0.05, 100*60, 0.01, 30, 5*60, 0.1, 10)
#flareProb, flareMeanDuration, flareMeanAmplitude, flareMeanLat, flareDurationScale, flareAmplitudeScale, flareLatScale):
 
# Assuming that the components are tidally locked we only need to subtend one hemisphere to get the incident luminosity at any given time
hemiLuminosity = sphere.get_luminosity_over_solid_angle(0, 0, omega=2*np.pi, time='*')

Evolving:   0%|          | 0/36000 [00:00<?, ?it/s]

snippet.python

fig, ax = plt.subplots(1,1,figsize=(10, 7))
ax.plot(times, hemiLuminosity)
ax.set_xlabel("Time [s]")
ax.set_ylabel("Luminsoity [erg s$^{-1}$]")

Text(0, 0.5, 'Luminsoity [erg s$^{-1}$]')

snippet.python

# The conversion from solar mass to g is done in the orbital params function since dsep tracks
# everyging in either solar units or log solar uints
a = orbital_params(basedMDwarfModel.mass, 0.075, period)
 
dt = (times[1] - times[0]) * 60
 
BrownDwarfTemps = list()
bdtemp = BDTEMP
bdmass = 0.075 * SMASS
 
# Here is a handwavey thing. I am basically assuming that all the energy
#  dumped onto the star stays in the atmosphere (and that the atmosphre is 
#  and idea monatomic gas (f=3))
N = estimateParticlesBBN(bdmass*BDATMMFRAC)
 
# this is used to figure out how much energy the star radiates away so
#  that some interneral energy fudge factor can be added to keep the
#  atmosphere from rapidly cooling to 0
initialCooling = bdtemp - cool_blackbody(BDRAD, bdtemp, bdmass, dt, 0)
 
# Due to the high temporal resolutions required to resolve flares and the long
# baseline needed to heat the star I run this section to just fine the mean 
# temperature change per orbit.
porbitID = 0
temps = list()
pOrbitTemp = bdtemp
for tid, (t, l) in enumerate(zip(times, hemiLuminosity)):
    orbitID = t // period
    if orbitID != porbitID:
        temps.append(bdtemp - pOrbitTemp)
        pOrbitTemp = bdtemp
        porbitID = orbitID
    incidentUponBD = incident_luminosity(l, sphere.radius, 1.001e+10, a)
    bdtemp = calculate_new_surface_params(incidentUponBD, bdtemp, bdmass, dt)
    BrownDwarfTemps.append(bdtemp)

snippet.python

line = lambda x, m, b: m*x + b
fit, ax = curve_fit(line, range(len(temps)), temps)
plt.plot(temps)
plt.plot(range(len(temps)), line(range(len(temps)), *fit))

[<matplotlib.lines.Line2D at 0x7f75cf909370>]

Now we will assume that this mean temperature change per orbit (on the hot side) holdes true for ever and evolve this over many many more orbits. We make the assumption that the hot side temperature is independent of the cold side temperature so we first calculate just the hot side temperature|. This allows us to more easily use a lookback time for the cold side temperature latter

snippet.python

meanTempChangePerOrbit = np.mean(temps)
 
morbits = 1e9
orbitJump = 5e2
 
orbits = np.arange(0, morbits, orbitJump)
 
hotSideTemp = np.zeros_like(orbits)
coldSideTemp = np.zeros_like(orbits)
 
hottemp = BDTEMP
coldtemp = BDTEMP
 
period = 6*3600 #seconds
longTimes = np.linspace(0, morbits*period, orbits.shape[0])
 
bdmass = 0.075 * SMASS
intrinsicCoolingCoef = 0.9
initialCooling = hottemp - cool_blackbody(BDRAD, hottemp, bdmass, period*orbitJump, 0)
internalEnergy = (3/2) * N * kb * initialCooling * intrinsicCoolingCoef
 
 
for oID, orbit in tqdm(enumerate(orbits), total=orbits.shape[0]):
    hottemp += (orbitID + 1) * meanTempChangePerOrbit
    hottemp = cool_blackbody(BDRAD, hottemp, bdmass, orbitJump * period, internalEnergy) 
    hotSideTemp[oID] = hottemp

  0%|          | 0/2000000 [00:00<?, ?it/s]

snippet.python

fig, ax = plt.subplots(1,1,figsize=(10, 7))
ax.plot(longTimes/(np.pi*1e7 * 1e6), hotSideTemp)
ax.set_xlabel("Time [Myr]", fontsize=27)
ax.set_ylabel("Hot Side Teff [K]", fontsize=27)

Text(0, 0.5, 'Hot Side Teff [K]')

snippet.python

# we then build an interpolator so that we can look back on the hot side temp from before
hst = interp1d(longTimes, hotSideTemp)

snippet.python

tauCO = 70*24*3600 # s (wright et al. 2018 -- calibrated for mid M dwarfs, likeley not correct for Brown Dwarfs)
offset = 300 # to allow for thermal transfer time # need to connect this to the convective overturn timescale or better yet dont be stupid about it
coldSideTemp = np.zeros_like(orbits) + BDTEMP
for oID, orbit in tqdm(enumerate(orbits[offset:]), total=orbits[offset:].shape[0]):
    delayedHotSideTemp = hst(oID*period)
    coldSideTemp[oID + offset] = delayedHotSideTemp
    coldSideTemp[oID + offset] = cool_blackbody(BDRAD, coldSideTemp[oID + offset], bdmass, orbitJump * period, internalEnergy) 

  0%|          | 0/1999700 [00:00<?, ?it/s]

snippet.python

fig, ax = plt.subplots(1,1,figsize=(10, 7))
ax.plot(longTimes/(np.pi*1e7*1e6), hotSideTemp, label="Hot Side")
ax.plot(longTimes/(np.pi*1e7*1e6), coldSideTemp, label="Cold Side")
ax.set_xlabel("Time [Myr]", fontsize=27)
ax.set_ylabel("Teff [K]", fontsize=27)
ax.legend()

<matplotlib.legend.Legend at 0x7f75cfd272b0>

snippet.python