© Y. Dallilar et al. 2022

1. Introduction

Synchrotron processes describe the energy loss of charged particles in the presence of a magnetic field. Although there are many applications of synchrotron processes in astrophysics, the most luminous synchrotron sources originate near extreme environments, such as black holes, and often in the form of jets. Measurement of the spectral energy distribution (SED) and its modelling allow extraction of the physical conditions in the vicinity of such sources.

A prime example is Sgr A^⋆, the nearest supermassive black hole located at a distance of 8.27 kpc towards the centre of our Galaxy (GRAVITY Collaboration 2019, 2021b; Do et al. 2019) with a mass of 4.3 × 10⁶ M_⊙. Recent interferometric observations with GRAVITY pinpoint the origin of Sgr A^⋆ flares to localised ‘hot-spots’ orbiting at several R_S (= 1.3 × 10¹² cm) in the accretion flow of the black hole, an interpretation also supported by polarisation signatures (GRAVITY Collaboration 2018, 2020a,b). However, possible emission mechanisms for Sgr A^⋆, or accreting black holes in general, are still subject to debate, particularly from near-infrared (NIR) to X-ray. Proposed scenarios include a mixture of synchrotron, synchrotron self-Compton (SSC), and inverse Compton (IC) models (Dodds-Eden et al. 2009; Eckart et al. 2012; Dibi et al. 2014; GRAVITY Collaboration 2021a). A variety of temporal models for the Sgr A^⋆ flares have also been developed with various approaches; see for example Yusef-Zadeh et al. (2008), Dodds-Eden et al. (2010), and Witzel et al. (2021).

Interpreting observations from state-of-the-art observational facilities necessitates flexible numerical utilities that can accommodate easy development of various source models to disentangle these emission scenarios and to explore the nature of these sources in detail. The Python package FLAREMODEL described in this paper is designed for this purpose¹. The code is developed around generic low-level utility functions which are written in C for performance. These functions constitute the core of FLAREMODEL. For convenience and accessibility, we developed the user interface in Python. The user interface includes exported utility functions and readily available source models as a guidance, some of which are described in this paper. Our code has also been used in GRAVITY Collaboration (2021a) to investigate emission mechanisms of a peculiar NIR/X-ray flare of Sgr A^⋆.

Thanks to the modular design of FLAREMODEL, it is also possible to build a variety of models applicable to astrophysical synchrotron sources with the help of these utility functions. This feature allows users to explore a wide selection of models in a short time and is one of the main strengths of our code in comparison to similar alternatives, such as AGNPY (Nigro et al. 2020). Options for spatial and temporal modelling are also available to a limited extent as described later in this paper.

Another distinguishing feature of this code is the support for parallel computations. In terms of overheads, OpenMP implementation in C is far more efficient than alternative implementations available in Python for parallel computations. Therefore, multi-threading is implemented at a lower level in C functions using OpenMP. There may be situations where higher level multi-threading is desired on the Python side but the efficiency is inherently hindered by the Python global interpreter lock (GIL). As the C functions in our code bypass GIL, this provides flexibility for developing source models targeting high-end computing hardware. As an experimental feature, we provide a Python interface to extend the functionality with algorithms targeting graphics processing units (GPUs) with OpenCL implementation. Currently, there is only built-in support for IC model calculations.

We further provide an SED fitting module built on top of the Python fitting package LMFIT (Newville et al. 2016). This package provides a user-friendly interface to a range of optimisation methods available within the Python ecosystem ranging from least-square minimisation to Markov-chain Monte Carlo methods. As FLAREMODEL opens many possibilities towards source modelling, we do not attempt to cover all the bases for fitting. Users can extend the fitting class following the examples in our documentation. Alternatively, users may opt for other fitting packages of their preference for which we do not provide explicit support and examples are not included in our documentation.

The organisation of this paper is as follows. In Sect. 2, we describe the structure of our code and physical processes implemented therein. Section 3 includes selected test problems and discusses selected cases for which we demonstrate this code.

2. Methodology

2.1. Electron distributions

Physical processes and interactions implemented in our code are limited to electrons. We provide a selection of built-in electron distributions with this code; thermal, kappa, and a family of power-law distributions. Power-law distributions are either in the form of a plain power law, a broken power law, or their exponential cutoff variants. The equations of the built-in electron distributions and their normalisation are detailed in Appendix A.

Our code includes a modular low-level C interface which is used to define built-in electron distributions. As utility functions in our code do not contain instructions specific to the electron distributions, the same interface can be used to add new electron distributions and these will be readily available for source models. Our code documentation covers the guidelines for defining new electron distributions in detail.

2.2. Synchrotron emissivity and absorption coefficients

We calculate the synchrotron emissivity and absorption coefficients (j_ν, α_ν) per steradian as follows,

$$\begin{array}{cc}\hfill j_{\nu }& =\frac{1}{4\pi }{\displaystyle {\int }_{{\gamma }_{min}}^{{\gamma }_{max}}{P}_{\mathrm{e}}{n}_{\mathrm{e}}(\gamma )\mathrm{d}\gamma ,}\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\alpha }_{\nu }& =\frac{1}{8\pi m_{\mathrm{e}}{\nu }^2}{\displaystyle {\int }_{{\gamma }_{min}}^{{\gamma }_{max}}{P}_{\mathrm{e}}(\frac{\mathrm{d}{n}_{\mathrm{e}}(\gamma )}{\mathrm{d}\gamma }-\frac{2n_{\mathrm{e}}(\gamma )}{\gamma })\mathrm{d}\gamma ,}\hfill \end{array}$$(2)

where m_e is the mass of an electron, n_e(γ) is the electron density, γ is the Lorentz factor, and ν is the frequency of the emitted/absorbed radiation. The power emitted by a single electron (P_e) is given by

$$\begin{array}{c}\hfill P_{\mathrm{e}}=\frac{\sqrt{3}{e}^{3}Bsin\varphi }{m_{\mathrm{e}}{c}^2}F(X).\end{array}$$(3)

Here, e is the charge of an electron, c is the speed of light, B is the magnetic field, and ϕ is the pitch angle. The F(X) function is defined as

$$\begin{array}{c}\hfill F(X)=X{\displaystyle {\int }_X^{\infty }{K}_{5/3}(x)\mathrm{d}x,}\end{array}$$(4)

where X = ν/ν_c and K(x)_5/3 is the modified Bessel function of order 5/3. The critical frequency (ν_c) is

$$\begin{array}{c}\hfill {\nu }_c=\frac{3}{4\pi }\frac{eB{\gamma }^{2}sin\varphi }{m_{\mathrm{e}}c}.\end{array}$$(5)

Evaluating the F(X) function is numerically expensive. Instead, we use pre-evaluated values of the function and perform interpolation at run-time in the range 10⁻⁴ < X < 10. Beyond these limits, we use the approximations from Rybicki (1979). The reader may refer to Blumenthal & Gould (1970) for the detailed formulation of synchrotron radiation.

There are two different methods implemented in our code; brute or userdist. The only difference between these two methods is the estimate of electron distribution derivative as required by Eq. (2). While the brute method uses the analytical form of a given electron distribution, the userdist method estimates the quantity from a user-defined numerical grid. Therefore, the brute method is more accurate and is the default method. The userdist method is originally motivated by non-trivial electron distributions arising from the temporal evolution of electrons, which we describe in Sect. 2.6.

Our default integration scheme relies on an equidistant grid in log(γ). Limits of the integration (γ_min, γ_max) and grid points per unit in log(γ) are user-defined. For some electron distributions, such as the thermal distribution, setting integration limits may be less intuitive. However, it is necessary within our formalism and implies truncated electron distributions at the limits. Computations involving power-law distributions adapt the integration limits from the distribution parameters (γ_min, γ_max) and the ones with exponential cutoff use the integration limits (γ_min, 10 × γ_max) instead. This only applies to the integration limits and does not necessarily reflect the normalisation of electron distributions, which are discussed in Appendix A.

For tangled magnetic fields, we approximate the angle-averaged computation with ϕ = arcsin(π/4). This is the preferred approach in some studies such as Dodds-Eden et al. (2010) for a faster computation by avoiding integration over the solid angle.

2.3. Inverse Compton scattering

Adapting the formalism of Blumenthal & Gould (1970), we define a Compton kernel (F_C),

$$\begin{array}{c}\hfill F_{\mathrm{C}}(ϵ,\gamma ,{ϵ}^{\prime })=\frac{3{\sigma }_{T}c}{4{\gamma }^2ϵ}[2qlnq+(1+2q)(1-q)+\frac{{({\mathrm{\Gamma }}_{\mathrm{e}}q)}^2(1-q)}{2(1+{\mathrm{\Gamma }}_{\mathrm{e}}q)}],\end{array}$$(6)

where Γ_e = 4ϵγ/m_ec² and q = ϵ′/Γ_e(γm_ec² − ϵ′), σ_T is the Thomson scattering cross-section. Here, ϵ and ϵ′ respectively refer to the incident and the scattering photon energies, that is, ϵ′=hν with h as the Planck constant. We then calculate the Compton emissivity (j_ν, IC) with a double integral along ϵ and γ,

$$\begin{array}{c}\hfill j_{\nu ,\mathrm{IC}}=\frac{h\nu }{4\pi }{\displaystyle {\int }_{ϵ}h{n}_{\mathrm{p}}(ϵ/h){\int }_{\gamma }{n}_{\mathrm{e}}(\gamma )F_{\mathrm{C}}(ϵ,\gamma ,{ϵ}^{\prime })\mathrm{d}\gamma \mathrm{d}ϵ,}\end{array}$$(7)

where n_p is the seed photon density.

The corresponding generic low-level utility function operates on a given arbitrary grid of seed photon density and electron distribution. While this is the most computationally expensive task in our code, it can be massively parallelised. Therefore, we also provide an OpenCL implementation using the PYOPENCL package², which executes the computations on a GPU and is hardware agnostic.

2.4. Ray-tracing

To simplify the radiative transfer calculations, either the slab approximation or a homogeneous sphere geometry is used for the analytical modelling of astrophysical synchrotron sources. Without these assumptions, radiative transfer calculations either become very complex or impossible. Ray-tracing is an alternative method for numerical calculations, for which imaginary rays passing through a source and exiting towards the line of sight are drawn. Specific intensities are calculated along these rays numerically. The specific intensities at the source surface are then used to estimate its brightness.

Our code provides a basic support for ray tracing to enable the modelling of synchrotron sources with complex geometries. The corresponding utility function is capable of parallel computation of multiple rays with the same length and fixed grid resolution, Δx. A single ray tracing has the form

$$\begin{array}{c}\hfill I_{\nu (k)}={\displaystyle \sum _{i=0}^k\mathrm{\Delta }x(j_{{\nu }_{(i)}}/{\alpha }_{{\nu }_{(i)}})\left(\prod _i^k(1-exp(-{\alpha }_{{\nu }_{(i)}}\mathrm{\Delta }x))\right),}\end{array}$$(8)

where I_ν(k) is the specific intensity at grid point k and frequency ν, and the boundary condition is set by I_ν(0) = 0. As such, this function stores the specific intensities along the rays although the stored data are not necessary to calculate synchrotron luminosity. Specific intensities on a grid can then be used to compute seed photon densities for SSC emission as described in Sect. 2.5.2.

2.5. Geometries

We provide two alternative geometries with our code. These demonstrate how to construct source models with the built-in utility functions.

2.5.1. Homogeneous sphere

A homogeneous sphere is the simplest assumption one can make on the geometry of compact synchrotron sources. It is therefore one of the most commonly used. As an advantage to analytical modelling, our implementation can compute synchrotron SEDs from a wide selection of electron distributions. Radiative transfer on a homogeneous sphere can be written as,

$$\begin{array}{c}\hfill L_{\nu }=4\pi {\displaystyle {\int }_0^R\frac{2\pi r{j}_{\nu }}{{\alpha }_{\nu }}[1-exp(-2{\alpha }_{\nu }\sqrt{R^2-r^2})]\mathrm{d}r,}\end{array}$$(9)

where R is the radius of the sphere³. We use the analytical form of this integral because its computation is faster. The luminosity integral reduces to

$$\begin{array}{c}\hfill L_{\nu ,\mathrm{thin}}=\frac{4\pi }{3}{R}^3(4\pi j_{\nu })\end{array}$$(10)

and to

$$\begin{array}{c}\hfill L_{\nu ,\mathrm{thick}}=4\pi \left(\frac{\pi j_{\nu }{R}^2}{{\alpha }_{\nu }}\right)\end{array}$$(11)

when the emission is optically thin and optically thick, respectively.

In this model, seed photon densities for SSC calculations are estimated using the equation

$$\begin{array}{c}\hfill n_{\mathrm{p}}(\nu )=\frac{L_{\nu }}{4\pi c{h}^2\nu R^2}.\end{array}$$(12)

We assume the medium is transparent to the IC emission. Hence, the IC luminosity (L_ν, IC) is equivalent to Eq. (10) replacing j_ν with j_ν, IC.

We include a Doppler beaming formalism shared by both geometrical models. If we define a Doppler factor,

$$\begin{array}{c}\hfill \delta \equiv {\mathrm{\Gamma }}^{-1}/(1-\beta sin\mathrm{\Phi })\end{array}$$(13)

with ${\mathrm{\Gamma }}^{-1}=\sqrt{1-{\beta }^2}$, for a source moving with velocity β(=v/c) at an angle Φ to the observer, the implementation is as follows,

$$\begin{array}{cc}\hfill L_{\mathrm{syn}}(\nu )& ={\delta }^n{L}_{\mathrm{syn}}^{\prime }({\nu }^{\prime }),\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill L_{\mathrm{SSC}}(\nu )& ={\delta }^n{L}_{\mathrm{SSC}}^{\prime }({\nu }^{\prime },n_{\mathrm{p}}^{\prime }),\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill L_{\mathrm{IC}}(\nu )& ={\delta }^n{L}_{\mathrm{IC}}^{\prime }({\nu }^{\prime },n_{\mathrm{p}}^{\star }),\hfill \end{array}$$(16)

where ν′=ν/δ and n is a geometry-dependent parameter; typically n = 3 for a point source. The quantities L′, ν′, and $n_{\text{p}}^{\prime }$ refer to the values in the rest frame of the source. The photon density of an external source ($n_{\text{p}}^{\star }$) needs to be supplied by the user because it is dependent on the modelling geometry. In our formulation, this quantity is in the rest frame and at the location of the main source. This is implemented in the form of a function decorator in Python, and therefore the definition can be adapted to different scenarios; for example, extra-galactic sources with a redshift.

2.5.2. Radial sphere

The treatment of inhomogeneous spherical sources is discussed in Gould (1979) and Marscher (1977), and a full theoretical framework is introduced in Band & Grindlay (1985) for power-law spheres, in which electron density and magnetic field radially scale in the form of a power law. Our implementation of a radial sphere is similar to those concepts but generalised to a spherical source with any given electron distribution, and it is possible to define an arbitrary radial scaling of electron density and magnetic field.

To estimate synchrotron emission from the sphere, we first calculate j_ν(r) and a_ν(r) on a one-dimensional regular grid where the number of grid points is specified by the user n_r = R/Δr. This is only valid for R_in < r < R where R_in and R are respectively the inner and the outer radii of the source. j_ν(r) and a_ν(r) are set to zero beyond these limits. We then interpolate j_ν and a_ν to a regular Cartesian grid of size (R, 2R) and of resolution Δr. In the end, we perform ray tracing with the utility function provided in our code; see Eq. (8). From here on, computing the synchrotron luminosity is straightforward with

$$\begin{array}{c}\hfill L_{\nu }=4\pi {\displaystyle \sum }{I}_{\nu }(r)2\pi r\mathrm{\Delta }r.\end{array}$$(17)

Unlike the formalism used in the homogeneous sphere model (Eq. (12)), this geometrical model provides proper numerical estimates of the seed photon densities within a spherical source. Even for a homogeneous sphere (physically, not referring to the corresponding model in Sect. 2.5.1), the formalism is not exact and seed photon densities additionally vary with radius. Here, we describe our methodology while the numerical result is explored in Sect. 3.2.

Our methodology is similar to the concept of Band & Grindlay (1985), in which they apply a ray-tracing solution for rays passing through a single point in radius from all directions. The specific intensities of each ray at that point are then normalised. In our case, we only apply ray tracing along a single direction as described above, and then we normalise specific intensities at grid points falling into a layer of finite size in radius.

Mathematically, this computation can be summarised as

$$\begin{array}{c}\hfill n_{\mathrm{p}}(\nu ,r)=\left(\frac{2\pi }{h^2\nu c}\right)\frac{{\sum }^{\mathcal{R}(r)}{I}_{{\nu }_{(i,j)}}sin{\theta }_{(i,j)}}{{\sum }^{\mathcal{R}(r)}sin{\theta }_{(i,j)}},\end{array}$$(18)

where ℛ(r) represents the grid points in the range r − Δr/2 < r_(i, j) < r + Δr/2. Finally, we sum the contributions of each layer to the total SSC emission using the computed n_p(ν, r)⁴.

2.6. Temporal evolution

We include a basic temporal model for high-energy electrons used in Dodds-Eden et al. (2010) to investigate the temporal evolution of a bright and complex multi-wavelength flare of Sgr A^⋆. With this formulation, it is possible to evolve and calculate custom electron distributions and/or to build time-dependent models simulating adiabatic/synchrotron cooling with particle injection, by writing

$$\begin{array}{c}\hfill \frac{\partial N_{\mathrm{e}}(\gamma ,t)}{\partial t}=c_{\mathrm{inj}}(t)N_{\mathrm{e}}^{\star }(\gamma )-\frac{N_{\mathrm{e}}(\gamma ,t)}{{\tau }_{\mathrm{esc}}}+\frac{\partial \gamma N_{\mathrm{e}}(\gamma ,t)}{\partial \gamma }.\end{array}$$(19)

The first term in the equation represents the particle injection into the system with normalised electron distribution $N_{\text{e}}^{\star }(\gamma )$ and of rate c_inj(t). The particle escape is given in the second term with a timescale of τ_esc. The third term simulates the cooling of electrons in the system, where ̇γ is the sum of contributions from adiabatic and synchrotron cooling. Instantaneous adiabatic losses when geometry is reduced to a sphere is given by

$$\begin{array}{c}\hfill -{\gamma }_{\text{adiab}}=\gamma \frac{\mathrm{d}ln(R(t))}{\mathrm{d}t}=\frac{\gamma c{\beta }_{exp}}{R(t)},\end{array}$$(20)

where β_exp is the expansion speed of the source in units of c, and R is the instantaneous radius of the source at time t. Angle-averaged synchrotron losses can be written as

$$\begin{array}{c}\hfill -{\gamma }_{\mathrm{syn}}={\gamma }^2\left(\frac{4{\sigma }_T}{3m_{\mathrm{e}}c}\right)\left(\frac{B^2}{8\pi }\right)=\frac{{\gamma }^2{B}^2}{7.74\times {10}^8}{\mathrm{s}}^{-1}.\end{array}$$(21)

As the temporal evolution of electrons can be formulated in many different ways, our utility function only performs one step integration from one state to the following state after a time-step of Δt. Appendix D details our numerical formulation. We provide Python interfaces for selected scenarios as guidance, such as the examples in Sects. 3.3 and 3.4.

The effects of IC processes are not included in this model. Therefore, there is a hidden assumption that IC cooling is negligible compared to other cooling processes. With a first-order approximation, the relative contribution of IC to synchrotron cooling can be estimated from the ratio of the seed photon energy density (u_p) to the magnetic energy density (u_B = B²/8π). This approach can serve as a validity check for the compatibility of a given scenario with this model.

In its most basic form, IC cooling can be similarly implemented by introducing an effective energy density (u_eff = u_p + u_B) to Eq. (21) at the cost of an additional computation time associated with the estimate of u_p, which needs to be updated at each integration step. However, there are additional concerns associated with it such as radiation pressure and multiple scatterings when IC processes are dominant. Implementing these physics would require much more complex formulations, which are beyond the scope of this study. If they are not implemented, there is no clear cut to when these physics will be problematic. Hence, we choose not to explicitly support IC cooling in a generic method with no restrictions.

3. Tests and examples

3.1. Synchrotron emission from a power-law sphere

To demonstrate our radial sphere implementation, we reproduce the solution of Band & Grindlay (1985) (in their Fig. 1b) corresponding to their power-law sphere formulation. Using the definition provided by these latter authors, the magnetic field in the sphere scales with B = B₀(r/R_in)^−m and m = 1; the electron density scales with n_e = n_e, 0(r/R_in)⁻ⁿ and n = 2. The index of the power-law electron distribution is p = 3 (p = 2α + 1). Finally, the radius is parameterised as x = R/R_in with x = 50. The choice of B₀ and n_e, 0 is not relevant for the solution as long as the effects of (γ_min, γ_max) are negligible in a selected frequency range.

The resulting SED shown in Fig. 1 has three distinct regimes from higher to lower frequencies: (1) The sphere is optically thin above ν_T, which we calculate numerically, and is the transition frequency to the optically thin emission. (2) The partially opaque regime can be interpreted as a superposition of SEDs corresponding to different layers along the radius. (3) The sphere is fully opaque below a certain frequency (ν ≲ 3ν_T × 10⁻²) and produces an SED with L_ν ∝ ν^2.5.

Fig. 1. Synchrotron SED using our radial sphere implementation constructed with parameters m = 1, n = 2, x = 50 and α = 1 as formulated for a power-law sphere in Band & Grindlay (1985). SEDs plotted in different colours refer to calculations at selected radial resolutions (R/ΔR) – 50, 100, and 200 grid points in radius respectively for x, 2x, and 4x. The black dashed line is the analytical form of the expected optically thin emission in Eq. (22). νT is the numerical estimate of the frequency where the SED peaks.

There is no exact analytical solution for this particular problem⁵. Band & Grindlay (1985) provide an analytical form describing the optically thin regime as

$$\begin{array}{c}\hfill L_{\nu }=\frac{4\pi R_{\mathrm{in}}^3}{3}\left(\frac{x^{3-s}-1}{3-s}\right)[4\pi j_{\nu }(R_{\mathrm{in}})],\end{array}$$(22)

where s = n + m(a + 1). This is the dashed line in Fig. 1.

One reason to present this particular example is to further discuss our ray-tracing solution. We provide our solution computed at different radial resolutions with multiples of x in the figure. While all three computations agree below ν_T, the reproducibility of the optically thin SED requires common sense as to the choice of radial resolution, because it is necessary to resolve the innermost region near R_in as shown in Eq. (22). The difference can be more dramatic with steeper slopes either in the magnetic field or in the electron density along the radius and a finer resolution may be required.

3.2. Synchrotron self-Compton emission from a homogeneous sphere

One of the challenging steps for accurate SSC calculations is the estimation of seed photon densities (from synchrotron emission) within a given source. A common approach is to use Eq. (12) as in our homogeneous sphere model for simplicity, but this is only valid at a radius beyond the surface of the source. Proper estimation of seed photon densities at a given location (r) requires summation of contributions from different locations (r ′) in an entire source. Hence, an exact analytical solution is not possible and numerical calculations are computationally expensive.

Briefly, we only discuss our results for a homogeneous sphere complementary to the description in Sect. 2.5.2. There are two limits to the photon densities:

• At frequencies where the sphere is opaque, the mean free path of photons is small compared to the size of the source. Hence, photon densities in this limit are constant except near the outer regions of the sphere.

• At frequencies where the sphere is optically thin, the radial profile is given by

  $$\begin{array}{c}\hfill \mathrm{\Phi }(x)=\frac{2}{3}(\frac{1-x^2}{2x}ln\left|\frac{1+x}{1-x}\right|+1),\end{array}$$(23)

  where x is the normalised radius r/R (Band & Grindlay 1985).

We show our results for selected frequencies in Fig. 2. Our calculations compare well with the expected radial profiles in both limits. This result alone may not be very intuitive. Therefore, we also compute the SSC spectrum for p of 2, 2.5, and 3 in Fig. 3. Compared to the approximation with Eq. (12), our code predicts 15%–20% excess SSC emission in the optically thin regime and about 50% excess in the optically thick regime with small variations at different p values. The difference can be incorporated in n_e and γ_min in this specific case. As a result, for a homogeneous sphere, this calculation is unnecessarily expensive: it roughly scales linearly with the radial resolution.

Fig. 2. Top: representative synchrotron SED where vertical dashed lines are equivalent to the colour-coded frequencies in the bottom plot. Bottom: numerical solution to the photon densities in a homogeneous sphere at the selected frequencies. The optically thin and thick limits are shown with the two dot-dashed lines.

Fig. 3. Synchrotron self-Compton SED for a homogeneous sphere with p values 2, 2.5, and 3 shown in different colours. Solid lines represent the SEDs with the traditional photon density estimates and dashed lines those with the numerical photon densities as described in the text. Bottom plot: difference between the two methods.

3.3. Adiabatic cooling

The adiabatic cooling model introduced in van der Laan (1966) is still one of the most commonly used models to describe the temporal evolution of compact synchrotron sources at longer wavelengths (typically from radio to sub-mm) due to its simplicity. In essence, the model is a recipe to produce synchrotron flares without injection of electrons to the system. As the source expands, the self-absorbed peak of the synchrotron emission propagates towards lower frequencies with the combination of decreasing electron density (∝1/R³) and magnetic field (∝1/R²). This evolution describes the rise of the flare. As the self-absorbed peak propagates through a given frequency, flux starts to drop due to adiabatic cooling, and therefore the decay of the flare sets in. Naturally, this formalism predicts a certain time-delay from shorter to longer wavelengths.

Aside from its simplicity and fast numerical evaluation, there are some inherent drawbacks to the model. The model can only describe the temporal evolution of the synchrotron SED arising from a plain power-law electron distribution and it is not possible to directly infer the physical properties of the source. In addition, physical implications of the magnetic field, such as synchrotron cooling, are not taken into account, whereas in our numerical implementation, the temporal model directly operates on the electron population itself and can simulate the adiabatic and synchrotron cooling simultaneously.

We set up a test problem using the parameters from Yusef-Zadeh et al. (2008) as a reference (p = 2, B₀ = 10 G, β_exp = 0.02c, R₀ = 2R_S of Sgr A^⋆ with R_S = 1.3 × 10¹² cm). We set γ_max of the initial electron distribution to a large number so it is not a factor in our calculations. We integrate over Eq. (19) with a time resolution of one second. It is not necessary to develop the full model here. We rather use the scaling relations from van der Laan (1966),

$$\begin{array}{cc}\hfill L_{\nu ,m}(R)& =L_{\nu ,m}(R_0){(R/R_0)}^{-(4p+6)/(p+4)},\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\nu }_m(R)& ={\nu }_m(R_0){(R/R_0)}^{-(7p+3)/(p+4)},\hfill \end{array}$$(25)

where L_ν, m(R) and ν_m(R) represent the self-absorbed peak of SED at a given radius. In Fig. 4, we plot results of our numerical integration at selected radii with synchrotron cooling enabled and disabled. We also overplot the evolution of the self-absorbed peak predicted from the van der Laan (1966) model. As seen in the figure, synchrotron cooling introduces a spectral cutoff at high energies which propagates to lower frequencies as the source expands. This is despite the choice of large γ_max, which has no effect here.

Fig. 4. Temporal evolution of synchrotron SEDs calculated at selected radii with the utility function discussed in Sect. 2.6. We refer to the text for further details. Coloured dashed lines refer to an evolution in which synchrotron cooling is assumed to be negligible, and this would be equivalent to the traditional van der Laan (1966) formalism. Solid lines represent the case in which synchrotron cooling is enabled. The black line is the expected evolution of the self-absorbed peak according to the van der Laan (1966) model.

An analytical approximation for the evolution of the high-energy cutoff is provided in Appendix C. The final result is sensitive to the initial magnetic field and implies an upper limit to the magnetic field where the van der Laan (1966) assumption is still valid. In other words, the underlying magnetic field needs to be large enough to produce the observed flux arising from synchrotron processes, but at the same time it needs to be small enough so that the impact of the synchrotron cooling to the system is negligible at the frequency range of interest. This aspect of the model limits its applicability to lower frequencies, typically radio to sub-mm range, where the effects of synchrotron cooling are less prominent.

3.4. Synchrotron cooling

There are two traditional textbook solutions to the synchrotron cooling: (1) A propagating break (γ_b) in the electron distribution during sustained particle injection, and (2) a propagating cutoff (γ_max) during limited or no injection (Longair 2011). In both cases, the formulation is similar (either for the break or the cutoff in the electron distribution respectively) with

$$\begin{array}{c}\hfill {\gamma }_{b,max}\sim 1.3\times {10}^5{(B/10\mathrm{G})}^{-2}{(t_{\mathrm{syn}}/1\min )}^{-1},\end{array}$$(26)

where t_syn is comparable to the injection timescale (t_inj, or other dynamic timescale in the system) for the former case and the time after the particle injection for the latter.

We proceed with a formulation somewhere between the two extremes to demonstrate both cases at once. We consider particle injection in the form of a power-law distribution with p = 2, γ_min = 50 and γ_max = 4 × 10⁵. The particle injection rate has a Gaussian profile with σ = 5 min, peaking at 3σ. In this setup, we can assume an effective injection timescale of t_inj ∼ 2.35σ/2. We fix the magnetic field to 20 G.

We plot the evolution of the electron distribution with this experimental setup in Fig. 5. While it is difficult to describe each time stamp with a standard electron distribution, we pick two cases where an approximation is possible. These refer to the time stamps; injection peak at 15 minutes and when synchrotron cooling ‘erases’ the history of the particle injection at t ∼ 30 min.

Fig. 5. Evolution of a power-law electron distribution under synchrotron cooling and Gaussian particle injection with parameters noted in the text. Time evolution of the electron distribution is presented as snapshots taken every 5 min. Grey dashed lines refer to the rising part of the evolution and solid lines to the decaying part. The red line represents the spectrum at the peak of injection profile (t = 15 min) whereas red dots are the approximation with a broken power-law distribution. The blue line is the electron distribution 15 min after injection peak (t = 30 min) and blue dots are the approximation with a simple power law.

The former case can be approximated with a broken power-law electron distribution. The break in the electron distribution is estimated from t_inj and γ_b using the arguments described above. Additionally, we use an exponential cutoff in the electron distribution setting γ_max = (γ_max)_inj/1.5. The latter case can be approximated with a simple power-law distribution. In this case, γ_max is chosen so that there is synchrotron cooling with an effective duration of 15 min, that is, the 15 min following the peak injection.

As demonstrated here, analytical approximations are possible for this particular setup by reducing the problem to selected time episodes. Similar arguments have been used either to describe flare SEDs linking certain source properties to spectral features (Ponti et al. 2017; GRAVITY Collaboration 2021a) or to extract physical properties of sources at certain time intervals (Dallilar et al. 2017). However, there are also studies providing evidence for multiple injection and cooling episodes or perturbations in physical properties during flaring episodes. For example, Dodds-Eden et al. (2010) investigated a bright Sgr A^⋆ flare with a rapid substructure and linked changes in the magnetic field to the substructure in the multi-wavelength light curves. Numerical tools such as those described here provide the potential to investigate the dynamical evolution of such events in more detail, which is useful for disentangling possible emission scenarios and understanding the history of particle acceleration and/or cooling.

4. Conclusion

In this paper, we present FLAREMODEL, an open-source Python package for one-zone modelling of astrophysical synchrotron sources. We discuss our methodology and the structure of our code along with selected examples to demonstrate the capabilities and options for future modelling efforts.

Instead of focusing on selected source models, we provide a set of generic low-level utility functions from which users can build their own models. This is a step towards modularity and ease of development in future modelling efforts while still keeping the duration of numerical calculations reasonably low with multi-threading support. As such, basic models can still be executed on low-end hardware but complex models designed for high-end hardware can take advantage of built-in multi-threading support. This approach distinguishes FLAREMODEL from similar alternatives.

We do not provide detailed benchmarks with this paper but only order-of-magnitude estimates for built-in geometrical models. We use the default grid settings. SEDs in each case are computed for 100 points in frequency on a workstation with an i7-8700 CPU using a single thread. For the homogeneous sphere model, SEDs are computed within 1 ms and 100 ms for synchrotron and SSC, respectively. For the radial sphere model, computations are performed in roughly 100 ms and several seconds for synchrotron and SSC, respectively. When SSC calculations are offloaded to the integrated graphics of this CPU, the execution time reduces to about 1 s.

Acknowledgments

SvF, and FW acknowledge support by the Max Planck International Research School. GP acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 865637).

References

Appendix A: Analytical formulation of implemented electron distributions

As described in the text, we provide a handful of useful electron distributions with our code and their analytical formulation is given here. The normalised thermal electron distribution is implemented as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\mathrm{th}}}{\mathrm{d}\gamma }=\frac{\gamma {({\gamma }^2-1)}^{1/2}}{{\mathrm{\Theta }}_e{K}_2(1/{\mathrm{\Theta }}_e)}exp(-\frac{\gamma }{{\mathrm{\Theta }}_e}),\end{array}$$(A.1)

where Θ_e(=k_BT/m_ec) is the dimensionless electron temperature and K₂ is the modified Bessel function of the second kind.

For non-thermal electron distributions, there are two main options. For a power-law electron distribution, we use the form

$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\mathrm{pl}}}{\mathrm{d}\gamma }=N_{\mathrm{pl}}{\gamma }^{-p},\text{for}{\gamma }_{min}\le \gamma \le {\gamma }_{max},\end{array}$$(A.2)

where p is the exponent of the distribution, and γ_min and γ_max are lower and higher energy limits, respectively. We define the broken power-law distribution as the combination of two power-law distributions with exponents p₁ and p₂, below and above the break in gamma (γ_b), respectively. This latter has the form,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\text{bpl}}}{\mathrm{d}\gamma }=N_{\mathrm{bpl}}\times \{\begin{array}{cc}{\gamma }^{-p_1},\hfill & \text{if}{\gamma }_{min}\le \gamma \le {\gamma }_b,\hfill \\ {\gamma }^{-p_2}{\gamma }_b^{p_2-p_1},\hfill & \text{if}{\gamma }_b<\gamma \le {\gamma }_{max}.\hfill \end{array}\end{array}$$(A.3)

The normalisation factors for these distributions are respectively

$$\begin{array}{c}\hfill N_{\mathrm{pl}}=\frac{p-1}{{\gamma }_{min}^{1-p}-{\gamma }_{max}^{1-p}},\end{array}$$(A.4)

and

$$\begin{array}{c}\hfill N_{\mathrm{bpl}}={[\left(\frac{{\gamma }_{min}^{1-p_1}-{\gamma }_b^{1-p_1}}{p_1-1}\right)+{\gamma }_b^{p_2-p_1}\left(\frac{{\gamma }_b^{1-p_2}-{\gamma }_{max}^{1-p_2}}{p_2-1}\right)]}^{-1}.\end{array}$$(A.5)

In addition, we provide an alternative form for both power-law electron distributions with an exponential cutoff at γ_max. Explicitly, this is implemented as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\text{expc}}}{\mathrm{d}\gamma }={\left(\frac{\mathit{dn}}{d\gamma }\right)}_{\mathrm{pl},\mathrm{bpl}}exp(-\frac{\gamma }{{\gamma }_{max}}),\text{for}{\gamma }_{min}\le \gamma .\end{array}$$(A.6)

In this case, we still use the corresponding normalisation factors of each distribution given in Eqs. (A.4) and (A.5). Although this formulation removes the necessity of an upper limit on γ, the synchrotron routines enforce an upper limit to integration at 10 × γ_max for practicality.

Lastly, we include the kappa electron distribution as described in Pandya et al. (2016). However, the normalisation of this electron distribution is not exact in their formulation. The error in selected cases is provided in their paper. We do not go into further details here.

Appendix B: Synchrotron emissivity and absorption coefficients

In Fig. B.1, we compare our numerical computations of synchrotron emissivity and absorption coefficients with the publicly available SYMPHONY (Pandya et al. 2016) code for thermal (Θ_E = 30) and power-law (p = 3, γ_min = 1, γ_max = 10³) distributions with B = 10 G and ϕ = π/4. Numerical integration over the electron distributions with discontinuities is problematic and is also the case with SYMPHONY. We can only avoid this complication by setting proper integration ranges with our brute method.

Fig. B.1. Numerical computation of synchrotron emissivity (black) and absorption coefficient (red) with the brute and userdist methods in our code. Thermal (top) and power-law (bottom) electron distributions are shown separately with the parameters specified in the text. The results with SYMPHONY are overplotted.

Hence, we use the convenient exponential cutoff option at γ_max in SYMPHONY and equivalent implementation in our code. This also motivates setting ${\gamma }_{min}(\sim \sqrt{\nu /{\nu }_c})$ low enough so that it is not a factor in this comparison. Our computations are well within a few percent in the range of the plots with the exception of the userdist method and then only for ν/ν_c < 10⁶. This originates from the estimate of dn(γ)/dγ in Eq. (2). This calculation is always more accurate using an analytical form of the distribution. Therefore, we promote the brute method as the default while the userdist method remains as an option for the electron distributions with no analytical form.

Appendix C: Synchrotron cooling in the presence of adiabatic expansion

When Eq. (21) is integrated under the assumption of a spherical expansion at a constant speed, i.e. R(t) = R₀ + βct, the result is

$$\begin{array}{c}\hfill {\gamma }_{max}(t)=\frac{3\times 7.74\times {10}^8}{{B_0}^2{t}_{exp}}{(1-{\left(\frac{R_0}{R(t)}\right)}^3)}^{-1}\sim \frac{2.3\times {10}^9}{{B_0}^2{t}_{exp}},\end{array}$$(C.1)

where t_exp = R₀/βc and γ_max approaches its asymptotic value for 1 < < R₀/R. For the choice of parameters in Sect. 3.3, t_exp is about 4.2 × 10³ seconds. With the critical frequency defined in Eq. (5), we can estimate the time evolution of ν_max as

$$\begin{array}{c}\hfill {\nu }_{max}(t)\sim 0.3\times {\nu }_c({\gamma }_{max},B,\pi /4)\sim 3\times {10}^{14}{\left(\frac{B_0}{10\text{G}}\right)}^{-3}{\left(\frac{R_0}{R(t)}\right)}^2\text{Hz.}\end{array}$$(C.2)

There are two main points to emphasise as a result of this formulation. Synchrotron cooling is effective only during the early stages of the expansion with an effective duration on the order of t_exp. This determines the high energy cutoff in the electron distribution. At the same time, magnetic flux conservation implies a rapidly decaying magnetic field which pushes ν_max originating from synchrotron cooling to even lower frequencies as the source expands.

Appendix D: Numerical scheme used in the temporal evolution

The temporal evolution model introduced in Sect. 2.6 is a two-dimensional problem. We solve the equation in two steps. Firstly, we estimate the quantities on the right-hand side of the equation in γ grid at time t. We then calculate electron distribution at t + Δt with the forward Euler method.

Going in the reverse order, the time integration is discretised in the form,

$$\begin{array}{c}\hfill N_{e,\mathrm{i}}^{+}=N_{e,\mathrm{i}}^{-}+(A_{\mathrm{i}}^{-}+B_{\mathrm{i}}^{-})\mathrm{\Delta }t,\end{array}$$(D.1)

where i denotes the grid points in γ. We use the notation of minus (−) and plus (+) respectively for quantities at t and t + Δt. The quantities $A_{\text{i}}^{-}$ and $B_{\text{i}}^{-}$ are separately written in an open form below,

$$\begin{array}{cc}\hfill A_{\mathrm{i}}^{-}& =c_{\mathrm{inj}}{N}_{e,\mathrm{i}}^{\star }-\frac{N_{e,\mathrm{i}}}{{\tau }_{\mathrm{esc}}},\hfill \end{array}$$(D.2)

$$\begin{array}{cc}\hfill B_{\mathrm{i}}^{-}& =\frac{{\gamma }_{\mathrm{i}+1}{N}_{e,\mathrm{i}+1}-{\gamma }_{\mathrm{i}}{N}_{e,\mathrm{i}}}{{\gamma }_{\mathrm{i}+1}-{\gamma }_{\mathrm{i}}}.\hfill \end{array}$$(D.3)

Here, we drop the plus and minus notation for readability, because all the quantities are at time t. The term $A_{\text{i}}^{-}$ is reserved for particle injection and escape, and $B_{\text{i}}^{-}$ for adiabatic and synchrotron cooling. Calculating ̇γ is straightforward with the sum of Eqs. (20) and (21).

In the end, with the assumption of a spherical source, the corresponding utility function updates the radius (R) and the magnetic field (B) in-place with the equations,

$$\begin{array}{cc}\hfill R^{+}& =R^{-}+c{\beta }_{exp}\mathrm{\Delta }t,\hfill \end{array}$$(D.4)

$$\begin{array}{cc}\hfill B^{+}& =B^{-}{(R^{-}/R^{+})}^2.\hfill \end{array}$$(D.5)

As a word of caution, Δt is the time-step of the integration, not simply any time interval. The utility function provided with our code can only perform one step integration as described in this section. The proper choice of Δt and the formulation of the temporal evolution remains at the discretion of the user.

All Figures

Fig. 1. Synchrotron SED using our radial sphere implementation constructed with parameters m = 1, n = 2, x = 50 and α = 1 as formulated for a power-law sphere in Band & Grindlay (1985). SEDs plotted in different colours refer to calculations at selected radial resolutions (R/ΔR) – 50, 100, and 200 grid points in radius respectively for x, 2x, and 4x. The black dashed line is the analytical form of the expected optically thin emission in Eq. (22). νT is the numerical estimate of the frequency where the SED peaks.

In the text

	Fig. 2. Top: representative synchrotron SED where vertical dashed lines are equivalent to the colour-coded frequencies in the bottom plot. Bottom: numerical solution to the photon densities in a homogeneous sphere at the selected frequencies. The optically thin and thick limits are shown with the two dot-dashed lines.
In the text

	Fig. 3. Synchrotron self-Compton SED for a homogeneous sphere with p values 2, 2.5, and 3 shown in different colours. Solid lines represent the SEDs with the traditional photon density estimates and dashed lines those with the numerical photon densities as described in the text. Bottom plot: difference between the two methods.
In the text

Fig. 4. Temporal evolution of synchrotron SEDs calculated at selected radii with the utility function discussed in Sect. 2.6. We refer to the text for further details. Coloured dashed lines refer to an evolution in which synchrotron cooling is assumed to be negligible, and this would be equivalent to the traditional van der Laan (1966) formalism. Solid lines represent the case in which synchrotron cooling is enabled. The black line is the expected evolution of the self-absorbed peak according to the van der Laan (1966) model.

In the text

Fig. 5. Evolution of a power-law electron distribution under synchrotron cooling and Gaussian particle injection with parameters noted in the text. Time evolution of the electron distribution is presented as snapshots taken every 5 min. Grey dashed lines refer to the rising part of the evolution and solid lines to the decaying part. The red line represents the spectrum at the peak of injection profile (t = 15 min) whereas red dots are the approximation with a broken power-law distribution. The blue line is the electron distribution 15 min after injection peak (t = 30 min) and blue dots are the approximation with a simple power law.

In the text

	Fig. B.1. Numerical computation of synchrotron emissivity (black) and absorption coefficient (red) with the brute and userdist methods in our code. Thermal (top) and power-law (bottom) electron distributions are shown separately with the parameters specified in the text. The results with SYMPHONY are overplotted.
In the text