Temporal evolution, two zones, cooling+acc+adb exp

import warnings
warnings.filterwarnings('ignore')

import matplotlib.pyplot as plt
import numpy as np

import jetset
print('tested on jetset',jetset.__version__)

tested on jetset 1.2.2

In this tutorial I show how to perform a full acc+radiative+adiabatic expansion simulation. To have full understanding of the analysis presented in this tutorial, it is advised to read the paper Tramacere et al (2022) [Tramacere2022].

We load the model of the flare simulated in :ref:temp_ev_two_zone_cooling_acc. And the we evolve the radiative region under the effect of radiative plus adiabatic cooling

from jetset.jet_timedep import JetTimeEvol

temp_ev_acc=JetTimeEvol.load_model('two_zone_rad_acc.pkl')

temp_ev_acc.show_model()

--------------------------------------------------------------------------------
JetTimeEvol model description
--------------------------------------------------------------------------------

physical setup:

--------------------------------------------------------------------------------

Table length=29

name	par type	val	units	val*	units*	log
delta t	time	5.000000e+01	s	0.00029979245799999996	R/c	False
log. sampling	time	0.000000e+00		None		False
R/c	time	1.667820e+05	s	1.0	R/c	False
IC cooling		off		None		False
Sync cooling		on		None		False
Adiab. cooling		on		None		False
Reg. expansion		off		None		False
Diff coeff		6.666667e-06	s-1	None		False
Acc coeff		4.000000e-05	s-1	None		False
Diff index		2.000000e+00		None		False
Acc index		1.000000e+00	s-1	None		False
Tesc acc	time	5.003461e+04	s	3.0	R_acc/c	False
Eacc max	energy	4.000000e+60	erg	None		False
Tesc rad	time	1.667820e+65	s	1e+60	R/c	False
Delta R acc	accelerator_width	5.000000e+14	cm	None		False
B acc	magnetic field	2.000000e-01	cm	None		False
R_rad rad start	region_position	5.000000e+15	cm	None		False
R_H rad start	region_position	1.000000e+17	cm	None		False
T_A0=1/ACC_COEFF	time	2.500000e+04	s	0.149896229	R/c	False
T_D0=1/DIFF_COEFF	time	1.500000e+05	s	0.899377374	R/c	False
T_DA0=1/(2*DIFF_COEFF)	time	7.500000e+04	s	0.449688687	R/c	False
gamma Lambda Turb. max		1.173358e+11		None		False
gamma Lambda Coher. max		1.173358e+10		None		False
gamma eq Syst. Acc (synch. cool)		7.832383e+05		None		False
gamma eq Diff. Acc (synch. cool)		1.309535e+05		None		False
T cooling(gamma_eq=gamma_eq_Diff)		1.477242e+05	s	None		False
T cooling(gamma_eq=gamma_eq_Sys)		2.469874e+04	s	None		False
T min. synch. cooling		1.934500e+02	s	None		False
L inj (electrons)	injected lum.	5.000000e+39	erg/s	None		False

model parameters:

--------------------------------------------------------------------------------

Table length=30

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
jet_time_ev	duration	time_grid	s	1.000000e+06	0.000000e+00	--	False	True
jet_time_ev	gmin_grid	gamma_grid		1.000000e+00	0.000000e+00	--	False	True
jet_time_ev	gmax_grid	gamma_grid		1.000000e+08	0.000000e+00	--	False	True
jet_time_ev	gamma_grid_size	gamma_grid		1.500000e+03	0.000000e+00	--	False	True
jet_time_ev	TStart_Acc	time_grid	s	0.000000e+00	0.000000e+00	--	False	True
jet_time_ev	TStop_Acc	time_grid	s	1.000000e+05	0.000000e+00	--	False	True
jet_time_ev	TStart_Inj	time_grid	s	0.000000e+00	0.000000e+00	--	False	True
jet_time_ev	TStop_Inj	time_grid	s	1.000000e+05	0.000000e+00	--	False	True
jet_time_ev	T_esc_acc	escape_time	(R_acc/c)*	3.000000e+00	--	--	False	True
jet_time_ev	Esc_Index_acc	fp_coeff_index		0.000000e+00	--	--	False	True
jet_time_ev	t_D0	acceleration_time	s	1.500000e+05	0.000000e+00	--	False	True
jet_time_ev	t_A0	acceleration_time	s	2.500000e+04	0.000000e+00	--	False	True
jet_time_ev	Diff_Index	fp_coeff_index	s	2.000000e+00	0.000000e+00	--	False	True
jet_time_ev	Acc_Index	fp_coeff_index		1.000000e+00	--	--	False	True
jet_time_ev	Delta_R_acc	accelerator_width	cm	5.000000e+14	0.000000e+00	--	False	True
jet_time_ev	B_acc	magnetic_field	G	2.000000e-01	0.000000e+00	--	False	True
jet_time_ev	E_acc_max	acc_energy	erg	4.000000e+60	0.000000e+00	--	False	True
jet_time_ev	Lambda_max_Turb	turbulence_scale	cm	1.000000e+15	0.000000e+00	--	False	True
jet_time_ev	Lambda_choer_Turb_factor	turbulence_scale	cm	1.000000e-01	0.000000e+00	--	False	True
jet_time_ev	T_esc_rad	escape_time	(R/c)*	1.000000e+60	--	--	False	True
jet_time_ev	Esc_Index_rad	fp_coeff_index		0.000000e+00	--	--	False	True
jet_time_ev	R_rad_start	region_size	cm	5.000000e+15	0.000000e+00	--	False	True
jet_time_ev	R_H_rad_start	region_position	cm	1.000000e+17	0.000000e+00	--	False	True
jet_time_ev	m_B	magnetic_field_index		1.000000e+00	1.000000e+00	2.000000e+00	False	True
jet_time_ev	t_jet_exp	exp_start_time	s	1.000000e+05	0.000000e+00	--	False	True
jet_time_ev	beta_exp_R	beta_expansion	v/c*	1.000000e+00	0.000000e+00	1.000000e+00	False	True
jet_time_ev	B_rad	magnetic_field	G	2.000000e-01	0.000000e+00	--	False	True
jet_time_ev	t_size	time_grid		2.000000e+04	0.000000e+00	--	False	True
jet_time_ev	num_samples	time_ev_output		5.000000e+02	0.000000e+00	--	False	True
jet_time_ev	L_inj	inj_luminosity	erg / s	5.000000e+39	0.000000e+00	--	False	True

here we set some relevant parameters that will be described in detail in the next version of the documentation

temp_ev_acc.plot_time_profile()

<jetset.plot_sedfit.PlotTempEvDiagram at 0x7fef1f7baf10>

Particle spectrum in the radiative region

p=temp_ev_acc.plot_tempev_emitters(region='rad',loglog=False,energy_unit='gamma',pow=0)
p.ax.axvline(temp_ev_acc.temp_ev.gamma_eq_t_A, ls='--')
p.ax.axvline(temp_ev_acc.temp_ev.gamma_eq_t_DA, ls='--')
p.setlim(x_max=1E7,x_min=1,y_min=1E-18,y_max=100)

SEDs in the radiation region

p=temp_ev_acc.plot_tempev_model(region='rad',sed_data=None, use_cached = True)
p.setlim(y_min=1E-18,x_min=1E7)

We generate a lightcurve in the range nu1=2.4E22 Hz, nu2=7.2E25 Hz, without the effect of the light crossing time, in the observer frame

lg=temp_ev_acc.rad_region.make_lc(nu1=2.4E22,nu2=7.2E25,name='gamma',eval_cross_time=False,delta_t_out=100,use_cached=True,frame='obs')

plt.plot(lg['time'],lg['flux'])
plt.xlabel('time (%s)'%lg['time'].unit)
plt.ylabel('flux (%s)'%lg['flux'].unit)

Text(0, 0.5, 'flux (erg / (cm2 s))')

We generate a lightcurve in the range nu1=2.4E22 Hz, nu2=7.2E25 Hz, with the effect of the light crossing time, in the observer frame

lg_cross=temp_ev_acc.rad_region.make_lc(nu1=2.4E22,nu2=7.2E25,name='gamma',eval_cross_time=True,delta_t_out=100,use_cached=True,frame='obs',cross_time_slices=100)

plt.plot(lg['time'],lg['flux'])
plt.plot(lg_cross['time'],lg_cross['flux'])

plt.xlabel('time (%s)'%lg['time'].unit)
plt.ylabel('flux (%s)'%lg['flux'].unit)

Text(0, 0.5, 'flux (erg / (cm2 s))')

lr_1=temp_ev_acc.rad_region.make_lc(nu1=1E10,name='1E10 Hz',eval_cross_time=False,delta_t_out=100,use_cached=True,frame='obs')
lr_2=temp_ev_acc.rad_region.make_lc(nu1=5E9,name='1E9 Hz',eval_cross_time=False,delta_t_out=100,use_cached=True,frame='obs')

plt.plot(lr_1['time'],lr_1['flux']/lr_1['flux'].max())
plt.plot(lr_2['time'],lr_2['flux']/lr_2['flux'].max())

plt.xlabel('time (%s)'%lr_1['time'].unit)

Text(0.5, 0, 'time (s)')

lr_1_cross=temp_ev_acc.rad_region.make_lc(nu1=1E10,name='gamma',eval_cross_time=True,delta_t_out=100,use_cached=True,frame='obs',cross_time_slices=100)
lr_2_cross=temp_ev_acc.rad_region.make_lc(nu1=5E9,name='gamma',eval_cross_time=True,delta_t_out=100,use_cached=True,frame='obs',cross_time_slices=100)

plt.plot(lr_1_cross['time'],lr_1_cross['flux']/lr_1_cross['flux'].max())
plt.plot(lr_2_cross['time'],lr_2_cross['flux']/lr_2_cross['flux'].max())

plt.xlabel('time (%s)'%lr_1_cross['time'].unit)

Text(0.5, 0, 'time (s)')

Expanding the radiative region

We now plug the radiative region from temp_ev_acc to new model with adiabatic expansion

the following two functions define an estimate of the total extent of the simulation to follow the expansion

def delta_t_est(t_exp,R0,beta_exp):
    return t_exp+R0/(beta_exp*3E10)

def t_dec_est(R0,a,beta_exp):
    return ((R0+beta_exp*3E10)*np.power(beta_exp*3E10,a))

we set the initial radius equal to the radius of the radiative region of the temp_ev_acc model

t_exp=1E7
beta_exp=0.3
R0=temp_ev_acc.rad_region.jet.parameters.R.val
duration=delta_t_est(t_exp,R0,beta_exp)+10*t_dec_est(R0,-1,beta_exp)

we build the temp_ev_expansion expansion model

from jetset.jet_timedep import JetTimeEvol
temp_ev_expansion=JetTimeEvol(jet_rad=temp_ev_acc.rad_region.jet,inplace=True,only_radiation=True,Q_inj=None)

temp_ev_expansion.rad_region.jet.nu_min=1E8
T_SIZE=np.int(duration/1000)
NUM_SET=np.int(T_SIZE)
NUM_SET=min(5000,NUM_SET)


temp_ev_expansion.parameters.TStart_Inj.val=-0
temp_ev_expansion.parameters.TStop_Inj.val=-0

temp_ev_expansion.parameters.duration.val=duration
temp_ev_expansion.parameters.T_esc_rad.val=1E60
temp_ev_expansion.parameters.Esc_Index_rad.val=0
temp_ev_expansion.parameters.t_size.val=T_SIZE
temp_ev_expansion.parameters.num_samples.val=NUM_SET


temp_ev_expansion.parameters.gmin_grid.val=1.0
temp_ev_expansion.parameters.gmax_grid.val=1E8
temp_ev_expansion.parameters.gamma_grid_size.val=1500

we set to 'on' the region expansion, and we set the relevant paramters

temp_ev_expansion.region_expansion='on'
temp_ev_expansion.parameters.t_jet_exp.val=t_exp
temp_ev_expansion.parameters.beta_exp_R.val = beta_exp
temp_ev_expansion.parameters.R_rad_start.val = R0

temp_ev_expansion.init_TempEv()
temp_ev_expansion.show_model()

--------------------------------------------------------------------------------
JetTimeEvol model description
--------------------------------------------------------------------------------

physical setup:

--------------------------------------------------------------------------------

Table length=12

name	par type	val	units	val*	units*	log
delta t	time	1.000008e+03	s	0.005995894232556255	R/c	False
log. sampling	time	0.000000e+00		None		False
R/c	time	1.667820e+05	s	1.0	R/c	False
IC cooling		off		None		False
Sync cooling		on		None		False
Adiab. cooling		on		None		False
Reg. expansion		on		None		False
Tesc rad	time	1.667820e+65	s	1e+60	R/c	False
R_rad rad start	region_position	5.000000e+15	cm	None		False
R_H rad start	region_position	1.000000e+17	cm	None		False
beta exp.	region_position	3.000000e-01	v/c	8993773740.0 cm / s	cm/s	False
T min. synch. cooling		1.934500e+02	s	None		False

model parameters:

--------------------------------------------------------------------------------

Table length=17

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
jet_time_ev	duration	time_grid	s	1.611112e+07	0.000000e+00	--	False	True
jet_time_ev	gmin_grid	gamma_grid		1.000000e+00	0.000000e+00	--	False	True
jet_time_ev	gmax_grid	gamma_grid		1.000000e+08	0.000000e+00	--	False	True
jet_time_ev	gamma_grid_size	gamma_grid		1.500000e+03	0.000000e+00	--	False	True
jet_time_ev	TStart_Inj	time_grid	s	0.000000e+00	0.000000e+00	--	False	True
jet_time_ev	TStop_Inj	time_grid	s	0.000000e+00	0.000000e+00	--	False	True
jet_time_ev	T_esc_rad	escape_time	(R/c)*	1.000000e+60	--	--	False	True
jet_time_ev	Esc_Index_rad	fp_coeff_index		0.000000e+00	--	--	False	True
jet_time_ev	R_rad_start	region_size	cm	5.000000e+15	0.000000e+00	--	False	True
jet_time_ev	R_H_rad_start	region_position	cm	1.000000e+17	0.000000e+00	--	False	True
jet_time_ev	m_B	magnetic_field_index		1.000000e+00	1.000000e+00	2.000000e+00	False	True
jet_time_ev	t_jet_exp	exp_start_time	s	1.000000e+07	0.000000e+00	--	False	True
jet_time_ev	beta_exp_R	beta_expansion	v/c*	3.000000e-01	0.000000e+00	1.000000e+00	False	True
jet_time_ev	B_rad	magnetic_field	G	2.000000e-01	0.000000e+00	--	False	True
jet_time_ev	t_size	time_grid		1.611100e+04	0.000000e+00	--	False	True
jet_time_ev	num_samples	time_ev_output		5.000000e+03	0.000000e+00	--	False	True
jet_time_ev	L_inj	inj_luminosity	erg / s	1.000000e+39	0.000000e+00	--	False	True

temp_ev_expansion.plot_time_profile()

<jetset.plot_sedfit.PlotTempEvDiagram at 0x7fef0d4593a0>

we set ``do_injection=False`` because we want only to evolve the particle already injected and evolved in the radiative region of the ``temp_ev_acc`` model

setting cache_SEDs_rad=True will generate and cache all the SED at any time of the NUM_SET. This will increase the computational time during the run. Anyhow, will speed up the computation of SEDs and light curves. Moreover, these SEDs will be saved in the model, and will be read if once you will load the model in the future.

temp_ev_expansion.run(cache_SEDs_rad=True,do_injection=False)

temporal evolution running

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

temporal evolution completed
caching SED for each saved distribution: start

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

caching SED for each saved distribution: done

we now evaluate light curves, and plot the combination of the flare and adiabatic expansion simulations, for both the radio and gamma

lr_1_exp=temp_ev_expansion.rad_region.make_lc(nu1=1E10,name='1E10 Hz',eval_cross_time=False,delta_t_out=100,use_cached=True,frame='obs')
lr_2_exp=temp_ev_expansion.rad_region.make_lc(nu1=5E9,name='1E9 Hz',eval_cross_time=False,delta_t_out=100,use_cached=True,frame='obs')
lr_1_exp['time']+=lr_1['time'][-1]
lr_2_exp['time']+=lr_2['time'][-1]

lg_exp=temp_ev_expansion.rad_region.make_lc(nu1=2.4E22,nu2=7.2E25,name='gamma',eval_cross_time=False,delta_t_out=100,use_cached=True,frame='obs')
lg=temp_ev_acc.rad_region.make_lc(nu1=2.4E22,nu2=7.2E25,name='gamma',eval_cross_time=False,delta_t_out=100,use_cached=True,frame='obs')
lg_exp['time']+=lg['time'][-1]

plt.plot(lr_1['time'],lr_1['flux']/lr_1_exp['flux'].max(),c='b')
plt.plot(lr_2['time'],lr_2['flux']/lr_2_exp['flux'].max(),c='g')

plt.plot(lr_1_exp['time'],lr_1_exp['flux']/lr_1_exp['flux'].max(),label='10 GHz',c='b')
plt.plot(lr_2_exp['time'],lr_2_exp['flux']/lr_2_exp['flux'].max(),label='1 GHz',c='g')
plt.plot(lg['time'],lg['flux']/lg['flux'].max(),c='purple',label='gamma')
plt.plot(lg_exp['time'],lg_exp['flux']/lg['flux'].max(),c='purple')
plt.xlabel('time (%s)'%lr_1['time'].unit)
plt.legend()

<matplotlib.legend.Legend at 0x7feee3d7a2e0>

we notice the two peaks in the radio lightcurves, due to transition of the SSA frequency generated by the expansion (see [Tramacere2022] for more details)

p=temp_ev_expansion.plot_tempev_model(region='rad',sed_data=None, use_cached = True,time_slice_bin=50)
p.setlim(y_min=1E-18,x_min=1E7)

from jetset.plot_sedfit import PlotSED

p=PlotSED(frame='obs',density=False)
p.resplot.remove()
skip_label=False

step=int(temp_ev_expansion.parameters.num_samples.val/50)

for i in  range(0,NUM_SET,step):
    t=temp_ev_expansion.rad_region.time_sampled_emitters._get_time_samples(time_slice=i)
    s=temp_ev_expansion.rad_region.get_SED(comp='Sum',time_slice=i,frame='obs',use_cached=True)
    s_sync=temp_ev_expansion.rad_region.get_SED(comp='Sync',time_slice=i,frame='obs',use_cached=True)
    s_IC=temp_ev_expansion.rad_region.get_SED(comp='SSC',time_slice=i,frame='obs',use_cached=True)

    if t[0][0]<temp_ev_expansion.parameters.t_jet_exp.val:
        c='C0'
    else:
        c='C1'
    label=None
    if i==0:
        label='pre expansion'
    if t[0][0]>=temp_ev_expansion.parameters.t_jet_exp.val and skip_label is False:
        label='expansion'
        skip_label=True
    p.add_model_plot(model=s,label=label,color=c,density=False,auto_label=False)

p.setlim(y_min=1E-18,x_min=1E7)