Thyr: A Volumetric RayMarching Tool for Simulating Microwave Emission
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 threedimensional flaring loops with spatially varying atmosphere and increased resolution in the lower corona and chromosphere. Thyr is modular and opensource, consisting of separate components to compute the thermal and nonthermal microwave emission coefficients and perform threedimensional radiative transfer (in local thermodynamic equilibrium). The radiative transfer integral is computed by a novel raymarching 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 highresolution performed with Thyr, with a spectral imaging analysis of differing regions. The highresolution 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.
keywords:
Sun: flares – Sun: radio radiation – radiative transfer – methods: numerical1 Introduction and Background
Strong increases in microwave emission are observed during solar flares. Gyrosynchrotron emission, originating from mildly relativistic nonthermal electrons spiraling around magnetic field lines, is responsible for the majority of the emission in the 1200 GHz range (Dulk, 1985; Bastian, 1998). Recent studies have also shown the importance of freefree 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 (Sunasastar) 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 ( GHz) for the looptop sources and at higher frequencies ( 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 SOL20170910 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 onedimensional 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 threedimensional 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 inbuilt tools for magnetic field extrapolation from photospheric magnetograms, and a chromospheric approximation computed using the semiempirical 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 Xray 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 twoloop 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 timedependent magnetohydrodynamic simulations. These loops were found to produce gyrosynchrotron emission with shortterm gradients of circular polarisation (crossloop polarisation gradients), that depend strongly on viewing angle, and would primarily be visible on a loop observed on the solar limb, clearly showing the threedimensional nature of the problem.
With the arrival of new and upgraded solar observatories, such as the Atacama Large Millimetricsubmillimetric 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 superflare 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, nonthermal spectrum towards higher frequencies (Giménez de Castro et al., 2009; Giménez de Castro et al., 2013), evidence for thermal (freefree) 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 (subTHz) (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 subTHz 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 110 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 nonthermal electrons penetrate this layer. However, if this plasma is heated to these greater temperatures, then the freefree 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 threedimensional 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 freefree 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 & PrekaPapadema (1984), Holman (2003), Simões & Costa (2006, 2010). Costa et al. (2013) and Cuambe et al. (2018) have both constructed large databases of threedimensional 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 threedimensional simulations of gyrosynchrotron emission from a dipole loop and has been used to investigate the effects of varying electron pitchangle 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 highresolution 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 : Directly determines the electron gyrofrequency.

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

Nonthermal electron density : 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 selfabsorption, via radiative transfer.

Nonthermal electron distribution : This function describes the distribution of nonthermal electron energies (here in terms of relativistic Lorentz factor ) and pitchangles – 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 , , and minimum and maximum cutoff energy values, and , respectively. The distribution of pitchangle 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 multipower law energy distributions and arbitrary distributions of pitchangle.

Viewing angle : 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 righthanded and the electric field vector rotates in the same direction as the electrons. In the opposite case the mode is called ordinary, or lefthanded. Following the convention of Klein & Trottet (1984) we use “” to refer to the extraordinary mode (righthanded polarisation for positive ) in equations, and “” to refer to the ordinary mode (lefthanded).
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 AppletonHartree equation (e.g. Stix (1962))
(1) 
where
The polarisation coefficient is then
(2) 
and is used to compute the spectral intensity of a single electron (Trulsen & Fejer, 1970)
(3) 
where
The electron distribution is described by the function such that
(4) 
This gives rise to the expression of the emission and absorption coefficients, and respectively, for a homogeneous ensemble of electrons with density and energy and pitchangle distribution
(5) 
(6) 
where
(7) 
The presence of the delta function in allows the integral over to be solved analytically for .
Following Ramaty (1969) we define and
(8) 
(9) 
where
Solving the integral over analytically and swapping the resultant integral over with the summation over by following the approach of Holman (2003) to equations (5) and (6) yields
(10) 
(11) 
where
and is a function that returns the floor of .
In reality, we do not always compute the summation up to as and may have converged for a smaller . This is determined by the relative change of and across successive iterations of the summation.
Additionally, ordinarymode emission may only be produced when and . Similarly, extraordinarymode emission requires
(12) 
and .
In our implementation the integrals over are computed using a GaussLegendre method. Additionally, significant speedups were obtained by using a lookup 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 prefetching 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.
2.2 Thermal Emission
In addition to the radio gyrosynchrotron emission parameters, we compute the thermal freefree 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).
(13)  
where is the thermal electron number density, is the number density of ion with charge , is the frequency, and is the temperature of the plasma. Herein we assume a uniform hydrogen plasma such that
(14) 
where is the proton density and is the ionised hydrogen density.
Now, by Dulk (1985) (and Kirchoff’s law) we have the freefree emission coefficient
(15) 
where is Boltzmann’s constant.
We follow the treatment of Kurucz (1970) computing the H opacity
(16) 
where is the neutral hydrogen density, , , and .
Finally, we have the thermal emission and absorption coefficients
(17)  
where is the Planck function. At solar gyrosynchrotron frequencies freefree 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 is then given by the radiative transfer equation
(18) 
where is the monochromatic optical length of the plasma at frequency along the photon’s path between emission and absorption, and is the monochromatic source function where .
Using the raymarching 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 (, , and respectively). The brightness temperature is computed for each pixel of each map as per the RayleighJeans approximation
(19) 
where
(20)  
(21) 
and the optical length travelled by a photon is related to the depth along the observer’s line of sight by , where is the sum of the absorption coefficients from different processes affecting the frequency . In these integrals, 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 .
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 quasilongitudinal and quasitransversal approximations.
Using the C7 semiempirical atmosphere (Avrett & Loeser, 2008) and the dipole model discussed in Sec. 3.1 with a looptop magnetic field strength of 100 G the quasilongitudinal approximation holds for outside the range at optical depth from freefree absorption. For cells where falls within the range we simply adjust the viewing angle to the closer edge of the range to retain the validity of the quasilongitudinal 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 400 km into the chromosphere. In regions deeper in the chromosphere, for viewing angles close to perpendicular, the quasilongitudinal 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 and in the frame of the wave
(22) 
In the reference frame of an observer the circular polarisation will be reversed such that .
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 raymarching 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 Raymarching
The standard approach for computing the emission map of a threedimensional volume is to use ray tracing, however this becomes significantly more computationally demanding as the number of discrete volume elements (voxels, assumed quasihomogeneous threedimensional cubic volumes) increases due to the number of rayvoxel 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.
Raymarching 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 raymarching 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 raymarching 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 threedimensional Cartesian grid and stored in a volume texture (threedimensional cuboidal array) before the raymarching 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 raymarching 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 axisaligned 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 raymarch. 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 subdomains. If a ray intersects with one of these higher resolution regions (determined by the intersection of the ray with the subdomain’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 raymarching approach will remain unchanged. The model is currently based around the concept of a single AABB containing the entire scene at lowresolution and a number of refined regions contained strictly within this AABB. The information from the lowresolution AABB is then ignored in the locations where a refined region is present. These AABBs serve as the boundaries of a twolevel threedimensional 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 raymarching code. It would equally be simple to extend the code to support a full multilevel 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 raytracing and raymarching. 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 lookedup wherever the blue point happens to be. When applying a traditional, simple, raytracing technique without additional voxel culling optimisations, the rayvoxel 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 raymarching 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 rayvoxel intersections, for the cases used in Thyr raymarching is typically 1.5 orders of magnitude faster than raytracing, 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 raymarching 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 stepsize (Thyr uses 0.1 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 twodimensional loop is simulated. The loop is observed from two different viewing angles . The specific intensity (sfu^{1}^{1}1solar flux units (1 sfu = Jy = W m Hz = erg s cm Hz)/sterad) is then computed for a telescope beam scanning across the loop with angle . 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 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 and 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 threedimensional simulation performed in Thyr.
The spatial offset of Thyr’s results in the 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 powerlaw nonthermal electron distribution, and the other atmospheric parameters specified by a semiempirical 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 postprocessing of results.
5.1 Semiempirical atmosphere and parameters
In our example simulations we use the C7 quiet Sun semiempirical atmosphere of Avrett & Loeser (2008). This atmosphere is extrapolated to coronal parameters, by extrapolating linearly in loglog 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 nonthermal electron distribution is embedded in the semiempirical atmosphere, with a density defined as cm in the corona, and cm at the top of the chromosphere. This increase of 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 s entering the chromosphere, which is consistent with observations (e.g. Simões & Kontar, 2013). These electrons follow a powerlaw distribution in energy, with a minimum cutoff energy of keV, a maximum cutoff of MeV, and a spectral index of . We assume that these electrons have an isotropic pitchangle distribution, but the code supports arbitrary pitchangle distributions.
As we descend into the chromosphere the nonthermal 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 where is the initial electron energy in keV and is in cm (TandbergHanssen & 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 cm () and a centretocentre footpoint separation of cm (). The radius at the looptop ( in the Kuznetsov et al. (2011) model) is 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 .
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 raymarching algorithm automatically transitions between the low and highresolution 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 highresolution 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 field. The effects of insufficient resolution within the simulation can be seen by comparing the simulations in Fig. 4. The lowresolution simulation of Fig 4a produces an lower brightness temperature than the highresolution 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 of these two simulations can be computed. In practical terms, this calculation is simply the sum of the flux density of each pixel (column and row ), obtained from their specific intensity over the pixel solid angle :
(23) 
The specific intensity of the pixels are directly found from their brightness temperature values via the RayleighJeans law:
(24) 
The resulting spectra are shown in Fig. 6, in sfu. Both of these plots show the characteristic “invertedv” 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 highresolution case, being spread across a much larger voxel in the lowresolution case. Although the spectrum of gyrosynchrotron radiation from a uniform source would present harmonics, they are not visible in the spatiallyintegrated spectra due to the overlap of harmonic structures from a spatiallyvarying field and atmospheric parameters, as has been previously discussed by Klein & Trottet (1984); Simões & Costa (2006). The spectra for the low and highresolution cases behave similarly at low frequencies, but differences emerge at higher frequencies. In the highresolution simulation there is a hardening of the spectrum of the ordinary and extraordinary modes of the gyrosynchrotron radiation at high frequencies. As these nonthermal modes dominate over the thermal emission the hardening in the total spectrum must be nonthermal 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 . 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 threedimensional 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 extraordinarymode emission is seen at less than 2 GHz as the extraordinarymode is unable to propagate due to the ratio of the plasma frequency to the gyrofrequency 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 highfrequency hardening originates.
The results presented here were simulated with an isotropic electron pitchangle 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 highfrequency 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 raymarching is performed in Lua. Plotting is done through our Plyght tool, also available on Github ^{2}^{2}2https://github.com/Goobley/Plyght, that allows any language with a socket API or foreign function interface to easily produce twodimensional plots using matplotlib (Hunter, 2007). The code also supports exporting all data to the widely supported comma separated variable (csv) format to allow for postprocessing and plotting in any language.
7 Discussion
The Thyr tool, its initial validation, and example simulations have now been presented. The highresolution 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 (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 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 smallscale 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 threedimensional 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 userspecified 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 nonthermal 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 freefree 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 highresolution – 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. Threedimensional 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/034066, 04/142485, 08/093392 and 2009/183867. We acknowledge the use of colourblind safe and printfriendly 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.
References
 Alissandrakis & PrekaPapadema (1984) Alissandrakis C. E., PrekaPapadema P., 1984, Astron. Astrophys., 139, 507
 Allred et al. (2015) Allred J. C., Kowalski A. F., Carlsson M., 2015, Astrophys. J., 809, 104
 Avrett & Loeser (2008) Avrett E. H., Loeser R., 2008, Astrophys. J. Suppl. Ser., 175, 229
 Bastian (1998) Bastian T. S., 1998, Annu. Rev. Astron. Astrophys., 36, 293
 Brajša et al. (2017) Brajša R., et al., 2017, Astron. Astrophys., 17
 Carlsson & Stein (1995) Carlsson M., Stein R. F., 1995, Astrophys. J., 440, L29
 Cohen (1960) Cohen M. H., 1960, Astrophys. J., 131, 664
 Costa et al. (2013) Costa J. E. R., Simões P. J. A., Pinto T. S. N., Melnikov V. F., 2013, Publ. Astron. Soc. Japan, 65, 1
 Cuambe et al. (2018) Cuambe V. A., Costa J. E., Simões P. J., 2018, Mon. Not. R. Astron. Soc., 477, 1512
 Dulk (1985) Dulk G. A., 1985, Annu. Rev. Astron. Astrophys., 23, 169
 Fleishman et al. (2015) Fleishman G., Loukitcheva M., Nita G., 2015, in Iono D., Tatematsu K., Wootten A., Testi L., eds, Astronomical Society of the Pacific Conference Series Vol. 499, Revolution in Astronomy with ALMA: The Third Year. p. 351 (arXiv:1506.08395)
 Fleishman et al. (2016) Fleishman G. D., Pal’shin V. D., Meshalkina N., Lysenko A. L., Kashapova L. K., Altyntsev A. T., 2016, ApJ, 822, 71
 Fletcher & Simões (2013) Fletcher L., Simões P. J. A., 2013, NRSO Rep., p. 6
 Fontenla et al. (2009) Fontenla J. M., Curdt W., Haberreiter M., Harder J., Tian H., 2009, ApJ, 707, 482
 Gary et al. (2018) Gary D. E., et al., 2018, Astrophys. J.
 Giménez de Castro et al. (2009) Giménez de Castro C. G., Trottet G., SilvaValio a., Krucker S., Costa J. E. R., Kaufmann P., Correia E., Levato H., 2009, Astron. Astrophys., 507, 433
 Giménez de Castro et al. (2013) Giménez de Castro C. G., Cristiani G. D., Simões P. J., Mandrini C. H., Correia E., Kaufmann P., 2013, Sol. Phys., 284, 541
 Gordovskyy et al. (2017) Gordovskyy M., Browning P., Kontar E., 2017, arXiv Prepr., 116, 1
 Graham et al. (2013) Graham D. R., Hannah I. G., Fletcher L., Milligan R. O., 2013, Astrophys. J., 767, 83
 Guidice (1979) Guidice D. A., 1979, Bull. Am. Astron. Soc., 11, 311
 Heinzel & Avrett (2012) Heinzel P., Avrett E. H., 2012, Sol. Flare Magn. Fields Plasmas, pp 31–44
 Heinzel et al. (2017) Heinzel P., Kleint L., Kasparova J., Krucker S., 2017, Astrophys. J., 48
 Holman (2003) Holman G. D., 2003, Astrophys. J., 586, 606
 Hudson et al. (1994) Hudson H. S., Strong K. T., Dennis B. R., Zarro D., Inda M., Kosugi T., Sakao T., 1994, Astrophys. J., 422, L25
 Hunter (2007) Hunter J. D., 2007, Comput. Sci. Eng., 9, 90
 Iwai et al. (2017) Iwai K., Loukitcheva M., Shimojo M., Solanki S. K., White S. M., 2017, Astrophys. J. Lett. Vol. 841, Issue 2, Artic. id. L20, 8 pp. (2017)., 841, L20
 Kaufmann et al. (2004) Kaufmann P., et al., 2004, Astrophys. J., 603, L121
 Kaufmann et al. (2008) Kaufmann P., et al., 2008, Soc. PhotoOptical Instrum. Eng. Conf. Ser., 7012, 70120L
 Kay & Kajiya (1986) Kay T. L., Kajiya J. T., 1986, ACM SIGGRAPH Comput. Graph., 20, 269
 Klein & Trottet (1984) Klein K. L., Trottet G., 1984, Astron. Astrophys., 141, 67
 Krucker et al. (2013) Krucker S., et al., 2013, Astron. Astrophys. Rev., 21
 Kundu et al. (2001) Kundu M. R., Nindos A., White S. M., Grechnev V. V., 2001, Astrophys. J., 20, 880
 Kuroda et al. (2018) Kuroda N., Gary D. E., Wang H., Fleishman G. D., Nita G. M., Jing J., 2018, Astrophys. J., 852, 32
 Kurucz (1970) Kurucz R. L., 1970, SAO Spec. Rep. 309, p. 286
 Kuznetsov et al. (2011) Kuznetsov A. a., Nita G. M., Fleishman G. D., 2011, Astrophys. J., 87, 1
 Lee & Gary (2000) Lee J., Gary D. E., 2000, Astrophys. J., 543, 457
 Lim et al. (1994) Lim J., Gary D., Hurford G. J., Lemen J. R., 1994, Astrophys. J., 430, 425
 Lüthi et al. (2004) Lüthi T., Magun A., Miller M., 2004, Astron. Astrophys., 415, 1123
 MacGregor et al. (2018) MacGregor M. A., Weinberger A. J., Wilner D. J., Kowalski A. F., Cranmer S. R., 2018, Astrophys. J. Lett., 855, L2
 Mrozek & Tomczak (2004) Mrozek T., Tomczak M., 2004, Astron. Astrophys., 415, 377
 Nakajima et al. (1985) Nakajima H., et al., 1985, Publ. Astron. Soc. Japan, 37, 163
 Nakajima et al. (1994) Nakajima H., et al., 1994, IEEE Proc., 82, 705
 Nindos et al. (2000) Nindos A., White S. M., Kundu M. R., Gary D. E., 2000, Astrophys. J., 533, 1053
 Nita et al. (2015) Nita G. M., Fleishman G. D., Kuznetsov A. A., Kontar E. P., Gary D. E., 2015, Astrophys. J., 799
 Osborne (2018) Osborne C. M. J., 2018, Goobley/Thyr2, doi:10.5281/zenodo.1483653, https://doi.org/10.5281/zenodo.1483653
 Ramaty (1969) Ramaty R., 1969, Astrophys. J., 158
 Simões & Costa (2006) Simões P. J. A., Costa J. E. R., 2006, Astron. Astrophys., 453, 729
 Simões & Costa (2010) Simões P. J. A., Costa J. E. R., 2010, Sol. Phys., 266, 109
 Simões & Kontar (2013) Simões P. J. A., Kontar E. P., 2013, Astron. Astrophys., 551, 10
 Simões et al. (2013) Simões P. J., Fletcher L., Hudson H. S., Russell A. J., 2013, Astrophys. J., 777, 1
 Simões et al. (2015) Simões P. J. A., Hudson H. S., Fletcher L., 2015, Sol. Phys., 290, 3625
 Simões et al. (2017) Simões P. J. A., Kerr G. S., Fletcher L., Hudson H. S., Guillermo Giménez De Castro C., Penn M., 2017, Astron. Astrophys., 605
 Stähli et al. (1989) Stähli M., Gary D. E., Hurford G. J., 1989, Sol. Phys., 120, 351
 Stix (1962) Stix T. H., 1962, The Theory of Plasma Waves
 TandbergHanssen & Emslie (1988) TandbergHanssen E., Emslie A. G., 1988, The Physics of Solar Flares
 Trottet et al. (2012) Trottet G., Raulin J. P., Giménez De Castro G., Lüthi T., Caspi A., Mandrini C. H., Luoni M. L., Kaufmann P., 2012, Energy Storage Release through Sol. Act. Cycle Model. Meet Radio Obs., 9781461444, 33
 Trottet et al. (2015) Trottet G., et al., 2015, Sol. Phys., 290, 2809
 Trulsen & Fejer (1970) Trulsen J., Fejer J. A., 1970, J. Plasma Phys., 4, 825
 Tzatzakis et al. (2008) Tzatzakis V., Nindos A., Alissandrakis C. E., 2008, Sol. Phys., 253, 79
 Valio et al. (2013) Valio A., Kaufmann P., Giménez de Castro C. G., Raulin J. P., Fernandes L. O., Marun A., 2013, Sol. Phys., 283, 651
 Varady et al. (2010) Varady M., Kašparová J., Moravec Z., Heinzel P., Karlický M., 2010, IEEE Trans. Plasma Sci., 38, 2249
 Wang et al. (1994) Wang H., Gary D. E., Lim J., Schwartz R. A., 1994, Astrophys. J., 433, 379
 Wang et al. (1995) Wang H., Gary D. E., Zirin H., Schwartz R. A., Sakao T., Kosugi T., Shibata K., 1995, Astrophys. J., 453, 505
 Wang et al. (2015) Wang Z., Gary D. E., Fleishman G. D., White S. M., 2015, Astrophys. J., 805, 1
 Wedemeyer et al. (2016) Wedemeyer S., et al., 2016, Space Sci. Rev., 200
 White et al. (2017) White S. M., et al., 2017, Sol. Phys., 292
 Wild & Hill (1971) Wild J. P., Hill E. R., 1971, Aust. J. Phys., 24, 43