In [ ]:
%matplotlib inline
from __future__ import print_function, division
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
import scipy.constants as sc
# class to read in dStar output files
from reader import dStarCrustReader, dStarReader
# plot format
charsize=14
major_ticklength=0.6*charsize
major_tickwidth=0.9
minor_ticklength=0.3*charsize
minor_tickwidth=0.7
rc('mathtext',**{'fontset':'stixsans'})
rc('font',**{'size':charsize,'sans-serif':'Bitstream Vera Sans'})
rc('axes',**{'titlesize':charsize,'labelsize':charsize})
rc('xtick',**{'major.size':major_ticklength,
'major.width':major_tickwidth,'minor.size':minor_ticklength,
'minor.width':minor_tickwidth,'labelsize':charsize})
rc('ytick',**{'major.size':major_ticklength,'major.width':major_tickwidth,
'minor.size':minor_ticklength,'minor.width':minor_tickwidth,'labelsize':charsize})
First, we define some functions: a routine to compute limits for a logarithmic axis; a routine to compute the thermal time \begin{equation} \newcommand{\Teff}{T_{\mathrm{eff}}^\infty} \tau = \frac{1}{4}\left[ \int \left(\frac{\rho C}{K}\right)^{1/2}\mathrm{d}r\right]^2; \end{equation} and a routine to compute the radius from the area.
In [ ]:
def thermal_time(c):
'''
Returns an array tau containing the thermal diffusion time to the surface to the location
of each crust zone.
c: of type dStarCrustReader
'''
sc.day
P = c.crust.pressure[0,:]
rho = c.crust.density[0,:]
g = c.crust.gravity[0,:]
K = c.crust.K[0,:]
Cp = c.crust.cp[0,:]
kernel = np.sqrt(rho*Cp/K)/rho/g
dP = np.zeros_like(kernel)
dP[1:] = P[1:]-P[0:-1]
dP[0] = P[0]
tau = np.zeros_like(kernel)
tau[0] = 0.25*(kernel[0]*dP[0])**2 / sc.day
for i in range(1,len(tau)):
tau[i] = 0.25*(np.dot(kernel[0:i],dP[0:i]))**2 / sc.day
return tau
def radius(area):
return np.sqrt(area/4.0/np.pi)*1.0e-5
Next we define functions to make common plots: the surface effective temperature $\Teff$ during cooling and theh temperature in the crust during outburst.
In [ ]:
MeVfm3 = 10.0*sc.eV*1.0e6/sc.fermi**3
def loglimits(x,delta=0.05):
'''
Returns two limiting values that enclose the array x with a padding
of delta (in log-space).
x: array-like
'''
lx = np.log10(x)
lmin,lmax = lx.min(),lx.max()
d = delta*(lmax-lmin)
return 10.0**(lmin-d),10.0**(lmax+d)
def grayscale(n,light_to_dark=True,scale=0.9):
'''
returns an array of grayscale values (length `n`) ranging from
`scale` to 0; if `light_to_dark=False`, then the order is reversed.
'''
if light_to_dark:
return [str(scale*(n-i-1)/(n-1)) for i in range(n)]
else:
return [str(scale*i/(n-1)) for i in range(n)]
def cooling_plot(h,ePhi):
# store cooling times
t = h.data.time
quiescent = np.where(t > 0.0)
t = t[quiescent]
Teff = 10.0**h.data.lg_Teff[quiescent]*sc.k/sc.eV*ePhi
plt.xlabel(r'$t$ [d]')
plt.ylabel(r'$T_{\mathrm{eff}}^\infty$ [K]')
plt.xlim(loglimits(t))
plt.semilogx(t,Teff,color='k')
plt.tight_layout()
def outburst_plot(c):
tout = np.where(c.times <= 0.0)[0]
clrs = grayscale(len(tout))
T = c.crust.temperature[tout,:]*1.0e-9
plt.xlim(loglimits(P))
plt.xlabel(r'$P\;[\mathrm{MeV\,fm^{-3}}]$')
plt.ylabel(r'$T$ [GK]')
for i in range(T.shape[0]):
plt.semilogx(P,T[i,:],color=clrs[i])
plt.tight_layout()
As a first case, we set the impurity parameter $Q=100$. Of course, for such a large $Q$ the meaning of impurity breaks down; it is closer to the mark to say that the lattice is very disordered (amporphous?) in this case. As a result, the thermal conductivity is low and the thermal diffusion is slow.
In [ ]:
c = dStarCrustReader('LOGS-Q100-H0',dt=10)
crust = c.crust
R = radius(crust.area[0,:])
P = crust.pressure[0,:]
drip = np.where(crust.Xn[0,:] > 0.0)
Pnuc = P/MeVfm3
ePhi = crust.ePhi[0,0]
Before analyzing the cooling, we plot the radius as a function of pressure and highlight the inner crust.
In [ ]:
xmin,xmax = Pnuc.min(),Pnuc.max()
dx = 0.05*(np.log10(xmax/xmin))
xmin,xmax = 10.0**(np.log10(xmin)-dx),10.0**(np.log10(xmax)+dx)
ymin,ymax = 11.0,11.85
plt.xlim(xmin,xmax)
plt.ylim(ymin,ymax)
plt.semilogx(Pnuc,R,color='k',label='outer crust')
plt.semilogx(Pnuc[drip],R[drip],linewidth=4,color='0.4',solid_capstyle='round',
label='inner crust')
plt.xlabel(r'$P\;[\mathrm{MeV\,fm^{-3}}]$')
plt.ylabel(r'$R = \sqrt{\mathcal{A}/(4\pi)}$ [km]')
plt.legend(loc='lower left',frameon=False,fontsize='small')
plt.tight_layout()
plt.savefig('radius.pdf',format='pdf',facecolor='none',edgecolor='none')
Now we read in the history file and plot the observed surface effective temperature $\Teff$.
In [ ]:
h = dStarReader('LOGS-Q100-H0/history.data')
cooling_plot(h,ePhi)
plt.savefig('cooling-Q100-H0.pdf',format='pdf',facecolor='none',edgecolor='none')
The cooling timescale is indeed long, $> 1000\,\mathrm{d}$. Let us check the thermal time.
In [ ]:
tau = thermal_time(c)
P = c.crust.pressure[0,:]/MeVfm3
plt.xlim(loglimits(P))
plt.xlabel(r'$P\;[\mathrm{MeV\,fm^{-3}}]$')
plt.ylabel(r'$\tau$ [d]')
plt.loglog(P,tau,color='k')
plt.tight_layout()
plt.savefig('tau-Q100-H0.pdf',format='pdf',facecolor='none',edgecolor='none')
Finally, let us plot the temperature evolution through the outburst.
In [ ]:
outburst_plot(c)
plt.savefig('profile-Q100-H0.pdf',format='pdf',facecolor='none',edgecolor='none')
tout = np.where(c.times <= 0)[0][0]
T0 = c.crust.temperature[tout,:]*1.0e-9
P0 = P
In [ ]:
h4 = dStarReader('LOGS-Q4-H0/history.data')
cooling_plot(h4,ePhi)
plt.savefig('cooling-Q4-H0.pdf',format='pdf',facecolor='none',edgecolor='none')
Let's compare the thermal diffusion time with that for the high-impurity crust.
In [ ]:
c4 = dStarCrustReader('LOGS-Q4-H0')
P = c4.crust.pressure[0,:]/MeVfm3
plt.xlim(loglimits(P))
plt.xlabel(r'$P\;[\mathrm{MeV\,fm^{-3}}]$')
plt.ylabel(r'$\tau$ [d]')
plt.loglog(P,thermal_time(c4),color='k')
plt.loglog(P,thermal_time(c),linestyle=':',color='k')
plt.tight_layout()
plt.savefig('tau-Q4-H0.pdf',format='pdf',facecolor='none',edgecolor='none')
As expected from the evolution of the surface temperature, the thermal diffusion time is an order of magnitude less in this case. Because the thermal diffusion time is now less than the duration of the outburst, the crust has time to thermally relax; compare the profile of temperature in crust with that for the high-impurity case.
In [ ]:
outburst_plot(c4)
plt.savefig('profile-Q4-H0.pdf',format='pdf',facecolor='none',edgecolor='none')
This explains why $\Teff$ is so constant during the first 100 days: the outer layers of the star are already thermally relaxed, so the flux is unchanging until the inner crust has time to evolve. To see this, let's plot the temperature as a function of $\tau$.
In [ ]:
tout = np.where(c4.times == 0.0)[0][0]
indcs = [tout,tout+1,tout+3,tout+6,tout+8]
Ts = c4.crust.temperature[indcs,:]*1.0e-9
clrs = grayscale(len(indcs))
times = c4.times[indcs]
plt.xlim(loglimits(tau,delta=0.2))
plt.xlabel(r'$\tau$ [d]')
plt.ylabel(r'$T$ [GK]')
for i in range(len(indcs)):
plt.semilogx(tau,Ts[i,:],color=clrs[i])
if i > len(indcs)-3:
tstr = r'$t = {0:0.0f}$ d'.format(times[i])
plt.annotate(s=tstr,xy=(tau[0],Ts[i,0]),verticalalignment='center',
horizontalalignment='right',color='k',fontsize='small')
plt.tight_layout()
plt.savefig('Ttau-Q4-H0.pdf',format='pdf',facecolor='none',edgecolor='none')
To get $\Teff$ to cool at early times, we need the temperature to decrease with depth in the outer crust. We achieve this by adding a heat source, of strength $L = 1.7\,\mathrm{MeV}\times \mathrm{d}\dot{M}/\mathrm{d}t$.
In [ ]:
h17 = dStarReader('LOGS-Q4-H1.7/history.data')
cooling_plot(h17,ePhi)
plt.savefig('cooling-Q4-H1.7.pdf',format='pdf',facecolor='none',edgecolor='none')
In [ ]:
c17 = dStarCrustReader('LOGS-Q4-H1.7')
outburst_plot(c17)
plt.savefig('profile-Q4-H1.7.pdf',format='pdf',facecolor='none',edgecolor='none')
At a time $t$ after the end of the outburst, the crust has cooled down to a depth where $\tau \approx t$. The evolution of the surface temperature $\Teff$ is therefore mapping out the crust temperature as a function of depth.
In [ ]:
tout = np.where(c17.times == 0.0)[0][0]
indcs = [tout,tout+1,tout+3,tout+6,tout+8,tout+12,tout+16]
Ts = c17.crust.temperature[indcs,:]*1.0e-9
clrs = grayscale(len(indcs))
times = c17.times[indcs]
plt.xlim(loglimits(tau))
plt.xlabel(r'$\tau$ [d]')
plt.ylabel(r'$T$ [GK]')
for i in range(len(indcs)):
plt.semilogx(tau,Ts[i,:],color=clrs[i])
if i > len(indcs)-5:
tstr = r'$t = {0:0.0f}$ d'.format(times[i])
plt.annotate(s=tstr,xy=(1.0,Ts[i,0]),verticalalignment='bottom',
color='k',fontsize='small')
plt.tight_layout()
plt.savefig('Ttau-Q4-H1.7.pdf',format='pdf',facecolor='none',edgecolor='none')
In [ ]: