Thyr: A Volumetric Ray-Marching Tool for Simulating Microwave Emission

TL;DR

Thyr is an open-source, modular tool that uses a novel ray-marching technique to simulate microwave emission from solar flares with high spatial resolution, aiding interpretation of observations and including chromospheric effects.

Contribution

The paper introduces Thyr, a new volumetric ray-marching tool for detailed 3D microwave emission modeling of solar flares, including chromospheric thermal emission and high-resolution effects.

Findings

High-resolution simulation reveals spectral hardening at higher frequencies.

Hardening originates from magnetic convergence at the chromosphere top.

Model can be used for coordinated flare observation analysis.

Abstract

Gyrosynchrotron radiation is produced by solar flares, and can be used to infer properties of the accelerated electrons and magnetic field of the flaring region. This microwave emission is highly dependent on many local plasma parameters, and the viewing angle. To correctly interpret observations, detailed simulations of the emission are required. Additionally, gyrosynchrotron emission from the chromosphere has been largely ignored in modelling efforts, and recent studies have shown the importance of thermal emission at millimetric wavelengths. Thyr is a new tool for modelling microwave emission from three-dimensional flaring loops with spatially varying atmosphere and increased resolution in the lower corona and chromosphere. Thyr is modular and open-source, consisting of separate components to compute the thermal and non-thermal microwave emission coefficients and perform three-dimensional radiative transfer (in local thermodynamic equilibrium). The radiative transfer integral is computed by a novel ray-marching technique to efficiently compute the contribution of many volume elements. This technique can also be employed on a variety of astrophysics problems. Herein we present a review of the theory of gyrosynchrotron radiation, and two simulations of identical flare loops in low- and high-resolution performed with Thyr, with a spectral imaging analysis of differing regions. The high-resolution simulation presents a spectral hardening at higher frequencies. This hardening originates around the top of the chromosphere due to the strong convergence of the magnetic field, and is not present in previous models due to insufficient resolution. This hardening could be observed with a coordinated flare observation from active radio observatories.

Equations

$$\begin{array}{rl}& n_{\pm }(\nu ,\theta {)}^2=1-\\ & \frac{2X(1-X)}{2(1-X)-Y^2{\sin }^2\theta \pm [Y^4{\sin }^4\theta +4Y^2(1-X{)}^2{\cos }^2\theta {]}^{1/2}}\end{array}$$

$$\begin{array}{rl}X& :=\frac{{\nu }_p^2}{{\nu }^2},\\ Y& :=\frac{{\nu }_B}{\nu },\\ {\nu }_p& :=\text{electron plasma frequency}=e\sqrt{\frac{n_p}{\pi m_e}},\\ {\nu }_B& :=\text{electron gyrofrequency}=\frac{eB}{2\pi m_{e}c},\\ e& :=\text{electron charge,}\\ m_e& :=\text{electron mass,}\end{array}$$

$$a_{\theta \pm }(\nu ,\theta )=\frac{2(1-X)\cos \theta }{-Y{\sin }^2\theta \pm [Y^2{\sin }^4\theta +4(1-X{)}^2{\cos }^2\theta {]}^{1/2}}$$

$$\begin{array}{rl}{ϵ}_{\pm }(\nu ,\theta ,\gamma ,\varphi )& =\frac{2\pi e^2}{c}{\nu }^2\frac{n_{\pm }^2}{1+a_{\theta \pm }^2}\\ & \cdot \sum _{s=-\infty }^{\infty }{\left[-\beta \sin \varphi J_s^{\prime }(x_s)+a_{\theta \pm }\frac{\cos \theta -n_{\pm }\beta \cos \varphi }{n_{\pm }\sin \theta }{J}_s(x_s)\right]}^2\\ & \cdot \delta [(1-n_{\pm }\beta \cos \varphi \cos \theta )\nu -s{\nu }_B/\gamma ]\end{array}$$

$$\begin{array}{rl}c& :=\text{speed of light,}\\ \gamma & :=(1-{\beta }^2{)}^{-1/2},\\ \beta & :=\frac{v}{c},\\ x_s& :=\frac{\gamma \nu }{{\nu }_B}{n}_{\pm }\beta \sin \varphi \sin \theta ,\\ J_s& :=\text{Bessel function of the first kind, of order s,}\\ & \text{and }{J}_s^{\prime }\text{ its derivative with respect to }{x}_s,\\ \delta & :=\text{Dirac delta function}.\end{array}$$

$$2\pi {\int }_1^{\infty }d\gamma {\int }_{-1}^{1}d(\cos \varphi )g(\gamma ,\varphi )=1.$$

$$j_{\pm }(\nu ,\theta )=2\pi n_e{\int }_1^{\infty }d\gamma {\int }_{-1}^{1}d(\cos \varphi )g(\gamma ,\varphi ){ϵ}_{\pm }(\nu ,\theta ,\gamma ,\varphi ),$$

$$\begin{array}{r}{\kappa }_{\pm }(\nu ,\theta )=\frac{2\pi n_e}{m_e{\nu }^2{n}_{\pm }}{\int }_1^{\infty }d\gamma {\int }_{-1}^{1}d(\cos \varphi ){ϵ}_{\pm }(\nu ,\theta ,\gamma ,\varphi )h_{\pm }(\nu ,\theta ,\gamma ,\varphi )\end{array}$$

$$h_{\pm }(\nu ,\theta ,\gamma ,\varphi )=\left[-\beta {\gamma }^2\frac{\partial }{\partial \gamma }\left(\frac{g(\gamma ,\varphi )}{\beta {\gamma }^2}\right)+\frac{n_{\pm }\beta \cos \theta -\cos \varphi }{{\beta }^2\gamma \sin \varphi }\frac{\partial }{\partial \varphi }g(\gamma ,\varphi )\right].$$

$$\begin{array}{rl}{G}_{\pm }(\nu ,\theta ,\gamma ,s)=& {\left[-\beta \sin {\varphi }_s{J}_s^{\prime }(x_s)+a_{\theta \pm }\left(\frac{\cos \theta }{n_{\pm }\sin \theta }-\beta \frac{\cos {\varphi }_s}{\sin \theta }\right)J_s(x_s)\right]}^2\\ & \cdot \frac{g(\gamma ,\varphi )}{\beta (1+a_{\theta \pm }^2)}\frac{2\pi }{\cos \theta }\frac{\nu }{{\nu }_B},\end{array}$$

$$\begin{array}{r}{H}_{\pm }(\nu ,\theta ,\gamma ,s)=G_{\pm }(\nu ,\theta ,\gamma ,s)\cdot h_{\pm }(\nu ,\theta ,\gamma ,{\varphi }_s)\end{array}$$

$$\begin{array}{rl}\cos {\varphi }_s& :=\frac{1-\frac{s{\nu }_B}{\gamma \nu }}{n_{\pm }\beta \cos \theta },\\ x_s& :=\frac{s{n}_{\pm }\beta \sin \theta \sin {\varphi }_s}{1-n\pm \beta \cos \theta \cos {\varphi }_s}.\end{array}$$

$$j_{\pm }(\nu ,\theta )=\frac{e^{3}B{n}_e}{m_e{c}^2}\sum _{s=s_{\text{min}}}^{s_{\text{max}}}{\int }_{{\gamma }_{\text{min}}}^{{\gamma }_{\text{max}}}d\gamma G_{\pm }(\nu ,\theta ,\gamma ,s),$$

$${\kappa }_{\pm }(\nu ,\theta )=\frac{4\pi e^2{n}_e}{B}\sum _{s=s_{\text{min}}}^{s_{\text{max}}}{\int }_{{\gamma }_{\text{min}}}^{{\gamma }_{\text{max}}}d\gamma H_{\pm }(\nu ,\theta ,\gamma ,s),$$

$$\begin{array}{rl}{s}_{\text{min}}& :=\lfloor \frac{\gamma \nu }{{\nu }_B}(1-n_{\pm }\beta \cos \theta )\rfloor +1,\\ s_{\text{max}}& :=\lfloor \frac{\gamma \nu }{{\nu }_B}(1+n_{\pm }\beta \cos \theta )\rfloor ,\\ {\gamma }_{\text{min,max}}& :=\frac{\frac{s{\nu }_B}{\nu }\pm n_{\pm }\cos \theta {\left[{\left(\frac{s{\nu }_B}{\nu }\right)}^2+n_{\pm }^2{\cos }^2\theta -1\right]}^{1/2}}{1-n_{\pm }^2{\cos }^2\theta },\end{array}$$

$$\nu >{\left({\nu }_p^2+\frac{1}{4}{\nu }_B^2\right)}^{1/2}+\frac{1}{2}{\nu }_B$$

$${\kappa }_{\nu ,\text{ff}}$$

$$\times \{\begin{array}{l}(17.9+\ln T^{3/2}-\ln \nu ),(T<2\times 10^5\text{K})\\ (24.5+\ln T-\ln \nu ),(T>2\times 10^5\text{K})\end{array}$$

$$\sum _i{Z}_i{n}_i=n_p=n_{\text{textsc}Hii}$$

$$j_{\nu ,\text{ff}}=\frac{2k_{B}T{\nu }^2{k}_{\nu ,\text{ff}}}{c^2}$$

$$\kappa_{\nu,\textsc{H${}^{-}$}}=\frac{n_{e}n_{\textsc{H}}}{\nu}(A_{1}+(A_{2}-A_{3}/T)/\nu)(1-e^{-h\nu/k_{B}T})$$

$$j_{\nu ,\text{therm}}$$

$${\kappa }_{\nu ,\text{therm}}$$

$$\frac{d{I}_{\nu }}{d{\tau }_{\nu }}=-I_{\nu }({\tau }_{\nu })+S_{\nu }({\tau }_{\nu }),$$

$$T_b=\frac{c^2}{2k_B{\nu }^2}{I}_{\nu }$$

$$I_{\nu }$$

$${\tau }_{\nu }$$

$$\frac{V}{I}=\text{sign}(\cos \theta )\frac{T_{b,+}-T_{b,-}}{T_{b,+}+T_{b,-}+T_{b,\text{therm}}}.$$

$$F_{\nu }=\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}{I}_{ij}{\Omega }_{\mathrm{pixel}}.$$

$$I_{ij}=\frac{2k_b{\nu }^2}{c^2}{T}_{b,ij}$$

Full text

††pubyear: 2018††pagerange: Thyr: A Volumetric Ray-Marching Tool for Simulating Microwave Emission–Thyr: A Volumetric Ray-Marching Tool for Simulating Microwave Emission

1 Introduction and Background

Strong increases in microwave emission are observed during solar flares. Gyrosynchrotron emission, originating from mildly relativistic non-thermal electrons spiraling around magnetic field lines, is responsible for the majority of the emission in the 1-200 GHz range (Dulk, 1985; Bastian, 1998). Recent studies have also shown the importance of free-free emission from thermal electrons at these temperatures (e.g. Heinzel & Avrett, 2012). Observations of the microwave emission produced during a solar flare can be used to characterise the properties of the accelerated electrons and magnetic field of the flaring region.

Current observations of solar flares are typically carried out by radio telescopes without imaging capabilities (Sun-as-a-star) and by the interferometer Nobeyama Radio Heliograph (NoRH, Nakajima et al., 1994), capable of producing images at 17 and 34 GHz with moderate spatial resolution (10*′′* and 5*′′*, respectively) and with a temporal resolution down to 1 s. Such observations are capable of resolving the radio flare sources, associated with footpoints (e.g. Tzatzakis et al., 2008), flaring loops (e.g. Kundu et al., 2001), or ribbons (e.g. Simões et al., 2013) These features are commonly associated with the geometry of magnetic field loops, filled with electrons accelerated during the energy release phase of flares.

The first microwave spectral imaging analysis of a flare, with high spectral resolution, was performed by Wang et al. (1994), using the Owens Valley Radio Observatory (OVSA, Lim et al., 1994) 1–18 GHz data, and showed that the microwave spectral peak occurs at lower frequencies ($\nu<10$ GHz) for the looptop sources and at higher frequencies ($\nu>10$ GHz)in the footpoints, following the strength of the magnetic field in a flaring loop. Another example of a microwave imaging spectroscopic analysis was presented by Wang et al. (1995), also using OVSA data. A similar study was recently performed by Gary et al. (2018) with the Extended Owens Valley Solar Array (EOVSA, Kuroda et al., 2018; Wang et al., 2015) observations of the SOL2017-09-10 X8.2 event and yielded results consistent with those of Wang et al. (1994) and Wang et al. (1995), and simulation results of Simões & Costa (2006). Nindos et al. (2000) used one-dimensional modelling to reproduce the spatially resolved emission of a flare loop in 5 and 15 GHz, showing the 15 GHz emission was produced only in the footpoints and the 5 GHz emission restricted to the upper loop. To fit the observations it was necessary to use a much stronger magnetic field in the photosphere than at the looptop (870 G in the footpoints and 280 G in the looptop). Analysis of similar events by Lee & Gary (2000) provided evidence for magnetic trapping that would restrict the 5 GHz emission to the upper loop while the 15 GHz emission comes from higher energy electrons that are able to escape mirroring effects.

The complexity of solar microwave simulations has increased significantly in recent years, primarily building on the three-dimensional modeling tool GS Simulator (now GX Simulator) (Kuznetsov et al., 2011; Nita et al., 2015). GX Simulator computes solar radio emission from a three dimensional model, with in-built tools for magnetic field extrapolation from photospheric magnetograms, and a chromospheric approximation computed using the semi-empirical atmospheres of Fontenla et al. (2009), based on the photospheric intensity and magnetogram values (Fleishman et al., 2015). This tool has been used to forward fit and reconcile observations of radio and hard X-ray emission with a unified electron population (Kuroda et al., 2018), and investigate the plasma heating mechanism during a “cold” flare event using Linear Force Free Field extrapolation on magnetogram data to reconstruct the observed two-loop geometry and explain the heating delay by electron trapping in the looptop (Fleishman et al., 2016). Gordovskyy et al. (2017) used GX Simulator to investigate the polarised microwave emission from relaxing twisted coronal loops based on time-dependent magnetohydrodynamic simulations. These loops were found to produce gyrosynchrotron emission with short-term gradients of circular polarisation (cross-loop polarisation gradients), that depend strongly on viewing angle, and would primarily be visible on a loop observed on the solar limb, clearly showing the three-dimensional nature of the problem.

With the arrival of new and upgraded solar observatories, such as the Atacama Large Millimetric-submillimetric Array (ALMA, Wedemeyer et al., 2016) and EOVSA that are providing higher spatial and spectral resolution it is essential to have detailed predictions and models for these observations. While solar observations with ALMA have started (White et al., 2017; Iwai et al., 2017; Brajša et al., 2017), the Sun is now at its period of minimum activity, and no flares have been detected with ALMA yet. However, ALMA’s capability for flare studies have been proven with the detection of a super-flare on Proxima Centauri (MacGregor et al., 2018).

Regular solar observations in the millimetric range have started in 1999 with the Solar Submillimeter Telescope (SST, operating at 212 and 405 GHz Kaufmann et al., 2008). Since then a number of solar flares have been detected, with the typical decreasing, non-thermal spectrum towards higher frequencies (Giménez de Castro et al., 2009; Giménez de Castro et al., 2013), evidence for thermal (free-free) emission (especially at high frequencies towards the infrared range) (Heinzel & Avrett, 2012; Simões et al., 2017; Trottet et al., 2012, 2015), and with an increasing spectrum component at millimetric wavelengths (sub-THz) (e.g. Kaufmann et al., 2004), which was also observed with the Köln Observatory for Submillimeter and Millimeters Astronomy (KOSMA, Lüthi et al., 2004). This sub-THz component still remains without a satisfactory physical explanation (Krucker et al., 2013).

Previous observations have shown that flaring ribbons and footpoints can reach temperatures in the range of 1-10 MK (Hudson et al., 1994; Mrozek & Tomczak, 2004; Graham et al., 2013; Simões et al., 2015). The chromospheric plasma is normally opaque to radio emission (in the GHz range), and it would therefore be irrelevant if non-thermal electrons penetrate this layer. However, if this plasma is heated to these greater temperatures, then the free-free opacity drops and the plamsa becomes optically thin to any gyrosynchrotron emission potentially produced in the chromosphere, as discussed in Fletcher & Simões (2013).

Herein we present a new tool, Thyr, to compute the microwave emission from a dipole loop of plasma, containing a spatially variable atmosphere. Emission maps are computed under the assumption of local thermodynamic equilibrium (LTE) in the three-dimensional dipole volume. Full spectral information is available for each simulated pixel and can also be integrated over larger regions.

Thyr’s novel feature is its ability to manually refine the simulation resolution in the chromosphere while maintaining a complete microwave radiation treatment, allowing it to resolve the much smaller scales over which the atmospheric parameters evolve helping to account for both the free-free and gyrosynchrotron emission from this thin layer. Our model also supports using arbitrarily varying atmospheric parameters, magnetic field configurations, and electron pitch angle distributions.

Thyr builds on several generations of gyrosynchrotron modelling tools, including, but not limited to, Klein & Trottet (1984), Alissandrakis & Preka-Papadema (1984), Holman (2003), Simões & Costa (2006, 2010). Costa et al. (2013) and Cuambe et al. (2018) have both constructed large databases of three-dimensional gyrosynchrotron simulations from a large parameter space to develop an understanding of solar bursts and attempt to infer flare parameters from observation respectively. GX Simulator (Nita et al., 2015) can produce three-dimensional simulations of gyrosynchrotron emission from a dipole loop and has been used to investigate the effects of varying electron pitch-angle distribution (Kuznetsov et al., 2011), but focuses primarily on coronal microwave emission. These tools, including Thyr, build on the analytic expressions describing gyrosynchrotron emission formulated by Ramaty (1969).

In this paper we describe the functioning of the Thyr tool, validate its output against simulations presented in Klein & Trottet (1984), and use new simulations demonstrating the importance of modelling the lower atmosphere at high-resolution to capture the fine structures of the chromosphere. The source code and examples are available on Github at http://github.com/Goobley/Thyr2 (Osborne, 2018).

2 Theory

2.1 Gyrosynchrotron Emission

Gyrosynchrotron emission depends on a large number of parameters of the emitting region. These are primarily (Ramaty, 1969; Stähli et al., 1989)

• •

Magnetic field strength $B$: Directly determines the electron gyrofrequency.

• •

Plasma density $n_{p}$: Describes the density of electrons in the thermal plasma and can have a large effect on measured emission due to free-free emission and absorption effects. Dense magnetised plasmas may also strongly suppress the gyrosynchrotron emission, an effect known as Razin suppression (Ramaty, 1969).

• •

Non-thermal electron density $n_{e}$: This parameter is largely responsible for the strength of the emission, as a higher electron density will produce more radiation. It also affects the importance of the gyrosynchrotron self-absorption, via radiative transfer.

• •

Non-thermal electron distribution $g(\gamma,\phi)$: This function describes the distribution of non-thermal electron energies (here in terms of relativistic Lorentz factor $\gamma$) and pitch-angles – the angle between the electron’s velocity vector and the magnetic field. It is often assumed that the electron energies follow a single power law distribution determined by their spectral index $\delta$, $f(\gamma)\propto(\gamma-1)^{-\delta}d\gamma$, and minimum and maximum cut-off energy values, $\gamma_{\mathrm{min}}$ and $\gamma_{\mathrm{max}}$, respectively. The distribution of pitch-angle $\phi$ has a strong effect on the angles radiation is observed at due to the strong beaming effect of the radiation from relativistic electrons. Thyr supports multi-power law energy distributions and arbitrary distributions of pitch-angle.

• •

Viewing angle $\theta$: The angle between the wave vector and the magnetic line has a strong effect on the observed radiation due to the polarisation of the radiation and the beaming of emission from relativistic particles.

As gyrosynchrotron emission is produced by electrons spiraling around the magnetic field, it is composed of two circularly polarised modes. We designate the mode extraordinary when the circular polarisation is right-handed and the electric field vector rotates in the same direction as the electrons. In the opposite case the mode is called ordinary, or left-handed. Following the convention of Klein & Trottet (1984) we use “$-$” to refer to the extraordinary mode (right-handed polarisation for positive $B$) in equations, and “$+$” to refer to the ordinary mode (left-handed).

The following treatment is based on that of Klein & Trottet (1984), developed from Ramaty (1969) and the corrections of Trulsen & Fejer (1970). The refractive index of the plasma is given by the Appleton-Hartree equation (e.g. Stix (1962))

[TABLE]

where

[TABLE]

The polarisation coefficient $a_{\theta}$ is then

[TABLE]

and is used to compute the spectral intensity $\epsilon$ of a single electron (Trulsen & Fejer, 1970)

[TABLE]

where

[TABLE]

The electron distribution is described by the function $g(\gamma,\phi)$ such that

[TABLE]

This gives rise to the expression of the emission and absorption coefficients, $j$ and $\kappa$ respectively, for a homogeneous ensemble of electrons with density $n_{e}$ and energy and pitch-angle distribution $g$

[TABLE]

[TABLE]

where

[TABLE]

The presence of the delta function in $\epsilon$ allows the integral over $\cos{\phi}$ to be solved analytically for $\theta\neq\pi/2$.

Following Ramaty (1969) we define $G_{\pm}$ and $H_{\pm}$

[TABLE]

[TABLE]

where

[TABLE]

Solving the integral over $\cos{\phi}$ analytically and swapping the resultant integral over $\gamma$ with the summation over $s$ by following the approach of Holman (2003) to equations (5) and (6) yields

[TABLE]

[TABLE]

where

[TABLE]

and $\lfloor x\rfloor$ is a function that returns the floor of $x$.

In reality, we do not always compute the summation up to $s_{\textrm{max}}$ as $G_{\pm}$ and $H_{\pm}$ may have converged for a smaller $s$. This is determined by the relative change of $G_{\pm}$ and $H_{\pm}$ across successive iterations of the summation.

Additionally, ordinary-mode emission may only be produced when $\nu>\nu_{p}$ and $n_{+}>0$. Similarly, extraordinary-mode emission requires

[TABLE]

and $n_{-}>0$.

In our implementation the integrals over $\gamma$ are computed using a Gauss-Legendre method. Additionally, significant speed-ups were obtained by using a look-up table for the common range of Bessel functions encountered during the simulation. This is very efficient because the computation of the Bessel function is by far the dominant computational cost in the gyrosynchrotron coefficients and similar regions of the Bessel function will be accessed sequentially, allowing the CPU’s caching and pre-fetching mechanisms to work efficiently. If a Bessel function outside the tabled space is requested, and the approximation is valid, then the expression from Wild & Hill (1971) is used.

The expressions for $G_{\pm}$ (8) and $H_{\pm}$ (8) do not hold for $\theta=90^{\circ}$, whilst it is possible to derive additional expressions for this case, there is little reason to, as it only holds for this singular case, and also for reasons discussed in Sec. 2.3.

2.2 Thermal Emission

In addition to the radio gyrosynchrotron emission parameters, we compute the thermal free-free emission and absorption coefficients, as well as the H- opacity, known to be relevant for submillimetric emission (Heinzel & Avrett, 2012; Simões et al., 2017).

For the free-free opacity we follow Dulk (1985), with the correction from Wedemeyer et al. (2016)

[TABLE]

where $n_{e}$ is the thermal electron number density, $n_{i}$ is the number density of ion $i$ with charge $Z_{i}$, $\nu$ is the frequency, and $T$ is the temperature of the plasma. Herein we assume a uniform hydrogen plasma such that

[TABLE]

where $n_{p}$ is the proton density and $n_{\textsc{Hii}}$ is the ionised hydrogen density.

Now, by Dulk (1985) (and Kirchoff’s law) we have the free-free emission coefficient

[TABLE]

where $k_{B}$ is Boltzmann’s constant.

We follow the treatment of Kurucz (1970) computing the H- opacity

[TABLE]

where $n_{\textsc{H}}$ is the neutral hydrogen density, $A_{1}~{}=~{}1.3727~{}\times~{}10^{-25}$, $A_{2}~{}=~{}4.3748~{}\times~{}10^{-10}$, and $A_{3}~{}=~{}2.5993~{}\times~{}10^{-7}$.

Finally, we have the thermal emission and absorption coefficients

[TABLE]

where $B_{\nu}(T)$ is the Planck function. At solar gyrosynchrotron frequencies free-free opacity always dominates over H$-$ emissivity, and therefore the effects of H$-$ emissivity were ignored during the numerical simulations.

2.3 Radiative Transfer

To compute the emission maps generated by this process we compute radiative transfer under the assumption of local thermodynamic equilibrium. The monochromatic observed intensity $I_{\nu}$ is then given by the radiative transfer equation

[TABLE]

where $\tau_{\nu}$ is the monochromatic optical length of the plasma at frequency $\nu$ along the photon’s path between emission and absorption, and $S_{\nu}$ is the monochromatic source function where $S_{\nu}=j_{\nu}/\kappa_{\nu}$.

Using the ray-marching algorithm of Section 3.2 these emission and absorption parameters are combined, from back to front along the observer’s line of sight, so as to produce a brightness temperature map for each radio mode and the thermal emission ($T_{b,+}$, $T_{b,-}$, and $T_{b,\textrm{therm}}$ respectively). The brightness temperature is computed for each pixel of each map as per the Rayleigh-Jeans approximation

[TABLE]

where

[TABLE]

and the optical length travelled by a photon $\tau_{\nu}$ is related to the depth $h$ along the observer’s line of sight by $d\tau_{\nu}=-k_{\nu}\,dh$, where $k_{\nu}$ is the sum of the absorption coefficients from different processes affecting the frequency $\nu$. In these integrals, $H$ is the photon emission point along the observer’s line of sight, and 0 is the observer’s location.

A total emission map is then computed from $T_{b,+}+T_{b,-}+T_{b,\textrm{therm}}$.

2.4 Circular Polarisation Degree

When computing the total emission it is sufficient to simply add the flux across the emitting modes, however, for investigating polarisation it is necessary to include propagation effects and the radiative transfer of all four Stokes parameters. Following Simões & Costa (2006) and Cohen (1960), from magnetoionic theory in solar conditions there are two regimes for the propagation of circular polarisation (Stokes V), the quasi-longitudinal and quasi-transversal approximations.

Using the C7 semi-empirical atmosphere (Avrett & Loeser, 2008) and the dipole model discussed in Sec. 3.1 with a looptop magnetic field strength of 100 G the quasi-longitudinal approximation holds for $\theta$ outside the range $(87^{\circ}-93^{\circ})$ at $\tau=1$ optical depth from free-free absorption. For cells where $\theta$ falls within the $(87^{\circ}-93^{\circ})$ range we simply adjust the viewing angle to the closer edge of the range to retain the validity of the quasi-longitudinal approximation. Unless a very highly anisotropic pitch angle distribution is used, that is peaked within this range, there are no significant differences in the intensities of the two modes due to the application of this approximation.

From the optical depth calculation we find that 10 GHz radiation does not penetrate through the transition region, 45 GHz radiation just enters the chromosphere, and 200 GHz radiation reaches $\sim$400 km into the chromosphere. In regions deeper in the chromosphere, for viewing angles close to perpendicular, the quasi-longitudinal approximation may fail, and one should be wary of interpreting the polarisation results from this region. The calculation of the total intensity is, however, unaffected by the use of this approximation.

Following the creation of emission maps we can compute the circular polarisation degree of the radiation. We follow the convention of Klein & Trottet (1984) and define for the Stokes $V$ and $I$ in the frame of the wave

[TABLE]

In the reference frame of an observer the circular polarisation will be reversed such that $V_{\textrm{obs}}=-V_{\textrm{wave}}$.

In Thyr the magnetic viewing angle term used for calculating the polarisation degree is simply computed as an average along the ray, and so caution should be taken with the interpretation of polarised radiation from overlapping sources with opposing magnetic field orientation. The viewing angles discussed in the previous sections are all treated without this approximation.

3 The Thyr Model

3.1 Magnetic field geometry

We adopt the same analytic dipole model as Kuznetsov et al. (2011, see Appendix) to describe the magnetic field geometry. This describes both the magnetic field at each point in the volume and the region contained within the dipole. This analytic dipole model is described in terms of the magnetic field at the centre of the looptop, the radius of the flux tube at the looptop, the height, and the submerged depth of the dipole. The magnetic dipole element at the base of the loop can also be inclined, and the ray-marching method allows the loop to be visualised from any angle which can be specified in terms of location on the Sun. A solar location is specified by four angles, solar latitude and longitude, tilt away from the local normal to the surface about an axis connecting the footpoints, and rotation about the local normal.

3.2 Ray-marching

The standard approach for computing the emission map of a three-dimensional volume is to use ray tracing, however this becomes significantly more computationally demanding as the number of discrete volume elements (voxels, assumed quasi-homogeneous three-dimensional cubic volumes) increases due to the number of ray-voxel intersection tests that need to be performed to determine whether a certain voxel interacts with a ray. This method can be optimised by using voxel culling methods, such as octrees, which describe layers of nested voxels, allowing many of the tests against the most refined voxels to be avoided. These optimisations are quite complex to implement, and given the simplicity of this problem are unnecessary.

Ray-marching instead assumes that the emission and absorption properties of the volume can be sampled continuously. Instead of computing the radiative transfer integral along a ray between voxel intersections whilst assuming homogeneous plasma parameters, when ray-marching the integral is computed without a grid, by moving a short step along the ray and sampling the local emission and absorption coefficients, whilst assuming homogeneity over this short step. If these steps taken along the ray are sufficiently short then this method will tend towards the true value of the integral.

Under the basic description of ray-marching given above, it is assumed that the emission and absorption coefficients of the plasma are defined continuously and can be sampled at any point. Due to the high cost of computing these coefficients for gyrosynchrotron radiation (Sec. 2.1) it is not feasible to recompute these at every step along the ray. The plasma emission and absorption coefficients for the two gyrosynchrotron modes and the thermal radiation are therefore computed on a discrete three-dimensional Cartesian grid and stored in a volume texture (three-dimensional cuboidal array) before the ray-marching procedure. The coefficients are also multisampled i.e. computed for multiple points within each voxel and then averaged to attempt to better capture the average plasma conditions than simply the conditions at the central point.

There are multiple methods for determining when the ray-marching step size needs to be decreased or can be increased. When the absorption and emission coefficients are continuously defined then metrics based on the local gradient can be used. In Thyr we choose to perform the refinement manually. The bounds of the volume texture containing the plasma coefficients define an axis-aligned bounding box (AABB) for the object stored in the texture. Then, using the simple but highly optimised, “slab” algorithm (Kay & Kajiya, 1986) to compute the intersection of rays with this AABB we have the start and end points between which to ray-march. This “slab” algorithm has been further optimised by relying on strict IEEE754 floating point behaviour. The manually specified regions of refinement are converted into AABBs contained within the initial AABB. The plasma coefficients are computed on a finer grid (with user specified refinement factor) within these sub-domains. If a ray intersects with one of these higher resolution regions (determined by the intersection of the ray with the sub-domain’s AABB) then the step size is decreased and the plasma coefficients are sampled from the finer grid when computing the radiative transfer integral.

Thyr uses three volume textures to store the coefficients for the dipole volume. The primary texture represents the entirety of the volume at the (lower) coronal resolution, whilst the other two are heavily refined on the footpoints of the dipole encompassing the photosphere, chromosphere, and the transition region as the atmospheric parameters and magnetic field change much more rapidly in this region. Whilst a dipole model is chosen for our example, any volumetric function or data (such as the results of a magnetogram extrapolation) can be used to fill the parameters of the texture, and the ray-marching approach will remain unchanged. The model is currently based around the concept of a single AABB containing the entire scene at low-resolution and a number of refined regions contained strictly within this AABB. The information from the low-resolution AABB is then ignored in the locations where a refined region is present. These AABBs serve as the boundaries of a two-level three-dimensional Cartesian grid hierarchy, with a uniform grid within each level. If multiple low resolution AABBs are desired this change would be relatively trivial, but using the convention of a single low resolution AABB allows any geometry to be traced by Thyr with no changes to the ray-marching code. It would equally be simple to extend the code to support a full multi-level grid hierarchy with further refined regions within the refined regions, if extreme resolution were required in some locations.

Fig. 1 shows (in two dimensions) the difference between the ray-tracing and ray-marching. The cost of finding the emission and absorption parameters for a point on the ray is constant if these parameters are described by a simple expression or discretised onto a known Cartesian grid, as these parameters can be computed or looked-up wherever the blue point happens to be. When applying a traditional, simple, ray-tracing technique without additional voxel culling optimisations, the ray-voxel intersection test must be performed against every voxel in the domain to determine the locations of the green points. This computation scales linearly with the number of voxels. Therefore, using the ray-marching approach, it is computationally tractable to sample the volume with a step size significantly smaller than the voxel side length. The integral then amounts to simply performing several elementary mathematical operations for each step to compute the closed form radiative transfer integral from a homogeneous slab, and looking up one value in a large array. Due to the cost of calculating the many ray-voxel intersections, for the cases used in Thyr ray-marching is typically 1.5 orders of magnitude faster than ray-tracing, even though the integral of radiative transfer is computed at significantly more points including a higher resolution region (orange box in Fig. 1).

An assumption that is implicit to the ray-marching approach is that it is acceptable to not sample voxels for which the space between the intersections is very small (e.g. just cutting through close to a corner), as their effect is insignificant. In practice with a sufficiently small step-size (Thyr uses 0.1$\times$ the local grid side length for the plasma coefficients) and reasonably smoothly varying emission and absorption coefficients this is not a problem as the effect of the contribution from this region is not significant. The multisampling of the coefficients helps improve the smoothness of the coefficients between voxels and guarantee that the parameters chosen accurately reflect the plasma contained within (assumed homogeneous).

4 Verification against Klein & Trottet (1984)

To verify Thyr’s correctness we use a problem presented in Klein & Trottet (1984). In this test case the gyrosynchrotron emission from a two-dimensional loop is simulated. The loop is observed from two different viewing angles $\Lambda$. The specific intensity (sfu111solar flux units (1 sfu = $10^{4}$ Jy = $10^{-22}$ W m*-2* Hz*-1* = $10^{-19}$ erg s*-1* cm*-2* Hz*-1*)/sterad) is then computed for a telescope beam scanning across the loop with angle $\chi$. In Thyr the rays along the observer’s line of sight are assumed parallel, and do not diverge from a point like the telescope in Klein & Trottet (1984). The coordinate $x$ is then the displacement in Thyr’s imaging plane in units of arcminutes on the surface of the Sun viewed from Earth.

Using Klein and Trottet’s parameters for loop with 3 G coronal magnetic field strength at $\Lambda=0^{\circ}$ and $\Lambda=50^{\circ}$ we have reproduced the emission at 320 MHz in Fig. 2. The results are very similar in shape and intensity. The slight differences in intensity can be accounted for by modifying the looptop magnetic field strength by less than 5%. Klein & Trottet (1984) define the looptop magnetic field strength at the outer boundary of the looptop, whereas in Thyr it is defined in the centre. This value had to be manually tuned to obtain a magnetic field strength similar to that used by Klein & Trottet (1984) in the slice taken from the three-dimensional simulation performed in Thyr.

The spatial offset of Thyr’s results in the $\Lambda=50^{\circ}$ case is due to Thyr performing the rotation but also translating the dipole so as to keep it centered in the field of view.

5 Example 3D Simulations

To demonstrate Thyr’s capabilities we performed a set of flare simulations using a simple dipole model for the structure of the emitting volume and the magnetic field, a power-law non-thermal electron distribution, and the other atmospheric parameters specified by a semi-empirical model. From these simulations we obtain brightness maps, spectra, and polarisation maps. Complete spectral data is available for every pixel and region of these maps, thus allowing us to perform a localised spectroscopic analysis (Section 5.5). We perform these simulations both with and without resolving the chromosphere in detail to show its role on emission at higher frequencies. These simulations can serve as examples for how to set up and run the code, in addition to basic post-processing of results.

5.1 Semi-empirical atmosphere and parameters

In our example simulations we use the C7 quiet Sun semi-empirical atmosphere of Avrett & Loeser (2008). This atmosphere is extrapolated to coronal parameters, by extrapolating linearly in log-log space the parameters from the top of the C7 model (at 47 Mm), until 150 Mm where the atmosphere is then set constant. Our atmosphere is plotted in Fig. 3.

A non-thermal electron distribution $g(E,\phi)$ is embedded in the semi-empirical atmosphere, with a density defined as $n_{e}=10^{6}$ cm*-3* in the corona, and $n_{e}=10^{8}$ cm*-3* at the top of the chromosphere. This increase of $n_{e}$ is caused by the decrease of the emitting volume due to the convergence of the dipole magnetic field. As can be seen Fig. 3, we have made the approximation to make this increase a step rather than scaling directly with the area of the dipole. This decision was made to ensure a simple example using a dipole magnetic field that produces chromospheric radio emission and to provide an electron number flux of $~{}10^{36}$ s*-1* entering the chromosphere, which is consistent with observations (e.g. Simões & Kontar, 2013). These electrons follow a power-law distribution in energy, with a minimum cut-off energy of $E_{\mathrm{min}}=10$ keV, a maximum cut-off of $E_{\mathrm{max}}=5$ MeV, and a spectral index of $\delta=3$. We assume that these electrons have an isotropic pitch-angle distribution, but the code supports arbitrary pitch-angle distributions.

As we descend into the chromosphere the non-thermal electron distribution is attenuated by balancing collisional losses against the column density traversed by the electrons. We use the approximate relationship that the critical stopping column density is $N_{\textrm{stop}}\approx 10^{17}E_{0}^{2}$ where $E_{0}$ is the initial electron energy in keV and $N_{\textrm{stop}}$ is in cm*-2* (Tandberg-Hanssen & Emslie, 1988). The lower energy bound of the power law is then proportionately increased as electrons below the stopping energy are cut off.

Finally, for these models we use a dipole with height of approximately $3.86\times 10^{9}$ cm ($53.3"$) and a centre-to-centre footpoint separation of $2.57\times 10^{9}$ cm ($35.5"$). The radius at the looptop ($\rho_{0}$ in the Kuznetsov et al. (2011) model) is $7.72\times 10^{8}$ cm. The looptop field is set to 100 G, yielding a strength of the order of 2000 G at the photosphere. The dipole here has a latitude of 30∘, longitude of 70∘, and rotation about the local normal of $-20^{\circ}$.

The simulations presented here are performed both with and without refined resolution in the chromosphere. In the low resolution simulation the voxels have a side length of 1800 km. In the high resolution simulation the voxels of the upper corona have a 450 km side length and the refined voxels in the lower atmosphere have a 75 km side length (refined six times). Thyr’s ray-marching algorithm automatically transitions between the low- and high-resolution regions with no ill effects due to the boundary. An additional high resolution simulation was also computed with voxels with 100 km side lengths, and no significant spectral difference was found between the two, so it was concluded that the 75 km resolution was sufficient to capture the details of the problem.

5.2 Emission Maps

Fig. 4 shows the total brightness maps at 35 GHz for the same loop simulated at different resolutions. The emission from the high-resolution loop at 2 GHz is shown in Fig. 5. This shows the complex nature of gyrosynchrotron emission at frequencies at which the plasma is optically thick. The high frequency emission maps (Fig. 4) have the most intense emission from the footpoints, associated with the strongest magnetic field values. In the low frequency emission maps (Fig. 5) the stripes are caused by the harmonic structure of the gyrosynchrotron emission at a single frequency for a spatially varying $B$ field. The effects of insufficient resolution within the simulation can be seen by comparing the simulations in Fig. 4. The low-resolution simulation of Fig 4a produces an lower brightness temperature than the high-resolution model shown in Fig. 4 due to the small region this high brightness originates from.

5.3 Spectra

Integrating over the brightness temperature maps, the flux density spectra $F_{\nu}$ of these two simulations can be computed. In practical terms, this calculation is simply the sum of the flux density of each pixel (column $i$ and row $j$), obtained from their specific intensity $I_{ij}$ over the pixel solid angle $\Omega_{\mathrm{pixel}}$:

[TABLE]

The specific intensity $I_{ij}$ of the pixels are directly found from their brightness temperature values $T_{b,ij}$ via the Rayleigh-Jeans law:

[TABLE]

The resulting spectra are shown in Fig. 6, in sfu. Both of these plots show the characteristic “inverted-v” shape of gyrosynchrotron emission, and the fact that the ordinary mode dominates at low frequencies and the extraordinary mode dominates at higher frequencies (Ramaty, 1969). The thermal emission is plotted, but is insignificant in the cases we are investigating here. The difference in the magnitude of the thermal emission is due to a region that is small and bright in the high-resolution case, being spread across a much larger voxel in the low-resolution case. Although the spectrum of gyrosynchrotron radiation from a uniform source would present harmonics, they are not visible in the spatially-integrated spectra due to the overlap of harmonic structures from a spatially-varying $B$ field and atmospheric parameters, as has been previously discussed by Klein & Trottet (1984); Simões & Costa (2006). The spectra for the low- and high-resolution cases behave similarly at low frequencies, but differences emerge at higher frequencies. In the high-resolution simulation there is a hardening of the spectrum of the ordinary and extraordinary modes of the gyrosynchrotron radiation at high frequencies. As these non-thermal modes dominate over the thermal emission the hardening in the total spectrum must be non-thermal in origin, as we will discuss in Section 7.

5.4 Polarisation Degree

Fig. 7 shows a map of the polarisation degree for the high resolution model at 17 GHz. The polarisation degree is given by the Stokes parameters $V/I$. The expected opposite polarisation in each footpoint is present. The asymmetry in peak polarisation degree between the two footpoints is purely an effect of viewing angle – the loop is fully symmetric. Klein & Trottet (1984) and Simões & Costa (2006) have previously shown the importance of the viewing angle on observed polarisation, and our results are in accordance. With a three-dimensional system there is a variety of situations that can produce significant differences between the intensity of footpoints, and the location of the change in polarisation direction. Polarisation data is an important component of radio observation, one that could be used to shed light on the magnetic field geometry (Simões & Costa, 2006), and an essential tool for diagnosing the electron energy distribution (Kuznetsov et al., 2011).

The structure of the polarised radiation shown in Fig. 7 is quite simple with a single transition between negative and positive polarisation. At lower frequencies where a large portion of the volume is optically thick and the complex “gyrostripe” patterns (see Fig. 5) are visible the polarisation patterns are also far more complex and follow these stripes.

In the simulations presented here a simple dipole magnetic field model is used. As the polarisation degree is strongly related to the direction of the magnetic field a more complex magnetic model can yield very different polarisation patterns (e.g. Gordovskyy et al., 2017). Thyr is capable of using an arbitrary magnetic field geometry and thus can be used to investigate cases with complex polarisation patterns, including fine structure if the resolution is increased in the region of interest.

5.5 Imaging Spectroscopy

In Fig. 5 we have marked a looptop and a footpooint region in teal and red, respectively. The central pixel is marked in teal for the looptop region, and multiple pixels are marked and numbered inside the footpoint region. The simulated spectral flux from the marked pixels and integrated over the regions is plotted in Fig. 8, where the colours in Fig. 8a indicate the total emission from the pixel of the same colour (not separated by emission mode). In the footpoint region no extraordinary-mode emission is seen at less than 2 GHz as the extraordinary-mode is unable to propagate due to the ratio of the plasma frequency to the gyro-frequency being too low (see (12) and Ramaty (1969)). It can be clearly seen that the peak emission occurs at a significantly lower frequency in the looptop than in the footpoint. Additionally, the peak flux emitted from the footpoint region is much larger than that from the looptop. The high frequency emission from the looptop drops off as expected for a homogeneous gyrosynchrotron source. This is not so for the footpoint region. Fig. 8b clearly shows that the footpoint is the origin of this hardening, which is logical as this effect is only seen in simulations with a heavily refined lower atmosphere. By investigating the spectra of the marked pixels within the footpoint region, we can see that the peak emission frequency continues to increase significantly with depth in the atmosphere, and that pixels 3, 4, and 5 lie within the compact region of the atmosphere from which the high-frequency hardening originates.

The results presented here were simulated with an isotropic electron pitch-angle distribution. As described in Section 5.1 the lower energy electrons are prevented from reaching the lower atmosphere as they are removed from the distribution above a critical column density. If a loss cone style distribution of pitch angles were also applied the low–energy end of the distribution could drop off faster in the atmosphere, further enhancing this high-frequency hardening in the footpooints (Simões & Costa, 2010; Kuznetsov et al., 2011). Arbitrary pitch angle distributions are supported by the gyrosynchrotron code in Thyr and this effect could therefore be easily investigated.

6 Code Description

Thyr requires a modern C++14 compliant compiler, Torch7 (built on LuaJIT), and Python 3 with the matplotlib package. The core calculation of the gyrosynchrotron components is computed in C++, whilst all set up, and ray-marching is performed in Lua. Plotting is done through our Plyght tool, also available on Github 222https://github.com/Goobley/Plyght, that allows any language with a socket API or foreign function interface to easily produce two-dimensional plots using matplotlib (Hunter, 2007). The code also supports exporting all data to the widely supported comma separated variable (csv) format to allow for post-processing and plotting in any language.

7 Discussion

The Thyr tool, its initial validation, and example simulations have now been presented. The high-resolution simulation presented here is computed at much higher resolution than is available with modern radio observatories. This fine structure remains nonetheless important due to the contribution of the radiation produced in these regions for the total beam power. This can be seen by looking at the image that an observatory such as Nobeyama would produce from the different models. These convolved maps are shown in Fig. 9. The maps shown in these figures appear very similar, but the peak brightness temperature is different, and the magnitude of this difference only increases at higher frequencies, as demonstrated by the spectra (e.g. Fig 6).

Using the spectral imaging techniques of Section 5.5 we have identified the location from which the hardening of the high frequency emission originates. Points 3, 4 and 5 from Fig 5 lie within a thin ($\sim$300 km) layer in the upper chromosphere and transition region. This corresponds to a layer within the region of 1.8–2.15 Mm on Fig. 3. The field chosen for this simulation was a dipole as this represents the simplest analytical case and a realistic flaring active region can be expected to have stronger magnetic convergence. In this simulation the average field along the line of sight for each of these pixels is $\sim$1.3 kG and they lie within a region of strong magnetic convergence. This produces the total hardening at high frequencies due to the superposition of the small-scale hardenings produced by emission from small regions with high magnetic field strength in the core of this thin slab.

Given its characteristic spectral signature, there is no need to spatially resolve this thin layer. If the emission from the footpoints cannot be resolved separately then a similar spectrum could also be obtained from a loop with different magnetic field strengths in the footpoints (e.g. with a significantly rotated magnetic moment). To verify if the emission behaves in this way increased spectral resolution at high frequencies is required. This can be achieved with today’s technology using simultaneous observations from Radio Solar Telescope Network (RSTN, operating up to 15.4 GHz, Guidice, 1979), Nobeyama Radio Polarimeters (NoRP, Nakajima et al., 1985), Nobeyama Radio Heliograph (NoRH, Nakajima et al., 1994) POlarization Emission of Millimeter Activity at the Sun (POEMAS, operating at 45 and 90 GHz Valio et al., 2013), ALMA (currently operating at 100 and 230  GHz Wedemeyer et al., 2016), and the Submillimeter Solar Telescope (SST, 212 and 405 GHz Kaufmann et al., 2008). The combination of these instruments could provide valuable information on the shape of the spectrum at higher frequencies and hence an estimate of the structure of the magnetic field in the atmosphere.

8 Conclusions

We have presented three-dimensional solar flare radio emission simulation tool Thyr, verification against Klein & Trottet (1984), and example simulations that demonstrate the importance of resolving the lower atmosphere to a much higher resolution than used in previous models. This tool is now available for use under the MIT license and can be downloaded at http://github.com/Goobley/Thyr2 (Osborne, 2018)

Thyr has the ability to simulate user-specified regions with increased accuracy, and we use this to increase the resolution in the lower corona, chromosphere and photosphere. By performing simulations with a high resolution lower atmosphere we have shown that a non-thermal hardening of the spectrum should be expected at higher frequencies from a dipole loop. Using a combination of RSTN, POEMAS, and ALMA, the existence of such a hardening could be investigated.

Recent studies have also argued for the importance of thermal free-free emission at higher frequencies (Simões et al., 2017; Heinzel & Avrett, 2012). Whilst this emission is not important in the model we have selected here, Thyr could well be used to investigate the parameters for which it becomes significant. The C7 atmosphere was chosen here due to its ubiquity and high-resolution – there is little reason to perform a high resolution simulation with a low resolution atmosphere! It is simple to investigate the influence of other atmospheres using Thyr and the files present in the github repository can serve as a base. For example, Thyr can accept snapshots of the flaring atmosphere generated by radiative hydrodynamic simulations, such as RADYN (Carlsson & Stein, 1995; Allred et al., 2015) and Flarix (Varady et al., 2010; Heinzel et al., 2017).

Thyr can also serve as a skeleton for future local thermodynamic equilibrium radiative transfer codes as it is simple to replace the geometry and/or use different emission and absorption coefficients. Three-dimensional modelling is especially useful in cases where the emission has a high angular dependence, such as the case with gyrosynchrotron here. We therefore hope that this code can be modified and be of use to the wider astronomical community.

Acknowledgements

The authors would like to thank Lyndsay Fletcher for helpful comments and suggestions. C.M.J.O is grateful for the financial support of the Royal Society of Edinburgh Cormack Undergraduate Vacation Research Scholarship (2016), and the Royal Astronomical Society Undergraduate Summer Bursary (2015) under which this project took shape. C.M.J.O also acknowledges financial support from STFC Doctoral Training Grant ST/R504750/1. P.J.A.S. acknowledges support from the University of Glasgow’s Lord Kelvin Adam Smith Leadership Fellowship. This work builds upon the results obtained from projects funded by FAPESP grants 03/03406-6, 04/14248-5, 08/09339-2 and 2009/18386-7. We acknowledge the use of colour-blind safe and print-friendly colour tables by Paul Tol (https://personal.sron.nl/~pault/). The authors would also like to thank the reviewer for detailed comments and suggestions for improvement.