Validation of the pp equilibrium against the Fokker-Plank equation solution#
In this tutorial we validate the integral solution for the \(e^{\pm}\) equilibrium, used for the pp jet against the Fokker-Plank equation solution implemented in the JetTimeEvol
class
(see the Hadronic pp jet model, for more details on the pp hadronic jet model)
Jet pp#
from jetset.jet_model import Jet
from jetset.jetkernel import jetkernel
from astropy import constants as const
from jetset.jet_emitters_factory import EmittersFactory
from jetset.jet_emitters import InjEmittersArrayDistribution
import numpy as np
import matplotlib.pyplot as plt
import jetset
print('tested on jetset',jetset.__version__)
tested on jetset 1.3.0rc7
j=Jet(emitters_distribution='plc',verbose=False,emitters_type='protons')
===> setting C threads to 12
j.parameters.z_cosm.val=z=0.001
j.parameters.beam_obj.val=1
j.parameters.gamma_cut.val=1000/(jetkernel.MPC2_TeV)
j.parameters.NH_pp.val=1
j.parameters.N.val=1
j.parameters.p.val=2.0
j.parameters.B.val=.5
j.parameters.R.val=1E18
j.parameters.gmin.val=1
j.parameters.gmax.val=1E4
j.set_emiss_lim(1E-60)
j.set_IC_nu_size(100)
j.gamma_grid_size=200
j.eval()
gmin=1.0/jetkernel.MPC2_TeV
j.set_N_from_U_emitters(1.0, gmin=gmin)
j.eval()
#j.show_model()
m=j.emitters_distribution.gamma_p>gmin
print('U N(p) p>1 TeV=%e erg/cm-3'%(jetkernel.MPC2*np.trapz(j.emitters_distribution.n_gamma_p[m]*j.emitters_distribution.gamma_p[m],j.emitters_distribution.gamma_p[m])))
U N(p) p>1 TeV=9.999992e-01 erg/cm-3
j.energetic_report(verbose=False)
%matplotlib inline
j.emitters_distribution.plot()
<jetset.plot_sedfit.PlotPdistr at 0x15357b7c0>
j.save_model('hadronic.pkl')
from jetset.jet_model import Jet
j=Jet.load_model('hadronic.pkl')
===> setting C threads to 12
setting up the JetTimeEvol model#
gamma_sec_evovled=np.copy(j.emitters_distribution.gamma_e)
n_gamma_sec_evovled=np.copy(j.emitters_distribution.n_gamma_e)
gamma_sec_inj=np.copy(j.emitters_distribution.gamma_e_second_inj)
n_gamma_sec_inj=np.copy(j.emitters_distribution.n_gamma_e_second_inj)
from jetset.jet_emitters_factory import EmittersFactory
from jetset.jet_emitters import InjEmittersArrayDistribution
q_inj=InjEmittersArrayDistribution(name='array_distr',emitters_type='electrons',gamma_array=gamma_sec_inj,n_gamma_array=n_gamma_sec_inj,normalize=False)
q_inj.parameters
name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
---|---|---|---|---|---|---|---|
gmin | low-energy-cut-off | lorentz-factor* | 1.000000e+00 | 1.000000e+00 | 1.000000e+09 | False | False |
gmax | high-energy-cut-off | lorentz-factor* | 1.836150e+07 | 1.000000e+00 | 1.000000e+15 | False | False |
Q | emitters_density | 1 / (s cm3) | 1.000000e+00 | 0.000000e+00 | -- | False | False |
None
%matplotlib inline
p=q_inj.plot()
p.ax.plot(gamma_sec_inj, n_gamma_sec_inj,'.',ms=1.5)
[<matplotlib.lines.Line2D at 0x154a84b20>]
from jetset.jet_timedep import JetTimeEvol
from jetset.jet_model import Jet
temp_ev=JetTimeEvol(jet_rad=j,Q_inj=q_inj,only_radiation=True,inplace=True)
===> setting C threads to 12
===> setting C threads to 12
/Users/orion/miniforge3/envs/jetset/lib/python3.10/site-packages/jetset/model_manager.py:158: UserWarning: no cosmology defined, using FlatLambdaCDM(name="Planck13", H0=67.77 km / (Mpc s), Om0=0.30712, Tcmb0=2.7255 K, Neff=3.046, m_nu=[0. 0. 0.06] eV, Ob0=0.048252)
warnings.warn(m)
temp_ev.Q_inj.parameters.Q.val
1
we use the acc region with escape time equal to radiative region
duration=5E9
duration_acc=0
T_SIZE=np.int32(2E6)
temp_ev.parameters.duration.val=duration
temp_ev.parameters.TStart_Inj.val=0
temp_ev.parameters.TStop_Inj.val=duration
temp_ev.parameters.T_esc_rad.val= 1
temp_ev.parameters.Esc_Index_rad.val=0
temp_ev.parameters.t_size.val=T_SIZE
temp_ev.parameters.num_samples.val=500
temp_ev.IC_cooling='off'
temp_ev.parameters.L_inj.val=0
temp_ev.parameters.gmin_grid.val=1.1
temp_ev.parameters.gmax_grid.val=5E7
temp_ev.parameters.gamma_grid_size.val=400
temp_ev.init_TempEv()
temp_ev.region_expansion='off'
temp_ev.show_model()
--------------------------------------------------------------------------------
JetTimeEvol model description
--------------------------------------------------------------------------------
physical setup:
--------------------------------------------------------------------------------
name | par type | val | units | val* | units* | log |
---|---|---|---|---|---|---|
delta t | time | 2.500000e+03 | s | 7.494811449999999e-05 | R/c | False |
log. sampling | time | 0.000000e+00 | None | False | ||
R/c | time | 3.335641e+07 | 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 | |||
Tesc rad | time | 3.335641e+07 | s | 1.0 | R/c | False |
R_rad rad start | region_position | 1.000000e+18 | cm | None | False | |
R_H rad start | region_position | 1.000000e+17 | cm | None | False | |
T min. synch. cooling | 6.190400e+01 | s | None | False | ||
L inj (electrons) | injected lum. | 7.490567e+38 | 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 | 5.000000e+09 | 0.000000e+00 | -- | False | True |
jet_time_ev | gmin_grid | gamma_grid | 1.100000e+00 | 0.000000e+00 | -- | False | True | |
jet_time_ev | gmax_grid | gamma_grid | 5.000000e+07 | 0.000000e+00 | -- | False | True | |
jet_time_ev | gamma_grid_size | gamma_grid | 4.000000e+02 | 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 | 5.000000e+09 | 0.000000e+00 | -- | False | True |
jet_time_ev | T_esc_rad | escape_time | (R/c)* | 1.000000e+00 | -- | -- | 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 | 1.000000e+18 | 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 | 5.000000e-01 | 0.000000e+00 | -- | False | True |
jet_time_ev | t_size | time_grid | 2.000000e+06 | 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 | 0.000000e+00 | 0.000000e+00 | -- | False | True |
p=temp_ev.plot_pre_run_plot(dpi=100)
p=temp_ev.plot_time_profile()
we perform the evoltion only for injection and cooling, withou acceleration region
temp_ev.run(only_injection=True,cache_SEDs_acc=False,do_injection=True,cache_SEDs_rad=False)
temporal evolution running
0%| | 0/2000000 [00:00<?, ?it/s]
temporal evolution completed
p=temp_ev.plot_tempev_emitters(region='rad',loglog=False,energy_unit='gamma',pow=0,plot_Q_inj=True)
p.ax.plot(gamma_sec_evovled,n_gamma_sec_evovled,'-',label='analytical solution',lw=4,color='gray',alpha=0.75,zorder=0)
p.ax.legend()
p.setlim(y_min=1E-30,y_max=1E-2)
m=n_gamma_sec_evovled>0
x_analytical=np.log10(gamma_sec_evovled[m])
y_analytical=np.log10(n_gamma_sec_evovled[m])
m=temp_ev.rad_region.time_sampled_emitters.n_gamma[-1]>0
x_num=np.log10(temp_ev.rad_region.time_sampled_emitters.gamma[m])
y_num=np.log10(temp_ev.rad_region.time_sampled_emitters.n_gamma[-1][m])
y_analytical_interp = np.interp(x_num, x_analytical,y_analytical, left=np.nan, right=np.nan)
m=~np.isnan(y_analytical_interp)
m=np.logical_and(m,x_num>0.25)
m=np.logical_and(m,x_num<6)
y_analytical_interp=10**y_analytical_interp[m]
x_out=x_num[m]
y_num=10**y_num[m]
d=np.fabs(y_analytical_interp-y_num)/y_num
assert(all(d<0.25))
plt.plot(x_out,d)
[<matplotlib.lines.Line2D at 0x15580e920>]