Temporal evolution, two zones, cooling+acc

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

This is a very preliminary documentation for the temporal evolution capabilities of jetset. Here we show how to create a decopuled radiative+acceleration region, and how to evolve the system in order to generate both particle spectra, SEDs, and lightcurves

To have full understanding of the analysis presented in this tutorial, it is advised to read the paper Tramacere et al. (2011) [Tramacere2011] for the understanding of stochastic and first order acceleration, and Tramacere et al (2022) [Tramacere2022] for the coupling of radiative and acceleration regions.

definition of the injected particle distribution (q_inj), and of the jet model for the radiative region

from jetset.jet_emitters_factory import InjEmittersFactory
from jetset.jet_model import Jet
jet_model=Jet()
q_inj=InjEmittersFactory().create_inj_emitters('pl',emitters_type='electrons',normalize=True)
q_inj.parameters.gmin.val=9
q_inj.parameters.gmax.val=10
q_inj.parameters.p.val=0.5

jet_model.parameters.beam_obj.val=30
jet_model.parameters.B.val=0.2
jet_model.parameters.z_cosm.val=0.03
jet_model.parameters.R.val=5E15

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

flare_duration=1.0E5
duration=flare_duration*10
t_D0=1.5E5
t_A0=2.5E4
T_esc_rad=1E60
L_inj=5.0E39
E_acc_max=4E60
Delta_R_acc_ratio=0.1
B_ratio=1.0
T_SIZE=2E4
NUM_SET=500
Diff_Index=2.0
Acc_Index=1.0

Here we instantiate the JetTimeEvol object, passing the radiative region jet model, and the injected particle class.

from jetset.jet_timedep import JetTimeEvol
temp_ev_acc=JetTimeEvol(jet_rad=jet_model,Q_inj=q_inj,inplace=True)

==> par: z_cosm from model: jet_leptonicacc_region linked to same parameter in model jet_leptonic

The IC cooling is switched off, as default, to make the process faster. to switch on the IC cooling temp_ev_acc.IC_cooling='on'

Now, we setup some relevant parameters

temp_ev_acc.rad_region.jet.nu_min=1E8
temp_ev_acc.acc_region.jet.nu_min=1E8
T_SIZE=np.int(T_SIZE)

if Delta_R_acc_ratio is not None:
    temp_ev_acc.parameters.Delta_R_acc.val=temp_ev_acc.parameters.R_rad_start.val*Delta_R_acc_ratio

T_esc_acc=t_A0/(temp_ev_acc.parameters.Delta_R_acc.val/3E10)*2



temp_ev_acc.parameters.duration.val=duration
temp_ev_acc.parameters.TStart_Acc.val=0
temp_ev_acc.parameters.TStop_Acc.val=flare_duration
temp_ev_acc.parameters.TStart_Inj.val=0
temp_ev_acc.parameters.TStop_Inj.val=flare_duration
temp_ev_acc.parameters.T_esc_acc.val=T_esc_acc
temp_ev_acc.parameters.T_esc_rad.val=T_esc_rad
temp_ev_acc.parameters.t_D0.val=t_D0
temp_ev_acc.parameters.t_A0.val=t_A0
temp_ev_acc.parameters.Esc_Index_acc.val=Diff_Index-2
temp_ev_acc.parameters.Esc_Index_rad.val=0
temp_ev_acc.parameters.Acc_Index.val=Acc_Index
temp_ev_acc.parameters.Diff_Index.val=Diff_Index
temp_ev_acc.parameters.t_size.val=T_SIZE
temp_ev_acc.parameters.num_samples.val=NUM_SET
temp_ev_acc.parameters.E_acc_max.val=E_acc_max
temp_ev_acc.parameters.L_inj.val=L_inj


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

temp_ev_acc.parameters.B_acc.val=temp_ev_acc.rad_region.jet.parameters.B.val*B_ratio
temp_ev_acc.init_TempEv()
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

temp_ev_acc.plot_time_profile()

<jetset.plot_sedfit.PlotTempEvDiagram at 0x7f9637ab46d0>

we do not want to evolve the particle in the ``jet_rad``, so we set ``only_injection=True``, and we set ``do_injection=True`` to injet the particle defined by ``q_inj``

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 read if you will reload the model in the future.

setting cache_SEDs_acc=True will generate and cache also the SEDs in the acceleration region.

only_injection=True
do_injection=True
plot_fit_model=True
plot_fit_distr=True
plot_emitters=True
plot_lcs=True
delta_t_out=1000
eval_cross_time=False
rest_frame='obs'
temp_ev_acc.run(only_injection=only_injection,
                do_injection=do_injection,
                cache_SEDs_acc=True,
                cache_SEDs_rad=True)

temporal evolution running

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

temporal evolution completed
caching SED for each saved distribution: start

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

caching SED for each saved distribution: done
caching SED for each saved distribution: start

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

caching SED for each saved distribution: done

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)

Particle spectrum in the acceleration region

p=temp_ev_acc.plot_tempev_emitters(region='acc',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-30,y_max=100)

SEDs in the acceleration 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)

SEDs in the acceleration region

p=temp_ev_acc.plot_tempev_model(region='acc',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')

lg

Table length=344

time	flux	R_blob	t_blob
s	erg / (cm2 s)	cm	s
float64	float64	float64	float64
0.0	0.0	5000000000000000.0	0.0
100.0	0.0	5000000000000000.0	2912.6213592233007
200.0	0.0	5000000000000000.0	5825.242718446601
300.0	4.4098133455386786e-86	5000000000000000.0	8737.864077669903
400.0	1.8338347214189153e-75	5000000000000000.0	11650.485436893203
500.0	4.619818537702253e-61	5000000000000000.0	14563.106796116504
600.0	4.074119989099911e-55	5000000000000000.0	17475.728155339806
700.0	4.480719706093064e-47	5000000000000000.0	20388.349514563106
800.0	3.859369921272289e-43	5000000000000000.0	23300.970873786406
...	...	...	...
33400.0	1.1745865571978132e-10	5000000000000000.0	972815.5339805826
33500.0	1.1673644350909526e-10	5000000000000000.0	975728.1553398059
33600.0	1.1601953475096658e-10	5000000000000000.0	978640.7766990291
33700.0	1.1530745106216505e-10	5000000000000000.0	981553.3980582524
33800.0	1.1460037834025703e-10	5000000000000000.0	984466.0194174757
33900.0	1.1389825349049919e-10	5000000000000000.0	987378.6407766991
34000.0	1.1320085680491812e-10	5000000000000000.0	990291.2621359223
34100.0	1.1250852647156014e-10	5000000000000000.0	993203.8834951456
34200.0	1.1182065068926023e-10	5000000000000000.0	996116.5048543689
34300.0	1.1113794126859807e-10	5000000000000000.0	999029.1262135921

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))')

We can save the model and reuse it later for plotting lightcurcves, SEDs, and electron distributions

temp_ev_acc.save_model('two_zone_rad_acc.pkl')

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

temp_ev_acc_1.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

p=temp_ev_acc_1.plot_tempev_model(region='rad',sed_data=None, use_cached = True)

lx=temp_ev_acc_1.rad_region.make_lc(nu1=1E17,nu2=1E18,name='X',eval_cross_time=False,delta_t_out=100,use_cached=True,frame='obs')
plt.plot(lx['time'],lx['flux'])
plt.xlabel('time (%s)'%lg['time'].unit)
plt.ylabel('flux (%s)'%lg['flux'].unit)

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