Overview of AM³

AM\(^3\) calculates the time-dependent evolution of the following particle species, assumed to be immersed in a homogeneous and isotropic region and magnetic field:

Particle	Designation in Python methods
\(\gamma\)	`photons`
\(\text{e}^+\) and \(\text{e}^-\)	`pairs`
\(\text{e}^-\)	`electrons`
\(\text{e}^+\)	`positrons`
\(\text{p}\)	`protons`
\(\text{n}\)	`neutrons`
\(\pi\) (both signs)	`pions`
\(\pi^+\)	`pion_plus`
\(\pi^-\)	`pion_minus`
\(\mu\) (both signs and handedness)	`muons`
\(\mu^-_\text{L}\)	`mu_minus_left`
\(\mu^-_\text{R}\)	`mu_minus_right`
\(\mu^+_\text{L}\)	`mu_plus_left`
\(\mu^+_\text{R}\)	`mu_plus_right`
\(\nu\) (all flavors)	`neutrinos`
\(\nu_{\mu}\)	`muon_neutrinos`
\(\bar{\nu}_\mu\)	`muon_anti_neutrinos`
\(\nu_\text{e}\)	`electron_neutrinos`
\(\bar{\nu}_\text{e}\)	`electron_anti_neutrinos`

The evolution of particle number density spectra \(n\) as a function of time \(t\) is calculated numerically by solving a system of time-dependent partial differential equations of the form

\[\partial_{t} n(E,t)=-\partial_{E} \dot{E}(E,t)n(E,t)-\alpha(E,t) n(E,t)+Q(E,t) \, .\]

This includes terms representing cooling [ \(\dot{E}(E,t)\) ], escape/sink [ \(\alpha(E,t)\) ] and injection/source [ \(Q(E,t)\) ]. For the computation the equations are re-written on a logarithmic energy grid, using \(x=\ln E/E_0\) with different energy scales that can be specified by the user (see below for the energy grid). Then, the above equation takes the form

\[ \partial_{t}n(x,t)=-\partial_{x}\left[A(x,t)n(x,t)\right]-\alpha(x,t)n(x,t)+\epsilon(x,t) \, , \label{equ:overalllogarithmicform} \]

with the substitutions

\[ x=\ln \frac{E}{E_0}, n(x)\equiv\frac{\partial^2 N}{\partial x \partial V} = E n(E), A(x)=\dfrac{\dot{E}(E)}{E}, \epsilon(x)=E Q(E), \alpha(x)=\alpha(E)\, .\]

Note that code-internally the signs of the cooling and skin terms are reversed.

The following physics processes are taken into account:

Process	Name in the Python methods	Corresponding terms in the equation system
Synchrotron	`syn`	\(Q_\gamma\), \(\alpha_\gamma\), \(\dot{E}\) for all charged particles
Inverse Compton	`compton`	\(Q_\gamma\), \(\alpha_\gamma\), \(\dot{E}\) for all charged particles
\(\gamma \gamma\) annihilation	`annnihilation`	\(\alpha_\gamma\), \(Q_{e^{\pm}}\)
Bethe-Heitler pair production	`bethe_heitler`	\(\dot{E}_p\), \(Q_{e^{\pm}}\)
Photo-pion production	`photopion`	\(\alpha_{\gamma}\), \(\alpha_p\), \(Q_{\pi^\pm}\), \(Q_\gamma\)
Proton-proton interactions	`proton_proton`	\(\alpha_p\), \(Q_{\pi^\pm}\), \(Q_\gamma\)
Expansion of the region.	`expansion`	\(\alpha\) for all particles
Adiabatic cooling	`adiabatic`	\(\dot{E}\) for all charged particles
Escape	`escape`	\(\alpha\) for all particles
Pion decay	`pion_decay`	\(\alpha_{\pi^\pm}\), \(Q_{\mu^\pm}\), \(Q_{\nu_\mu}\)
Muon decay	`muon_decay`	\(\alpha_{\mu^\pm}\), \(Q_{e^{\pm}}\), \(Q_{\nu_\mu}\), \(Q_{\nu_e}\)
Injection	`injection`	\(Q_{e^{-}}\), \(Q_{p}\)

Notes: (1) Due to their short lifetime, neutral pions are assumed to decay instantaneously. For photo-pion production we hence don’t list a source term for neutral pions but instead a photon source term. (2) Injection here refers to the build-in template injection. Arbitrary injection is possible by passing arrays.

For details on the implementation of the different processes we point to Klinger et al. (2023) [arXiv e-print 2312.13371] <https://arxiv.org/abs/2312.13371>_. The flexibility of AM\(^3\) is implemented by switches which allow to turn on/off processes and sub-processes. We provide a comprehensive list here.

Basic code structure

Because the Python interface allows the user to fully interact with the code features needed for performing siulations, the internal C++ structure is mostly relevant at the developer level. The C++ structure underlying each Python method can be traced back by consulting the Python interface file, under libpython/python_interface.cc.

In the following, we assume that AM\(^3\) has been compiled and imported as a Python module using import am3.

Running AM\(^3\) is possible through an user-interface class called AM3 which communicates with the underlying C++ structure. This class is also mirrored in the python interface. At the beginning of each run it is thus neccesarry to initialize an AM3 object:

am3 = am3.AM3()

Internally, each physics process is implemented in a separate C++ module, contained in .cc files such as BetheHeitler.cc. The solver methods are stored in Solver.cc, including the inverse Compton and synchrotron radiation losses of charged particles. In addition, SimulationManager.cc is the core class managing the simulation, PhysicsHandler is a helper class to collect all physics processes, RunParams holds the simulation parameters and AM3Arrays holds the arrays containing the particle number densities and other grid related information during the simulation.

Energy Grids

AM\(^3\) works with four different energy grids, all equally spaced in logarithmic space (natural base \(e\)). While the grid spacing is pre-defined to \(\Delta \ln (E/E_0) = 0.1\) in log-space (accessible through am3.get_log_energy_spacing()), the grids can be initialized at the beginning of the run specifying the minimum photon and neutrino energies, as well as the maximal energy of all particles:

am3.update_energy_grid(minimum photon energy [eV], minimum neutrino energy [eV], maximum energy of all particles [eV])

The proton, neutron, pion and muon grid has a minimum energy of \(1\) GeV (corresponding approximately to the proton rest mass). The electron grid has minimum energy \(m_e c^2\). If the grids are not explicitly initialized, they take the default values.

Photon grid. Access in eV with am3.get_egrid_photons().
Default: 620 bins, from 85 MHz up to 270 EeV.
“Leptonic” grid (electrons and positrons). Access values in eV with am3.get_egrid_lep().
Default: 340 bins, from \(m_e c^2\) up to 270 EeV.
“Hadronic” grid (protons, neutrons, pions, and muons). Access values in eV with am3.get_egrid_had().
Default: 314 bins, from \(m_p c^2\) up to 40 ZeV or 40,000 EeV.
“Neutrino” grid. Access values in eV with am3.get_egrid_neutrinos().
Default: 464 bins, from 300 eV up to 39 ZeV.

For all grids, the grid length may be accessed with am3.get_grid_size_<grid_type>, and the conversion between eV and the grid index with am3.eV_to_<grid_type>_grid_index, where <grid_type> can again be photons, lep, had and neutrinos.

Simulation Parameters

Energy-independent escape timescale \(t_\mathrm{esc}\). The escape timescale will typically represent the light-crossing time of the region. This parameter should be set before the beginning of the simulation using am3.set_escape_timescale(t_in_seconds). This defines an energy-independent escape timescale for all particle species. After this, custom escape timescales can be set for each individual particle species, including energy-dependent escape timescales, using the am3.set_t_<particle>_escape functions described below.
At any point during runtime, the am3.set_escape_timescale function can be called any number of times to change the escape timescale at that point in the simulation. This resets all particle escape timescales to the new energy-independent value 1./am3.get_escape_timescale(). Therefore, this operation overrides previous calls to the am3.set_t_<particle>_escape functions described below.\ Switch dependence: If the escape switch is off (see section on switches below), physical escape is not included in the solver, and calling am3.set_escape_timescale will not affect the simulation. Access current value in seconds with am3.get_escape_timescale().\
Energy-dependent escape timescales \(t_\mathrm{esc}(E)\). After calling am3.set_escape_timescale, all particle species are given an energy-independent escape timescale. After this, the escape timescale of each individual species can be changed using am3.set_t_<particle>_escape(timescale_array), where <particle> can be photon, pair, proton, neutron, pion, muon, or neutrino. The array timescale_array must be given in seconds and be of the same size as the energy array of the specific particle species. For example, am3.set_t_pair_escape(np.full(am3.get_egrid_lep().size, 1e6)) sets an energy-independent escape timescale of \(10^6\) s for electrons and positrons.\ The am3.set_t_<particle>_escape functions can be called at any point and any number of times during the simulation.
Switch dependence: If the escape switch is off (see section on switches below), physical escape is not included in the solver, and calling am3.set_t_<particle>_escape will not affect the simulation. Access array of current escape timescales in seconds with am3.get_t_<particle>_escape().\
Expansion timescale \(t_\mathrm{exp}\). This affects (1) the expansion of the volume of the simulated region with time, which leads to a gradual dilution of the particle densities of all species; (2) energy losses by charged particles due to adiabatic cooling, due to their electromagnetic coupling to the expanding plasma.
It should be set before the beginning of the simulation with am3.set_expansion_timescale(t_in_seconds). It may also be changed at any point during runtime.
Switch dependence: volume expansion and adiabatic cooling of charged particles are only included in the solver when the respective switches are on (see section on process switches below). If both switches are off, the expansion timescale will not play a role in the simulation. Access current value in seconds with am3.get_expansion_timescale().\
Magnetic field strength \(B\). This turbulent magnetic field is considered homogeneous and isotropic.
It should be set before the beginning of the simulation with am3.set_mag_field(b_in_gauss). It may also be changed at any point during runtime with the same set function. Every time this function is called, the kernels that depend on the magnetic field are automatically re-calculated with the updated magnetic field value. Access current value in Gauss using am3.get_mag_field()\
Density of target protons for pp interactions \(n_\mathrm{pp}\). This is considered as a homogeneous and isotropic distribution of non-relativistic thermal protons, and is only relevant when pp interactions are turned on (am3.set_process_pp(1), see List of switches).
If accounting for proton-proton interactions, this parameter should be set before the beginning of the simulation with am3.set_pp_target_proton_density(protons_per_cm3). It may also be changed at any point during runtime with the same set function. Access current value in protons per cubic centimeter using am3.get_pp_target_proton_density() (default: 0.0).\
Time step of the numerical solver \(\Delta t\). This controls the precision of the solver in the time domain.
The default value is a rather conservative \(10^{-3}t_\mathrm{esc}\).
It may be adjusted before the beginning of the simulation via am3.set_solver_time_step(step_in_seconds). It may also be changed at any point during runtime, which allows to implement also non-linear time-stepping such as a logarithmically growing time step.
For a typical AGN simulation, this may typically be increased for the sake of performance. On the contrary, for very extreme environments with high optical thickness, the time step may have to be reduced to avoid imprecisions derived from hardness problems. Access current value in seconds with am3.get_solver_time_step().\

For all set functions introduced above the current settings may be retrieved replacing set with get in the function name.

Accessing the current particle density spectra

The current particle density spectra can be accessed in the form of Numpy arrays in units of \(\text{cm}^{-3}\), using

am3.get_<particle>()

where <particle> represents one of the Python particle designations provided at the beginning of this overview. For example, am3.get_photons() returns the current photon density spectrum per cubic centimeter. Before the first step of the time-dependent solver is performed, all particle densities are initialized to zero.

The particle densities can also be set at any moment during runtime using

am3.set_current_densities_<particle>(array_particles_per_cm3)

Note, however, that this does not simulate particle injection, since it does not implement any source term in the equations. Instead, it sipmply re-defines the density spectrum of a given particle species at the current time and it overwrites any spectrum self-consistently emitted until that point. For particle injection methods, see below.

Accessing particle emission from a specific process

When switching on am3.set_process_parse_sed(True), it is possible to access the current density spectrum of secondaries produced by specific processes, following the structure:

am3.get_<paritcle>_<process>() # array with particles per cubic centimeter

For pions, muons and neutrinos, <process> can be photopion, or proton_proton. For electron-positron pairs, <process> can be injected, annihilation, bethe_heitler, photopion, or proton_proton. For photons, we implemented the following options:

Command	`<particle>` options	`<process>` options	Comment
`am3.get_photons()`			Total photon density, including external photon fields.
`am3.get_photons_except_injected()`			Total density of non-thermal photons, but not those injected by the user.
`am3.get_photons_<particle>_<process>()`	`injected_electrons`, `annihilation_pairs`, `bethe_heitler_pairs`, `photo_meson_pairs`, `pions`, `muons`	`syn` (synchrotron), `compton` (inverse Compton), `syn_compton` (sum of both)
`am3.get_photons_pi0_decay_<process>()`	-	`photopion` `proton_proton`
`am3.get_<neutrino_species>_<process>()`	`neutrinos`, `muon_neutrinos`, `muon_anti_neutrinos`, `electron_neutrinos`, `electron_anti_neutrinos`,	`photopion` `proton_proton`

For example, by calling

am3.get_photons_photo_meson_pairs_compton() # array in photons per cubic centimeter 

we retrieve an array with the density spectrum of inverse Compton radiation emitted by electron/positron pairs produced via the photo-meson process. For neutrinos we call

am3.get_neutrinos_proton_proton() # array in neutirnos per cubic centimeter 

to get neutrinos of all flavours from proton-proton interactions.

Accessing interaction timescales

Interaction timescales can be accessed with the following syntax:

am3.get_t_<particle>_<process>() # array with timescales in seconds

The possible <particle> and <process> terms are listed below (note the name pair for electrons or positrons).

`<particle>` options	`<process>` options	Definition
`photon`	`ssa`	Attenuation due to synchrotron self-absorption
	`compton`	Upscattering due to inverse Compton
	`escape`	Physical escape
	`annihilation`	Annihilation into pairs
	`photopion`	Loss timescale due to photopion interaction
	`emission_from_syn`	Total synchrotron emission timescale \(n/Q\)
	`emission_from_ic`	Total inverse Compton emission timescale \(n/Q\)
	`emission_from_pi0`	Total emission timescale from neutral pion decay \(n/Q\)

`pair`	`syn`	Synchrotron interaction timescale
	`compton`	Inverse Compton interaction timescale
	`acceleration`	Acceleration timescale (see below)
	`adiabatic`	Adiabatic cooling timescale
	`escape`	Sink timescale due to phsyical escape
	`emission_from_annihilation`	Time scale of pair emission from photon-photon annihilation \(n/Q\)
	`emission_from_injection`	Time scale of pair emission from direct injection \(n/Q\)
	`emission_from_bethe_heitler`	Time scale of pair emission from Bethe-Heitler \(n/Q\)
	`emission_from_photopion`	Time scale of pair emission from photopion interactions \(n/Q\)
	`emission_from_proton_proton`	Time scale of pair emission from proton-proton interactions \(n/Q\)

`proton`	`syn`	Synchrotron interaction timescale
	`compton`	Inverse Compton interaction timescale
	`acceleration`	Acceleration timescale (see below)
	`adiabatic`	Adiabatic cooling timescale
	`escape`	Sink timescale due to phsyical escape
	`bethe_heitler`	Proton Bethe-Heitler interaction
	`photopion`	Proton photopion interaction
	`pp`	Proton-proton interactions

`pion`	`syn`	Synchrotron interaction timescale
	`compton`	Inverse Compton interaction timescale
	`escape`	Sink timescale due to phsyical escape
	`decay`	Sink timescale due to decay

`muon`	`syn`	Synchrotron interaction timescale
	`compton`	Inverse Compton interaction timescale
	`escape`	Sink timescale due to phsyical escape
	`decay`	Sink timescale due to decay

`neutron`	`escape`	Sink timescale due to phsyical escape

`neutrino`	`escape`	Sink timescale due to phsyical escape

As we can see, in addition to the processes that can be explicitly simulated, it is possible to retrieve the acceleration timescale for electrons and protons for purposes of plotting and comparing. It is defined as (note the definition of \(\eta\) as an inverse efficiency here)

\[t^{-1}_{\text{acc}} (E) = \frac{\eta e c B}{E}\]

for a particle of energy \(E\) and electric charge \(e\) in a magnetic field \(B\). The acceleration efficiency \(\eta\) can be set through am3.set_electron_accel_factor(1/eta) (for leptons) and am3.set_proton_accel_factor(1/eta) (for hadrons). From this, code-internally also the maximum energy can be determined, balancing losses and acceleration, using am3.estimate_max_electron_energy() and am3.estimate_max_proton_energy(). The current values for the acceleration factors may be retrieved replacing set by get, while the calculated maximal energies are obtained through am3.get_estimated_electron_emax() and am3.get_estimated_proton_emax()

Particle injection

To inject a given particle distribution, use the respective “set injection” methods:

am3.set_injection_rate_<particle>(array_in_particle_number_per_second_per_cm3)

where <particle> represents one of the names provided at the beginning of this overview. Generic names like pairs, muons, and pions are not available here, since the specific species must be specified (e.g. am3.set_injection_rate_pi0 or am3.set_injection_rate_mu_plus_right). The array size must correspond to the size of the energy grid of the respective species.

The current injection rates can be accessed through

am3.get_injection_rate_<particle>()

which returns an array containing the number of particles injected per cubic centimeter per second. This does not include injection due to particle interactions, but refers only to the injection rates set by the user.

In gerenal, this represents external injection of particles. In the case of electrons and protons, it emulates pre-acceleration, while in the case of photons, it emulates an external photon field, such as from a broad line region surrounding an AGN jet (see blazar example).

Injection rate arrays are initialized to zeros, and can be set prior to the beginning of the time-dependent simulation or at any point during the simulation. Once an injection spectrum is defined using the above method, an injection term is aded to the equations. Particles of that species will be injected into the system at the given rates at every time stap until the end of the simulation, or until the same method is called again with a different input.

Injection of protons and electrons using template broken power-law spectra

Additionally to the above-mentioned injeciton methods, cosmic-ray electrons and protons can also be injected in the form of a broken-power-law spectrum by means of template functions:

am3.set_powerlaw_injection_parameters_electrons(
  volume_in_cm3, 
  injection_luminosity_in_erg_per_s,
  emin_in_eV, 
  ebreak_in_eV,
  emax_in_eV,
  index_below_break, 
  index_above_break,
  cutoff_steepness)

and similarly for protons. For example, to inject a simple \(E^{-2}\) power-law proton spectrum ranging from cold protons up to 1 PeV with a simple exponential cutoff, we would use

light_speed = 3e10 # cm/s
volume = 4/3. * np.pi * (am3.get_escape_timescale * light_speed) ** 3 # cm
am3.set_powerlaw_injection_parameters_electrons(
  volume, # cm
  1e45,   # erg/s
  1e9,    # eV
  1e9,    # eV
  1e15,   # eV
  2.0,  
  2.0,
  1.0     # cutoff steepness
)

This method can be called at any time during the simulation. To use the internal estimation of the maximum energies for the template power-law spectra (which will overwrite your specified emax_in_eV), set am3.set_estimate_max_energies(1). Note that during the first timestep, AM\(^3\) will use the user-specified maximal energy, as the cooling rates are not yet computed internally.

Important: (1) Using the template for injection overwrites any injection of the same particle species set previously either using this method or am3.set_injection_rate_<particle>(). (2) Internally, AM\(^3\) only uses densities, not volumes. If the volume for the injection is expanding over time (and you specified the corresponding expansion timescale \(t_\mathrm{esc}\)), make sure to update the volume for the injection at each computation step!

Custom terms

As described above, AM\(^3\) solves a coupled set of transport equations with three energy-dependent terms per species/equation: (1) the injection term \(Q\), which can be set to arbitrary distributions as described above, (2) the escape term \(\alpha = 1/t_\mathrm{escape}\), which can be set to arbitrary distributions as described above, and (3) the continuous cooling term \(A=\dot{E}/E\) corresponding to advection in energy. In order to allow for full flexibility, we added the option to add an energy-dependent user specified \(A\) term as well, via an extra cooling timescale \(t_\mathrm{extra cooling} = 1/A\). It defaults internally to \(A=0\), corresponding to an infinitely large timescale, and has no switch to be turned on. It can be set via set_t_<species>_extra_cooling for species proton, neutron, pion, muon, neutrino, pair, photon, where pion, muon, neutrino, pair set the same timescale for all flavours/helicities/anti-particles. For an example see this notebook. We note that too spiky timescale might lead to numerical inaccuracies, and that strongly increasing timescales in \(A\) may challenge the solver, as it starts picking different solver regimes from the highest to lowest energies.

Performance profiline

To measure the run times for specific processes, set am3.set_profile_timing(1) (to retrieve the current setting, replace set by get). If profile timing is on, the run times for 36 sub steps of a single calculation timestep can be obtained through am3.get_profiled_times() and their labels through am3.get_profiled_time_labels().

This method can be called at any time during the simulation. It stores only the run times of the last time step.

An example of this procedure can be found in this example (see bottom plot).

Word of caution

Radiation energy density.: In principle, it is possible to reach a regime of catastrophically high radiation densities where the solver actually crashes due to numerical instability. The typical signature of this unstable regime is an exponential photon buildup followed by the replacement of densities and interaction rates with nans. This may be triggered by high injection luminosities of cosmic-ray electrons, cosmic-ray protons, or external photons, coupled with a strong magnetic field, or a low escape rate that leads to unphysical particle pileup over time. AM\(^3\) has no inbuilt sanity check for the selected simulation parameters, the injected particle rates, or the current radiation density levels. It is up to the user to verify at the end of the simulation that the densities and interaction rates are finite so as to confirm the validity of the predicted emission.
Solver time step.: As discussed under Simulation Parameters, the timestep of the numerical solver can be adjusted by the user for efficiency purposes. For example, in our AGN simulaton example, the system is optically thin and a timestep of 10% of the dynamical timescale was sufficient to capture even the fastest high-energy interactions. The numerical solver may not capture the process correctly if the time step selected is several orders of magnitude larger than the interaction timescales. Therefore, when simulating very optically thick environments, an adequately small time step has to be adopted by the user. AM\(^3\) does not internally check the relationship between the selected solver time step and the interaction timescales, nor the potential innaccuracies that can arise from setting time steps that are too large. It is up to the user to test the accuracy of the results, for example by checking that the predicted emission spectrum remains the same when a smaller solver time step is selected.
Flexibility of the model. One of the useful features of this release of AM\(^3\) is its generality regarding the geometry and nature of the environment being simulated. This empowers the user to build their own astrophysical model in a modular fashion; and with this great power comes also great responsibility. Users are free to apply AM\(^3\) to sources and environments that go beyond those tested so far, but they should be mindful of the underlying physics and critically evaluate the results produced. The AM\(^3\) team disclaims any responsibility for unintended or unphysical outcomes, including those resulting from oversights, errors, or misinterpretation in the application of the code. The user is responsible for comprehensively testing and validating the code when applying it to their specific purpose.

Deprecated functions

am3.set_escape_fraction_<particle>() and am3.get_escape_fraction_<particle>(), where <particle> can be photons, electrons, protons, neutrons, pions, muons, or neutrinos.
In conjunction with am3.set_escape_timescale, am3.set_escape_fraction_<particle>() allows the user to set an energy-independent escape rate for a given particle species. The absolute rate of physical escape is then given in \(s^{-1}\) by \(\eta_{X, esc}/ t_\mathrm{esc}\), which can be retrieved with am3.get_escape_fraction_<particle>() / am3.get_escape_timescale() .
The default value is 1.0 for all species. The set function can be called at any point during runtime.
Since version 1.2, it is recommended to use am3.set_t_<particle>_escape to directly set the escape rate in \(s^{-1}\). To set an energy-independent escape rate, a constant array can be used.\
am3.set_escape_fraction_charged_particles(fraction)and am3.set_escape_fraction_neutral_particles(fraction) can be used to set the same energy-independent escape fraction (see above) for all charged and neutral particles, respectively.
Since version 1.2, it is recommended to use am3..set_t_<particle>_escape for each particle species to set their escape rate in \(s^{-1}\). To set an energy-independent escape rate, a constant array can be used.\