Tracking solar radio bursts using Bayesian multilateration

Solar radio bursts (SRBs), are emitted by electrons propagating through the corona and interplanetary space. Tracking such bursts is key to understanding the properties of accelerated electrons and radio wave propagation as well as the local plasma environment that they propagate through. Here, we present a novel multilateration algorithm called BayEsian LocaLisation Algorithm (BELLA). In addition, apparent SRB positions from BELLA are compared with comparable localisation methods and the predictions of solar wind models. BELLA uses Bayesian inference to create probabilistic distributions of source positions and their uncertainties. This facilitates the estimation of algorithmic, instrumental, and physical uncertainties in a quantitative manner. We validated BELLA using simulations and a Type III SRB observed by STEREO A/B and Wind. BELLA tracked the Type III source from $\sim$ 10--150 $R_{sun}$ (2-0.15 MHz) along a spiral trajectory. This allowed for an estimate of an apparent solar wind speed of $v_{sw} \sim$ 400 km s$^{-1}$ and a source longitude of $\phi_0 \sim$ 30deg. We compared these results with well-established methods of positioning: Goniopolarimetric (GP), analytical time-difference-of-arrival (TDOA), and Solar radio burst Electron Motion Tracker (SEMP). We found them to be in agreement with the results obtained by BELLA. Additionally, the results aligned with solar wind properties assimilated by the Heliospheric Upwind Extrapolation with time dependence (HUXt) model. We have validated BELLA and used it to identify apparent source positions as well as velocities and densities of the solar wind. Furthermore, we identified higher than expected electron densities, suggesting that the true emission sources were at lower altitudes than those identified by BELLA, an effect that may be due to appreciable scattering of electromagnetic waves by electrons in interplanetary space.


Introduction
Solar radio bursts (SRBs) have been used to track energetic particles from the Sun for decades (e.g.Reiner et al. 1998;Cecconi 2007;Jebaraj et al. 2020).The most abundant of these SRBs are Type IIIs (Reid & Ratcliffe 2014), which are generated by mildly relativistic beams of electrons travelling along open magnetic field lines that then interact with the local plasma and generate Langmuir waves (Lorfing & Reid 2023).These electrostatic waves decay to produce radio wave radiation at the plasma frequency or its second harmonic (Ginzburg & Zhelezniakov 1958;Reid & Ratcliffe 2014).Type III radio bursts have proven to be a productive instrument to study the properties of the solar corona (Reid & Ratcliffe 2014) and the physical mechanisms that drive solar energetic particles (SEP) into the heliosphere (e.g.Reiner et al. 1998).Interferometric telescopes such as the Low-Frequency Array (LOFAR; van Haarlem et al. 2013), which operates at 10-250 MHz, allow for relatively accurate tracking of radio sources at < 3R ⊙ .For example, Morosan et al. (2014) tracked a Type III SRB at 30-90 MHz using LOFAR tied-array observations with high temporal (∼ 30 ms) and spectral (12.5 kHz) resolution.Furthermore, LOFAR facilitates interferometric imaging in the manner carried out by Maguire et al. (2021), where a Type II radio source is imaged to a height of ∼ 0.5 R ⊙ .However, due to the Earth's atmospheric cutoff, lower frequency emission from higher altitudes cannot be observed using ground-based telescopes, and thus space-based antennas are required at these frequencies (≤ 10 MHz).
Observations at less than 10 MHz are particularly well suited to studying radio sources at greater distances from the Sun.With multiple spacecraft, source positions can be triangulated using a variety of methods.A well-established method of triangulation is the goniopolarimetric (GP) method, or direction finding, (Manning & Fainberg 1980) in which the Poynting k vectors of the radio waves observed by different spacecraft are back propagated until they intersect (e.g.Reiner et al. 1998).An example of this method is described in Magdalenić et al. (2014), who used observations from the Waves instrument on board the Wind spacecraft (Wind/Waves; Bougeret et al. 1995) and the Waves instruments on board the Solar Terrestrial Relations Observatory (STEREO/Waves; Bougeret et al. 2008) to triangulate a Type II associated with a coronal mass ejection (CME).Weber et al. (1977); Steinberg et al. (1984) showed earlier time delay measurements with two spacecraft as an alternative to the GP method.With the addition of new spacecraft, the time-difference-of-arrival (TDOA) method was developed (Alcock 2018) for the purpose of tracking radio sources as shown in Badman et al. (2022), where a Type III radio burst was successfully tracked from 0.1-16 MHz using the Parker Solar Probe (PSP/FIELDS; Bale et al. 2016), STEREO A/Waves, and Wind/Waves.An alternative method of positioning using arrival times is the Solar radio burst Electron Motion Tracker (SEMP; Zhang et al. 2019), which fits a Type III to a Parker spiral using a forward modelling method.However, SEMP is constrained by the use of density models and an assumed constant solar wind speed and field line.Another innovative method of positioning is shown in Reiner et al. (2009) and Musset et al. (2021), where they used the intensity of the Type III radio sources in order to calculate the directivity of the radio emissions and density models to obtain the radial distance from the Sun.Lastly, Chen et al. (2023) have combined simulations of density fluctuations with several of these methods to account for scattering of the radio waves as they escape their source.
These localisation methods have various advantages and disadvantages as a result of their assumptions.For example, GP has the advantage that it requires a minimum of two spacecraft in order to obtain a solution as opposed to TDOA, which requires three spacecraft.However, GP is limited, as it requires extensive knowledge of antenna geometry, spacecraft potential, and spacecraft attitude state; for example, while Wind is spinning, STEREO is three-axis stabilised, requiring different techniques to obtain the Poynting vectors (Manning & Fainberg 1980;Krupar et al. 2012).It also has the disadvantage that when a pair of spacecraft is at a large angle of separation (e.g.∼ 180 • ), the Poynting vectors might point at each other, and the method therefore does not yield a solution.An example of this case was noted by Jebaraj et al. (2020), where the STEREO A/Waves-STEREO B/Waves pair was not used due to the unreliability of the results as a consequence of the large angular separation of the spacecraft.In contrast, multilateration works best at large separation angles (see Sect. 3).Badman et al. (2022) showed a longitudinal displacement close to the Sun on the order of the spread in positions from the time resolution error between the results obtained by their LOFAR interferometric results and the multilateration performed using PSP, STEREO A, and Wind.This discrepancy may likely be attributable to a number of factors, physical or systematic.The free-streaming assumption is intrinsic to the TDOA method via the assumption of light travelling at the speed of light (c), and as a consequence, the most relevant source of physical error could be radio scattering (Chen et al. 2023;Kontar et al. 2019).Other sources of errors could be the low temporal resolution of the spacecraft, in particular the Wind/Waves instrument, which provides data at a 1-min cadence; the close proximity of PSP to the Sun at perihelion (∼10 R ⊙ ), or a poor spacecraft configuration (discussed in Sect.3).Badman et al. (2022) estimated errors via repeating the TDOA process a specific number of times with different timestamps selected from an uncertainty range informed by instrument resolution and captured only one of the error sources listed above.A more elegant solution is the application of Bayes' theorem for the purpose of applying a statistical method that would allow for a quantitative analysis of the multiple error sources.
In this paper, we present the BayEsian LocaLisation Algorithm, (BELLA; Canizares 2023)1 a novel method of tracking the apparent positions of SRB sources.In Sect.2, we introduce Bayesian multilateration and describe its implementation.In Sect.3, we perform a simulation for the purposes of validation and choosing a suitable candidate for a use case.In Sect.4, we characterise a Type III radio burst, and finally in Sect.5, we compare the results obtained by BELLA with other well-established positioning methods and a solar wind model.We also discuss some of the potential applications of BELLA for the estimation of solar wind velocities and electron densities along the Parker spiral.

Methods
The chosen method of localisation for this study is multilateration due to its simplicity.Multilateration (also known as trilateration if only three receivers are used) uses timestamps from a number of different receivers to position the source of an emission using geometry.There are two main types of multilateration: Time-of-arrival (TOA) and time-difference-of-arrival, or TDOA, (see Fig. 1).The main difference between these two methods lies in the time difference (∆t) required.The TOA method is the method of multilateration used by global positioning systems such as Galileo and GPS, and it uses a ∆t n0 between the emission time t 0 and the received time t n , where n = 1, 2, 3, ..., depending on the number of receivers.The TDOA method, on the other hand, does not require the time of emission t 0 but only the time received by the different receivers, and thus the ∆t i j required is defined as ∆t i j = t i − t j , where i and j are subscripts for each of the receivers.This time difference (∆t) has implications in the geometrical analysis to detect the source of the emission where TOA is governed by the solution to the parametric equation of n circles with radius and the equation of the circle where (x n , y n ) corresponds to the coordinates of each of the spacecraft receivers.The TDOA method's geometrical approach is slightly different in that it requires solving the solution to the parametric equation of n number of hyperbolas as shown in the appendix of Badman et al. (2022) and available at Badman (2023).The disadvantage of this approach is that hyperbolas contain two branches, and therefore false positives may appear, but this can be mitigated by making further assumptions to choose one solution.Due to the fact that TOA yields a single solution, as opposed to multiple ones; for its mathematical simplicity; and since it can be easily extended to utilise multiple spacecraft, we chose TOA for the purposes of this work.

Bayesian multilateration
Bayesian statistics for multilateration is widely used in different fields, from communications where it is used to track Radiofrequency identification (RFID) nodes (Zhou & Shi 2009;Sanpechuda & Kovavisaruch 2008) to biology, where Reinwald et al. (2021) used pressure waves and seismographs to track the movement of elephants in Africa.Here, we apply the same Bayesian concepts and techniques used by Reinwald et al. (2021) and introduced by Speagle (2019) to track solar radio bursts in the solar corona.
The posterior of a Bayesian statistics positional problem is defined as the probability of finding the position of a source (x) and the speed of signal propagation (v) given an observable ∆t: (3) According to Bayes' theorem, The conditional probability may be written in terms of the priors Rearranging Eq. 5 and substituting into Eq. 4 results in where each term is defined as: -P(x, v | ∆t) is the posterior distribution or the probability of finding the source at position x with a speed of propagation v given an observable ∆t.-P(∆t | x, v) is the likelihood function or the probability of observing a ∆t given a position x and a speed of propagation v.This likelihood function is our physics model, and it is defined by the type of multilateration one wants to apply.In this case, TOA is governed by ∆t = d/v, where d is the distance between the source and the receiver.-P(x | v) is the prior distribution of the source position given a speed of propagation v. -P(v) is the prior distribution of the speed of propagation of the photons.The main assumption made by traditional multilateration methods such as TOA and TDOA is that the speed of propagation of the electromagnetic radiation is c, which is valid asymptotically far from the radio source but may break down close to the emission site where the emitted radiation is at a frequency similar to the local plasma frequency (e.g.Thejappa & MacDowall 2010;Kontar et al. 2019).The Bayesian solver does not make this assumption and allows for a distribution of the speed of propagation.The objective of this method is to accommodate the possibility that the ray from source to observer may have a longer than expected time of arrival due to scattering or refraction near the source but without needing to introduce detailed ray-tracing simulations (e.g.Musset et al. 2021;Chen et al. 2023).This can be modelled as a truncated normal distribution with a maximum of c to restrict all solutions to physical solutions, or it can be modelled as a normal distribution with c as the maximum, allowing for the Markov chain Monte Carlo (MCMC) solver (see Sect. 2.2) to test if the solution is physical or not.
Convergence is expected at v ≈ c or v < c for a physical solution.We note this implies that the path of the light ray is still along a straight line.-P(∆t) is the evidence (i.e.observed) probability shown as a distribution of the time difference between the time of emission t 0 and the time of arrival t 1 .The time of arrival is observed by the spacecraft, and the time of emission is fitted by the MCMC solver, which finds a time of emission that would be consistent for all spacecraft.

Implementation of Bayesian multilateration
The Bayesian multilateration (see Sect. 2.1) was performed with the PyMC (Salvatier et al. 2016;Wiecki et al. 2023) package in Python.The PyMC package allows for Bayesian inference based on probabilistic models and creates an environment where prior distributions can be defined intuitively within a model container that performs Bayesian statistics automatically.BELLA makes use of the normal distribution function embedded in PyMC as its primary source of prior distributions.The probability density function (PDF) of a normal distribution is where σ is the standard deviation and µ is the mean.In order to perform Bayesian inference, a PyMC model context manager was created.The priors were then defined as follows: -The prior for the coordinate space of the source positions P(x | v) was implemented as a stochastic random variable sampled from a normal distribution with µ = 0 and σ = hR ⊙ , where h is an educated guess of the radial distance of the emission to the Sun in the range [-215, 215] R ⊙ .This guess does not affect the posterior distribution, but it does aid in achieving convergence.The term h was chosen to be onefourth of the space domain in one dimension, which is ∼ 80 R ⊙ .-P(v) is the speed of propagation prior defined as a stochastic random variable sampled from a normal distribution with µ = c and σ = 0.1c.-P(∆t) is the probabilistic evidence implemented as a deterministic variable obtained by sampling from a t 0 uniform distribution with limits much larger than the light travel time corresponding to the coordinate space limits (for example, t lims = ±24 x 60 x 60 s) and subtracting t 0 from the expected time of arrival ∆t = d/v, where d is the distance between the location of the spacecraft x sc and the prior for source positions x as d = |x sc − x|, and v is sampled from the speed of propagation prior.
We note that the speed of propagation in the model is not c but is determined by the distribution of propagation speeds P(v).This is one of the advantages of this method, as it no longer assumes that the speed of propagation is v = c.We also emphasise that the prior P(v) was allowed to offer sampling where v > c.This was done for the purpose of testing if the Bayesian solver converges to a propagation speed that is physical (i.e.v ≤ c ).Thus, if the maximum of the posterior distribution shows that v > c, the results are no longer considered physical, and further investigation is necessary.Another benefit of this method is that it is fully compatible with any number of spacecraft as long as n ≥ 3 (and without the need to alter the code).
Once the priors are defined in the PyMC model context manager, the likelihood function P(∆t | x, v) is defined as an observed stochastic variable that samples from a normal distribution with µ = ∆t (calculated earlier) and σ = t cadence (instrument cadence).Then, an observed parameter t toa (the time of arrival as measured by each spacecraft at a given frequency) is used by the PyMC environment for validation.
With the parameter space set up, the PyMC environment was then ready for sampling.The PyMC package uses by default the No-U-Turn Sampler (NUTS; Hoffman et al. 2014), which is a self-tuning sampling method that prevents the sampling from retracing its own steps and thus helps efficiently explore the target distribution.The sampler was set up in the PyMC context manager to perform four different sequences of samples, also known as chains.Each chain was initialised using 2000 tuning samples, and then an additional 2000 samples were drawn to achieve a singular stable solution, also known as convergence.When convergence was achieved, the posterior distributions (P(x, v | ∆t)) were found to be approximate and discrete normal distributions.The mean of the distributions could then be defined as µ = [x x , x y ].This is the most probable source location with 1σ uncertainty for each coordinate.The area of confidence was obtained by generating an ellipse with the axes σ = [∆x x , ∆x y ] as read off the resulting distributions.same data: the left column is the distribution of the sampled values, and the right column is the value of the sample values as a function of index.From top to bottom, we found the posterior distribution of the x coordinate, the posterior distribution of the y coordinate, the distribution of the speed of propagation, and the distribution of the time of emission with respect to the observed time.Throughout this work, we use the peak of these distributions as the most probable location of the source and 1σ as the uncertainty of the distribution.The peak and σ values are highlighted with vertical lines, as seen in Fig. 2 (left).The distribution of the speed of propagation also displays a red line corresponding to the v = c limit.This was used as a quick view warning to test the v ≤ c criterion.

Simulations and method validation
In order to first validate the method and its capabilities, a simulation was set up under the controlled conditions of a known set of source positions and a known propagation speed (i.e.no refractive or scattering effects).In these simulations, we recovered the ground truth signal to a better or worse extent depending on the spacecraft configuration.This allowed us to test for the contribution of the spacecraft configuration to the uncertainty of the localisation results.Figure 3 and Fig. 4 show the results of these simulations, with the first figure showing four different instances of the simulation and the latter providing a summary.In order to see the effects of instrumental cadence on the performance of the simulations, all cases were simulated twice: once at a cadence of 30 s (blue) and once at a cadence of 60 s (red).The background maps were all simulated at 60-s cadence.In addition to validating the methodology, these simulations allowed for the selection of a suitable Type III SRB candidate by discarding spacecraft configurations that would guarantee uncertain results.The simulation was composed of the following: A three-spacecraft system analogous to the real STEREO A/B and Wind configuration was established.For simplicity, the Ahead and Behind spacecraft were kept at a constant distance of ±1% of 1 au, respectively.The Wind spacecraft was kept stationary at L1, while the Ahead and Behind spacecraft were synchronously separated from Wind in ±10 • increments, respectively.In order to investigate the effect of spacecraft geometry and instrument time resolution on the inherent uncertainty of the method, a mesh-grid of 500 x 500 nodes was established, and each node This simulation was performed assuming different cadences, which also had an effect on the performance of the multilateration, but it was not as dominant as the angle of separation.These simulations suggest that BELLA will have a minimal contribution to the uncertainty of the results if the angle of separation is in the range of ±[45 • , 165 • ].This range is defined as the 'sweet spot' region and is highlighted in green.See Fig. 3 for frames of this simulation.
was treated as an emission source.By assuming free streaming and a propagation speed of v = c, synthetic arrival times were calculated as t = d/c for each source-spacecraft pair.A small amount of Gaussian noise (on the order of <1%) was applied to the data.If there were regions where the uncertainty was above a user-specified threshold despite applying the smallest amount of noise, then the configuration was rejected.The user-specified threshold was 25 R ⊙ , as this is the light travel distance for an instrument with a 60-s cadence (typical instrument time resolution).The noisy data was used as an input for the Bayesian procedure described in Sect.2.2.Posterior distributions for the detected position at each node gave corresponding 1σ uncertainties in the x and y directions.The maximum value of the x and y uncertainties was taken as the value for the uncertainty map colour.As shown in Fig. 3, it is clear that when two or more spacecraft are close together, a source anywhere in the ecliptic plane has a large uncertainty (yellow regions), but when they are equidistantly spaced, the uncertainty is minimised (blue regions).A test burst was also simulated by selecting a number of simulated sources along a notional Parker spiral and performing the same procedure at two different cadences (60 s and 30 s).By recovering this known ground truth, we validated the method, and the validation allowed for testing of the effects of the in-strument cadence on the performance of the multilateration.The test burst also allowed for a visual depiction of the discrepancy between the multilaterated source (shown as red and blue stars) and its true counterpart (orange circles) as well as showing the difference in the x and y uncertainties.Figure 4 shows the average uncertainties for the full test burst in solid lines, and the maximum and minimum uncertainties of the nodes are shown as dashed lines.The 'sweet spot' regions are highlighted in green and are the regions where both the x and y uncertainties are under a user-defined threshold.Moreover, these regions served as a starting point for finding suitable candidates.The threshold in Fig. 4 was defined as the largest of the midpoints of the decay and rise in uncertainty regions of the 60-s cadence simulation.This corresponds to the ±[45 • , 165 • ] bands.Figure 4 also shows that for the ideal case of nearly equidistant spacecraft, the ground truth is retrieved accurately, meaning at the expected uncertainty due to the light travel time corresponding to the cadence.However, as n number of spacecraft move closer together, they start to behave as a unique receiver.This is particularly evident for the 180 • case, where the Ahead and Behind spacecraft are aligned with Wind.In this particular case, the system becomes one dimensional, and the X coordinate uncertainty remains stable as a minimum, while the Y coordinate becomes unstable, and the uncertainty becomes a maximum.The reason for this is that in a one-dimensional system, only two spacecraft are needed to multilaterate a source, and therefore the X coordinate still meets the minimum requirements for multilateration.

Observations
To validate our method and test it on a real event, we examined in detail an individual Type III burst that was clearly observed by multiple well-separated spacecraft (see Fig. 5).The event chosen for this study occurred on 7 June 2012 from 19:30 This makes this region a highly plausible Type III source region given the proximity of open field lines to active region loops.As we show, BELLA allows for a strong inference to back up this hypothesis for our studied burst.
The Type III radio burst was observed by Wind from 0.06 to 13.825 MHz, by STEREO A from 0.15 to 16 MHz, and by STEREO B at a range similar to that of STEREO A but with signs of solar disk occultation, as the higher frequencies of the burst show a decay in brightness.This particular event was chosen for a number of reasons.The location of the spacecraft is optimal for multilateration techniques.The STEREO spacecraft are at ±116 • from the Wind spacecraft, which resulted in nearequidistant separation between all three spacecraft.At this separation, Sect. 3 showed that the intrinsic uncertainty from BELLA is minimal and expected to be dominated by the cadence of the instrument and physical errors.The Type III radio burst is isolated.There is no other bursts or activity that could potentially cause confusion or uncertainty to the edge of the burst.The burst is bright, and the edges are well defined.The low cadence of the The event was chosen due to its relative brightness with respect to the background.The burst was observed to be an isolated event and lasted approximately 20 mins.The frequency range of the burst located by each of the spacecraft was from 13 MHz to 0.1 MHz; however, STEREO B showed a poor S/N at the higher frequency range, presumably due to solar disk occultation.The white line is the fit of the Type IIIs used for the multilateration in the range of 2 MHz to 0.15MHz.
instruments was a limiting factor, and therefore, the burst picked for the validation of BELLA was chosen upon visual inspection to be as clean as possible and with as high S/N as it was possible.All spacecraft are at ∼1 au from the Sun, simplifying the system analysed.It is yet unknown whether the effects of spacecraft distances from the Sun has an effect on the performance of the algorithm.As it is beyond the scope of this initial test of BELLA, we will explore how spacecraft heights, and other properties, can affect the performance of the algorithm in future studies.
The leading edge of the Type IIIs was chosen as the determining feature to obtain time measurements for the multilateration, as it corresponds to the shortest possible path of the radiation from the source to the receiver.Determining a consistent point for the rising time of a Type III light curve for a particular frequency poses a number of challenges when the time resolution of an instrument is as low as 60 s, which is the case of the Wind/Waves instrument, and 38.05 s, as in the case of STEREO/Waves.In order to obtain a consistent rise time for the light curve of a Type III, we fit a Gaussian-Hermite (GH) model (Van Der Marel & Franx 1993) from the Python package Kapteyn (Terlouw & Vogelaar 2014).The GH model is an asymmetric Gaussian-like function that derives values for the skewness and kurtosis of the PDF. Figure 7 shows different examples of the result of fitting a GH model using the kmpfit module of the package for a number of different light curves at varying frequencies.Figure 7 (a) is an example of a low-frequency light curve, and Fig. 7 (b) is an example of a high-frequency light curve.Figure 7 (c) is an example of a light curve with a poor S/N.The criteria for the rise time of the light curve was defined as 1σ to the left of the skewed distribution obtained from the GH fitting parameters.
Time stamps for the multilateration were extracted from the exact same frequencies.In order to account for the difference in frequency channels between STEREO/Waves and Wind/Waves as well as the low temporal resolution of these instruments, we fit a time evolution function to the rise times detected using the GH method.This time evolution function was parameterised using a polynomial of the form where t is the rise time; f is the frequency; and a 2 , a 1 , and a 0 are coefficients that represent the curvature, the drift rate, and initial position of the burst, respectively.These fits can be observed in Fig. 5 as white lines and show the data extracted for the multilateration.The white dashed lines represent the edges from where the MCMC sampled for the time prior of the Bayesian calculation (see Sect. 2.2).In addition, we directly interpolated individual rise times in order to assess the impact of the time evolution fitting (see Fig. 5).

Results and discussion
The results of BELLA are shown in Fig. 8.The background map was obtained by simulating emission from every pixel in the map, and it displays the contribution of BELLA to the uncertainty of the results.The blue regions of the background map are areas where BELLA has a minimal effect in the uncertainty of the output positions and the yellow regions are regions where BELLA cannot accurately position a source, giving results similar to Fig. 3 (d).The large trianguloid shape of the blue region is a direct consequence of the optimal position of the spacecraft for that particular date.The spacecraft were equally spaced, and therefore a symmetrical output was expected.Simulations of a similar burst with varying spacecraft positions shown in Figs. 3  and 4 show that the contribution of BELLA to the uncertainty of the measurements due to the position of the spacecraft is negligible for the position of the Ahead and Behind spacecraft during this event.
The map was generated using simulated data with a cadence of 38 s, which corresponds to a light travel time uncertainty of ∼ 15 R ⊙ .This simulation differs from that in Fig. 3, as the spacecraft were chosen directly from their heliographic locations on the date and time of the burst.From the map, we observed that the whole burst occurred in a region of < 15 R ⊙ in uncertainty.This means that BELLA has a negligible contribution to the uncertainty of the results, and any source of uncertainty is dominated by the cadence of the instruments and other (mostly physical) sources of uncertainty.This is consistent with the dynamic spacecraft simulation shown in Fig. 3. Overlaid on this map is the BELLA output from the TOA multilateration performed on the Type III radio burst characterisation from Fig. 5.The burst is colour-coded by the frequency of emission (right-hand colour bar).The ellipses around each of the source points are the 1σ uncertainties obtained from the BELLA posterior distributions.Also overlaid on this map is the output from the TDOA multilateration using the method from Badman et al. (2022), and shown in red is the output from the Solar ra- dio burst Electron Motion Tracker (SEMP, Zhang et al. 2019) at 400 km s −1 ± 20 km s −1 solar wind speed.The three methods of positioning were found to be in agreement with each other within 1σ uncertainty.We observed that the TDOA method yielded results nearly identical to the TOA method used by BELLA.This is in agreement with Kaune (2012), which showed that Monte Carlo simulations of TOA and TDOA have no discernible results in accuracy.The results indicate that the TOA solution con-  2022) and SEMP have also been included for the purpose of validating the results obtained by BELLA.Parker spirals show the results from fitted BELLA points (black) and SEMP (red) at 400 km s −1 , with dashed lines indicating the uncertainty at ±20 km s −1 .We found that there are some systematic differences between TDOA/BELLA (ϕ 0 ∼ 30 • ) and SEMP (ϕ 0 ∼ 20 • ), but they are consistent within the error bars that BELLA allows to be quantified.verged to the same solution that TDOA produces geometrically and that the posteriori signal propagation distribution is close to v = c.We note that the line of sight of the radio waves propagation is for the most part unobstructed by the lower layers of the solar corona, which makes this burst ideal for the assumption of v = c.The purpose of this study is to validate BELLA as a method for positioning the source of observed radio bursts.This convergence yields two possibilities: a) the effects of refraction and scattering can be neglected for this particular burst or b), in this case, BELLA has converged on an 'apparent' position of the burst shifted by propagation effects (Chen et al. 2023;Kontar et al. 2023).The effects of refraction and scattering will be the subject of further investigation, but as discussed later in this section, we find evidence that (b) is more likely for this burst.Additionally, Level 3 (L3) GP inversions from STEREO data available at Krupar et al. (2022) are shown in Fig. 9.The GP triangulation was performed using the time of peak flux of the Type III as opposed to the leading edge because the GP inversions at the leading edge are not stable.BELLA, on the other hand, is designed to utilise the leading edge of the Type III, making the two independent methods complementary to each other.A clear trend was observed between the results obtained by the triangulation and by the multilateration, with only a small number of points lying outside of the 1σ area of uncertainty.The agreement in orientation and frequency distribution along the spatial domain serves to validate BELLA as a suitable method of localisation.
Figure 8 also shows the results from SEMP, which assumed a pre-specified constant solar wind speed along the Parker spiral.
The SEMP longitude ϕ 0 parameter was found to be ∼ ϕ 0 = 20 • .A Parker spiral, shown as a black dotted line in the figure, was also overlaid on the map.The parameters for this spiral were obtained by assuming a Parker spiral to be an Archimedean spiral: where r is the distance from Sun, r 0 is the height at which the spiral starts, v sw is the velocity of the solar wind, Ω ⊙ is the angular velocity of the Sun, ϕ is the longitude, and ϕ 0 is the longitude at the base of the spiral.Equation 9shows that an Archimedean spiral is linear in polar coordinates.Using the slope and intercept of the linearised sources in polar coordinates shown in Fig. 10, we derived estimates of the solar wind speed and longitude as v sw = 400 ± 20 km s −1 and ϕ 0 = 30 • ± 2 • .This implies that there was only a discrepancy of ∼ 10 • between the two methods.Comparing these results with Fig. 6, we observed a number of active regions in the 20 • to 30 • longitude, suggesting that these ARs are the source of the emission.Despite this uncertainty, both methods still suggest the coronal hole and active region boundary as the likely Type III source region.Closed active region loops may have access to neighbouring open field lines in order to give injected electron beams access to the solar wind.Figure 8 (c) also shows the results of performing the Bayesian multilateration of the Type III SRB data obtained using the GH algorithm but without the time evolution fit (Eq.8).The difference between Figs. 8 (b) and (c) show that the fitting of Eq. 8 as a time evolution function to the Type III SRBs is responsible for the smooth Parker spiral seen in Fig. 8.This is consistent with the literature that shows that Type III exciters follow a Parker spiral (Reiner et al. 1998;Reiner et al. 2009), and therefore Eq. 8 is presented as a characterisation of the To further verify the plausibility of the inferred trajectory in the context of a realistic solar wind, Fig. 11 shows the BELLA multilateration results overlaid on a solar wind model from the Heliospheric Upwind Extrapolation with time dependence (HUXt; Owens et al. 2020;Barnard & Owens 2022).This model solution was obtained by assimilating the in situ solar wind speeds observed by ACE, STEREO A, and STEREO B to provide a data-constrained solution to the solar wind structure (Lang et al. 2021;Lang & Owens 2019).This data assimilation was performed over a 27-day window centred on the time of the Type III burst.
The Parker spiral obtained from the BELLA multilateration showed agreement with the HUXt output, suggesting that the electron beam of the emission follows a Parker spiral with a solar wind velocity of 480 km s −1 .The HUXt output shows a region of fast solar wind caused by the coronal hole CH1 shown in Fig. 6.This is consistent with the location of the cluster of active regions in Fig. 6 (b), which is adjacent to the eastern region of the CH1 coronal hole region.This further suggests the interface between these active regions and the coronal hole are the origin of the electron beam.
Lastly, Fig. 12 shows the results from plotting the radial distance from the Sun obtained from the BELLA posterior distributions with respect to the frequency of the emission.These frequencies are compared with the plasma frequencies obtained from density models Saito et al. (1977) and Leblanc et al. (1998) at 1 au and a 1/R 2 fundamental and harmonic projection of STEREO A in situ electron density data.We observed that the slope of the BELLA electron density follows an expected 1/R 2 trend.However, the densities obtained by BELLA compared to the different density models were seen to be significantly larger than expected for either the fundamental or harmonic emission extrapolated inwards from the 1-au density measurements and from the typical density models.With the STEREO A in situ data, which were collected over a Carrington rotation period (28 days) and centred on the time of the burst, we could state that the observed source locations are consistent with harmonic radiation along a path at the upper end of the distribution of 1-au measured densities for this interval, but this would require the trajectory to lie along a 'much more dense than usual' region, which is statistically unlikely given that studies such as Steinberg et al. (1984) have suggested the opposite.This is, however, consistent with Chen et al. (2023) and Kontar et al. (2023), which state that radio wave scattering has an effect on the apparent position of SRB exciters, making them appear to be further than their true position.This would imply that the distribution of the speed of propagation v converging at v = c is symptomatic of Type III SRB multilateration being susceptible to converging on the apparent location of the source as lensed by scattering effects.This effect will be the subject of further investigation.

Conclusions
We have presented BELLA, a novel method of positioning Type III radio burst exciters that uses Bayesian statistics as a gateway to perform a detailed error analysis of the uncertainties associated with traditional localisation methods.In order to validate BELLA, a simulation was set up under the controlled conditions Fig. 12. Plasma frequency as a function of heliocentric distance observed by the BELLA multilateration of the Type III SRB and compared against a 1/R 2 projection of STEREO A in situ electron density data (presented as histograms).Also shown are the harmonic and fundamental values given by the Saito and LeBlanc density models at 1 au.BELLA electron densities are higher than the expected electron densities obtained from models.The highest outliers of the harmonic emission derived from STEREO A in situ data are consistent with BELLA multilaterated sources.An alternative reason for these higher than expected electron densities is scattering (Chen et al. 2023;Kontar et al. 2023).
of a known set of source positions and known propagation speeds.BELLA later recovered the ground truth positions to a better or worse extent depending on spacecraft configuration.We also validated BELLA against a real event, and therefore a Type III radio burst detected by Wind and the STEREO pair was analysed and compared with three other methods of positioning (TDOA Badman et al. (2022); GP Krupar et al. (2012); and SEMP Zhang et al. (2019)) as well as with the HUXt solar wind model (Owens et al. 2020;Barnard & Owens 2022), which we constrained by assimilating the available in situ observations.In order to prepare the data for the Bayesian multilateration, we fitted a GH polynomial to each of the light curves using 1σ as the rise time detection trigger.These rise time positions were then used to obtain a time evolution function that allowed us to perform the BELLA multilateration despite instrumental discrepancies in cadence, frequency, and resolution as well as differences in frequency channels between Wind/Waves and STEREO/Waves.Despite these instrumental differences, the three spacecraft used in this study have relatively similar spectrogram and orbital specifications.We focused on the highly symmetric case of the STEREO A/B and Wind constellation.Future work will include a less idealised constellation placement and the use of PSP and Solar Orbiter (SolO; Marsch et al. 2005;Müller et al. 2020) data as well as using more than three receivers to obtain full 3D information.
In summary, the results obtained from this study are the following: 1. Simulations showed that BELLA successfully recovers ground truth data when the spacecraft are evenly separated.We showed that these simulations are capable of distinguishing between regions where the innate multilateration uncertainties are dominant and regions where these innate uncertainties are under the expected instrumental uncertainties.2. We localised a Type III radio burst exciter.Results from the BELLA posterior distributions were found to be in agree-ment with the output of the analytical TDOA multilateration.BELLA currently uses TOA positioning, and we showed that, in our case, TOA converged on the same apparent source location that TDOA retrieves.This agrees with Kaune (2012), suggesting that there is no difference in performance between TOA and TDOA.BELLA was also shown to be in agreement with STEREO L3 GP inversions, showing a clear trend between the two and with SEMP, as the localisation of the SEMP method only showed a difference of ∼ 10 • from the BELLA Parker spiral.3. BELLA does not make the v = c assumption a priori, but we found that the MCMC sampling converged on v ≈ c.This result implies one of two possible outcomes: a) The effects of refraction and scattering did not play a major role in this event, or b) only the apparent location of the burst source can be multilaterated.Chen et al. (2023) showed that the latter is likely the case, and scattering will therefore be the subject of further investigation.4. The morphology of the Parker spiral was analysed for the purpose of obtaining an evolution of the solar wind speeds.We found that the multilaterated positions follow an Archimedean Parker spiral of v sw = 400 ± 20 km s −1 and ϕ 0 = 30 • ± 2 • , which is consistent with the locations of a Coronal hole and the active regions observed on the same date of the event.These source regions were at low latitudes, therefore meaning our 2D treatment of the source location was appropriate.5. Output from the HUXt wind model was also compared to the result obtained by BELLA, showing that the BELLA points put the bursts in a solar wind speed of around 480 km s −1 .However, the curvature of the track and the model Parker spiral suggest this may have been slightly higher than the true solar wind speed at that location.This v sw is 80 km s −1 higher than the results obtained by fitting an Archimedean Parker spiral to the points obtained by BELLA.6. Density models were compared with the plasma density obtained from the frequencies of the Type III SRB and the radial distance from the Sun.We found that the slope of the linearised density model was in agreement with density models from the literature, following a roughly 1/R 2 heliospheric trend.The BELLA sources were found to be located at a significantly higher altitude than implied by density models or in situ density measurements at 1 au.This could be explained either by propagation along an anomalously high density filament and harmonic emission or by the localised source being an 'apparent' source shifted outwards in radius via scattering (Chen et al. 2023;Kontar et al. 2023).Assuming the latter conclusion, it is interesting to consider the implication that the trajectory still follows a feasible Parker spiral and 1/R 2 density trend.This may suggest that the lensing effect of scattering may produce a self-similar transformation on the source trajectory and would suggest the inverse transformation may not require a full ray-tracing simulation (Kontar et al. 2023) and could instead be much simpler.
BELLA is a novel method of multilateration that provides an in-depth analysis of the localisation uncertainties, allowing the user to distinguish between the innate multilateration uncertainties and the instrumental and physical uncertainties.BELLA is available as an open source tool for those that want to use it.It is fully capable of operating with more than three spacecraft in 2D.The 3D capabilities of BELLA still have not been developed, and they will be the subject of further work but are a relatively simple addition to the framework outlined in this paper.This will be of special interest when SolO has reached its high orbital inclination θ ∼ 30 • (Müller et al. 2020), resulting in a moderate angular separation between receivers and consequently reducing the uncertainty in the Z coordinate.

Fig. 1 .
Fig. 1.Time-of-arrival (TOA) method and time-difference-of-arrival (TDOA) method.(a) The TOA method uses the time difference between the time of emission and the time of arrival to obtain parametric equations of a circle, and the location of the source is the point of intersection.(b) The TDOA method uses the time difference between each pair of receivers to obtain parametric equations of a hyperbola, and the point of intersection of the hyperbolas is the location of the source.The uncertainty of a single measurement, shown by the thickness of a path, is governed by the cadence of the instrument.Spacecraft configuration determines how these areas overlap.If the overlapping is perpendicular, the area of uncertainty is minimised.If the overlapping is tangential, the area of uncertainty is maximised.

Fig. 2 .
Fig. 2. Example trace plots of the posterior distributions from BELLA at 0.154 MHz.Sampling chains are shown in the right column, while the corresponding distributions are shown in the left column.Vertical lines on the distributions show the peak and 1σ range of a combined distribution.Panels (a) and (b) are X and Y coordinate distributions, respectively, where the peaks are the coordinates of the source of the emission and 1σ shows the area of confidence.Panel (c) is the speed of propagation of the emission, v.The speed of propagation was allowed to be greater than c for the purposes of testing if the distribution converged at v = c or v < c.A red line is displayed as a quick view of the v = c limit.Panel (d) shows the time of emission t0 of the burst in seconds from midnight.The peak of the t0 distribution shows the most likely time of the emission, and 1σ shows the region of uncertainty for this time of emission.
Figure 2 shows an example of the trace plots generated after one multilateration.Each of the colours in the figure represents a different chain.The columns in the figure represent the exact

Fig. 3 .
Fig. 3. Simulated radio sources and their corresponding BELLA detection at (a) ±10• , (b) ±40• , (c) ±90• , and (d)±170• of Ahead-Behind separation with respect to Wind.Orange dots are the true location of the sources, and stars and ellipses are their corresponding detection and uncertainties, respectively.At each spacecraft configuration, simulations were performed twice with different cadences: 30 s (blue) and 60 s (red).The background maps were generated at a 60-second cadence and show a general overview of the multilateration uncertainty.Uncertainty of the measurements is dominated by the cadence of the instruments and the location of the spacecraft.The best performance of the multilateration technique occurs when the separation angles between all spacecraft are maximised and widely distributed.

Fig. 4 .
Fig. 4. Simulation showing the expected uncertainties when performing multilateration on a STEREO and Wind system.The solid line represents the average 1σ uncertainty of the sources in the simulated burst, and the dashed line represents the maximum and minimum 1σ uncertainties of the positions of the simulated bursts.The angle of separation of the Ahead and Behind spacecraft of the STEREO system with respect to Wind shows a strong contribution on the uncertainty of the localisation method when two or more spacecraft are close to each other.This simulation was performed assuming different cadences, which also had an effect on the performance of the multilateration, but it was not as dominant as the angle of separation.These simulations suggest that BELLA will have a minimal contribution to the uncertainty of the results if the angle of separation is in the range of ±[45 • , 165 • ].This range is defined as the 'sweet spot' region and is highlighted in green.See Fig.3for frames of this simulation.
UT and was observed by Wind (Wind/Waves, Bougeret et al. 1995), STEREO A, and STEREO B (STEREO/Waves, Bougeret et al. 2008)).The datasets were obtained from Bougeret et al. (2021) and Krupar et al. (2022).During the date of the event, the coronal hole identification via a multi-thermal emission recognition algorithm (CHIMERA, Garton et al. 2018) reported a large coronal hole in the north-western quadrant of the Sun spanning ∼ 30 • , labelled as CH1 in Fig. 6 (a). Figure 6 (b) shows over 10 active regions as observed by the SOLIS Vector Spectro Magnetograph(SOLIS VSM; Henney et al. 2006)) and labelled by the National Oceanic and Atmospheric Administration (NOAA).

Fig. 5 .
Fig. 5. Type III radio burst observed by Wind (left), STEREO A (middle), and STEREO B (right) on 7 June 2012.The event was chosen due to its relative brightness with respect to the background.The burst was observed to be an isolated event and lasted approximately 20 mins.The frequency range of the burst located by each of the spacecraft was from 13 MHz to 0.1 MHz; however, STEREO B showed a poor S/N at the higher frequency range, presumably due to solar disk occultation.The white line is the fit of the Type IIIs used for the multilateration in the range of 2 MHz to 0.15MHz.

Fig. 6 .
Fig. 6.Coronal holes on 7 June 2012 obtained by CHIMERA.Fast solar wind regions observed in Fig. 11 are consistent with the location of CH1.Magnetogram data from SOLIS show 10+ as the number of NOAA active regions, four of which are adjacent to the East of CH1.This is later shown to be consistent with the location of the triangulated Type III beam.

Fig. 7 .
Fig. 7. Edge detection from the GH Type III detection algorithm.Blue dots show the light curve data, the yellow dashed line shows the background level derived from the GH Type III fitter, the cyan dashed line shows the GH fitting, and the red star and dashed line show the detection point used to obtain the rise time of the burst.a) Example of a high S/N light curve fit using GH polynomials.b) Example of a fit with a low signal data count.The low data count had no noticeable effect on the performance of the algorithm to obtain the rise time of the light curve.c) Example of a light curve with a poor S/N.The GH algorithm was capable of fitting a Type III light curve profile despite the poor S/N.The detection was observed to be on the order of one cadence from the rise time location based on visual confirmation of the data points in the dynamic spectra plot in Fig. 5.

Fig. 8 .
Fig. 8. BELLA multilateration of the Type III SRB.(a) Top view of the ecliptic plane showing the location of the spacecraft and the multilaterated burst.A background map generated by simulated data tests the performance of BELLA for this particular spacecraft configuration at the location of the sources; (b) zoom-in of the burst; (c) same zoom-in, but the multilateration was performed on the output of the GH polynomial without the final fit of the Type III.We found that BELLA showed results with a confidence area of ∼ 15 − 20 R ⊙ in diameter.The TDOA results from the method in Badman et al. (2022) and SEMP have also been included for the purpose of validating the results obtained by BELLA.Parker spirals show the results from fitted BELLA points (black) and SEMP (red) at 400 km s −1 , with dashed lines indicating the uncertainty at ±20 km s −1 .We found that there are some systematic differences between TDOA/BELLA (ϕ 0 ∼ 30 • ) and SEMP (ϕ 0 ∼ 20 • ), but they are consistent within the error bars that BELLA allows to be quantified.

Fig. 9 .
Fig. 9. STEREO A/B L3 GP inversions(Krupar et al. 2022) overlaid on Fig. 8 (a).The GP points show a clear trend in agreement with BELLA in both orientation and frequency distribution (as indicated by the colour scales).The GP results have a larger spread and make interpretation of an underlying spiral more difficult.All localisation methods are shown to be in agreement, illustrating that complementary but independent methods give approximately the same answer, which further validates the BELLA methodology.

Fig. 10 .
Fig. 10.Localisation of the Type III radio burst in polar coordinates.An ideal Parker spiral is shown as a linear function in polar coordinates where the slope is proportional to the speed of the solar wind.A linear piecewise fit was performed in order to obtain the evolution of the wind speed along the source of the emission.

Fig. 11 .
Fig. 11.Results from the HUXt solar wind model compared with the sources localised by BELLA.(a) Top down ecliptic view of the multilaterated burst and HUXt results.(b) Zoom-in of the region of interest.The HUXt model shows a region of fast solar wind (v sw > 700 km s −1 ) on the +/+ quadrant attributed to the coronal hole CH1 (see Fig. 6 (a)).East of this high speed region is the location of the multilaterated SRB.We found that the SRB sources follow the 480 km s −1 field line with precision (< 5 R ⊙ ) despite the fact that the HUXt model and the BELLA multilateration are independent of each other.The location of the sources with respect to the region of high solar wind speed is consistent with the location of a cluster of active regions seen in Fig. 6 (b).