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.3.0rc9
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')
===> setting C threads to 12
===> setting C threads to 12
temp_ev_acc.show_model()
--------------------------------------------------------------------------------
JetTimeEvol model description
--------------------------------------------------------------------------------
physical setup:
--------------------------------------------------------------------------------
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:
--------------------------------------------------------------------------------
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 0x11259a4d0>
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 / (s cm2))')
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 / (s cm2))')
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
Important
if you use a jet
model with R
depending (i.e. you used JetBase.make_conical_jet()
) to perform a temporal evolution (in the JetTimeEvol
class), the dependencies on R
will be removed, and to have R
dependent on the position across the jet axis, use the parameter beta_exp_R
in the JetTimeEvol
instead. In the next release a more flexible and direct approach will be provided.
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.int32(duration/1000)
NUM_SET=np.int32(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
===> setting C threads to 12
===> setting C threads to 12
===> setting C threads to 12
===> setting C threads to 12
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:
--------------------------------------------------------------------------------
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:
--------------------------------------------------------------------------------
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 0x15156bf40>
Note
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 0x10d0dfb80>
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)