Investigating Extreme Scattering Events by Volumetric Ray-tracing

Extreme scattering events (ESEs) are observed as dramatic (>50%) drops in flux density that occur over an extended period of weeks to months. Discrete plasma lensing structures are theorized to scatter the radio waves produced by distant sources such as pulsars, causing the signature decrease in flux density and characteristic caustic spikes in ESE light curves. While plasma lens models in the extant literature have reproduced key features of ESE light curves, they have all faced the problem of being highly overdense and overpressured relative to the surrounding interstellar medium by orders of magnitude. We model ESEs by numerically ray tracing through analytic, volumetric plasma lens models by solving the eikonal equation. Delaunay triangulation connecting the rays approximates the wave front, generating a mapping from the observer plane to the source plane to account for multiple imaging. This eikonal method of ray tracing is tested against known analytic solutions and is then applied to a three-dimensional Gaussian-distributed electron volume density lens and a filament model inspired by Grafton et al. We find convergence of our numerical results with established analytic solutions, validating our numerical method, and reproduce ESE-like light curves. Our numerical ray-tracing method lends itself well to exploring the lensing effects of volumetric turbulence as well as sheet-like lenses, which is currently in progress.


INTRODUCTION
Distant, highly energetic objects such as pulsars and quasars are luminous and emit strongly in the radio regime.Typically, the mechanism by which the radio radiation arises comes from synchrotron radiation.Pulsars and quasars possess extremely strong magnetic fields; ∼ 10 9 G for pulsars, and up to ∼ 10 6 G in quasar jets.These strong magnetic fields accelerate free electrons, which emit radiation in the radio regime.Radiation is influenced by free electrons in plasma as it propagates through space.The influence on the propagation of radiation at frequency f through cold plasma is characterized by an index of refraction (Lorimer & Kramer 2004) (1) The plasma frequency is where e is the elementary charge, n e is the electron number density, and m e is the electron mass.
The phase velocity is modified by n as v = c/n where c is the speed of light in vacuum (n = 1).In the interstellar medium (ISM), generally f p ≪ f and so n ≈ 1 − f 2 p /2f 2 .Since n ≲ 1 within cold plasma, the phase velocity exceeds c.
In the early 80's Fiedler et al. (1987) measured the light curves of 36 extragalactic sources with the Green Bank interferometer at 2.7 and 8.1 GHz noting, "unusual minima" (Fiedler et al. 1987) in these sources.The most dramatic dimming was seen in both the 2.7 GHz and 8.1 GHz light curves of quasar QSO (quasi-stellar object) 0954+658.
This extreme scattering event (ESE) occurred over a period of approximately 80 days as opposed to more typical variations for this source, which occurred on the timescale of days.The ESE was characterized by its prolonged flux density minimum of approximately 60 days, which were flanked by two dramatic caustic spikes (focusing of light rays resulting in multiple imaging) at the beginning and end of the event each spanning approximately 15 days each.Fiedler et al. (1987) ruled out the possibility that this event was due to intrinsic variability since models based on an evolving emitter show flux density variations at a lower frequency lagging those at higher frequencies.However, this ESE showed simultaneous variations in 2.7 GHz and 8.1 GHz, atypical of its intrinsic flux variations which lag by 50 days between the two wavebands, which suggested it was due to an occulter of astronomical unit (AU) size at kiloparsec (kpc) distances.Variations and spikes in the light curve were thought to be caused by irregularities in ionised gas density and substructure.From the flux density profile of the ESE, assuming a spherical plasma lens implied the occulter had electron densities n e ∼ 4 × 10 4 cm −3 , much larger than the diffuse ISM.These small plasma lenses are a type of tinyscale ionised structure (TSIS).Such small-scale structures in the ionised ISM may reach sizes as small as 70 -100 km (Rickett 1977;Armstrong et al. 1995

, for example).
There have since been a handful of ESEs observed (Stanimirović & Zweibel 2018).It should be noted that the exact definition of an ESE is unclear; the precise limiting characteristics of an ESE have not been mutually agreed upon between authors.Collectively however, it seems ESEs are characterized by their significant drop (> 50%) in radio (∼ 1 GHz) flux density over the course of weeks to months.Such dramatic flux density variation and long duration distinguishes itself from typical scintillation due to inhomogeneities of the electron density in the interstellar plasma which varies only a few percent over hours or days.
An early, well-known model of a plasma lens was explored by Clegg et al. (1998).The Clegg et al. (1998) lens was modeled by a Gaussian-distributed electron column density.Clegg et al. (1998) found their model was able to produce some notable features of ESE light curves including a defocused region (the minima), and caustic spikes (Clegg et al. 1998, Figure 2).The refractive nature of plasma lenses are described by Equation 1 when the propagating wavelength is much smaller than the spatial variation of n e .The radiation is refracted toward regions of lower phase speed and lower n e .Multiple imaging was predicted by the model and later confirmed by an observation of QSO 2023+335 by Pushkarev et al. (2013).By considering the Fiedler et al. (1987) observations of QSO 0954+658, Clegg et al. (1998) found the plasma lens was 0.38 AU in size with electron density 10 5 cm −3 .This corresponds to pressures of approximately 10 6 − 10 9 K cm −3 , ". . .well in excess of the average ISM pressure of roughly 4000 K cm −3 " (Kulkarni & Heiles 1988).Such highly over-pressured structures suggest they should be transient in nature or would otherwise need to be situated within a high pressure environment.
Geometrical arguments have been proposed to explain the excessively large inferred density and pressures of plasma lenses.Not long after Fiedler et al.'s (1987) ESE discovery, Romani et al. (1987) suggested ionised filamentary structures or plasma sheets as possible plasma lenses.If a cylindrical filament is viewed end-on, then a high column density is observed with a small transverse size.Thus, an illusion is created where if the observer assumes the object is spherical (as Clegg et al. (1998) did), then the density and inferred pressure would be greatly overestimated.Similarly, plasma sheets viewed edge-on or many plasma sheet layers could provide the high observed column densities.Such plasma sheets could exist in the ionization fronts or cooling instabilities of old supernova remnants where small scattering angles could occur through hundreds of face-on sheets, on average (Romani et al. 1987).The most extreme scattering, however, would occur due to corrugations or inhomogeneities in the sheets' structure (Romani et al. 1987;Pen & Levin 2014).Similarly, ionised shells of molecular clouds could provide the necessary scattering conditions (Boldyrev & Königl 2006).In the Boldyrev & Königl (2006) model, the scattering would be dominated by a few strong scattering events as opposed to the culmination of multiple smaller scattering events.
Unaware of extreme scattering by plasma lenses, Goldreich & Sridhar (2006) proposed sheet-like structures primarily aligned along the line-of-sight were responsible for extreme scattering in the Galactic Centre.Motivated by Goldreich & Sridhar (2006), Pen & King (2012) model plasma lenses as under-dense current sheets ". . .as a generic solution to strong interstellar scattering" including ESEs.Interestingly, these under-dense sheets act as convergent lenses as opposed to divergent.A real-time survey with the Australia Telescope Compact Array surveying ∼ 1000 AGNs over 4−8 GHz identified an ESE in the direction of PKS 1939-315 in 2014(Bannister et al. 2016).A measurement of the electron column density as a function of time N e (t) for this event was achieved.Since the peak of the N e (t) profile corresponded with the minimum of the flux density profile, it was concluded that the ESE was caused by a divergent lens.Such evidence contradicts the Pen & King (2012) model.
ESEs are distinguished by their dramatic decrease in flux density with time.Lazio et al. (2000) conducted multi-epoch Very Long Baseline Interferometry (VLBI) observations of PKS 1741-038 during an ESE.The source appeared to increase in angular size during the event.By the conservation of surface brightness of a purely refractive lens model, an increase in angular size should correspond to an increase in brightness, not the decrease that is characteristic of ESEs.Lazio et al. (2000) explained their observations using both turbulent and refractive effects.Pen & King (2012) suggested that unresolved multiple images could explain the apparent increase in angular size.
Other "extreme scattering" phenomena such as intraday variability (IDV) of radio quasars, parabolic arcs in pulsar secondary spectra (Stinebring et al. 2001), and Galactic Centre scatter-ing (Goldreich & Sridhar 2006) have been observed.All such phenomena are suspected to arise from the scattering of radiation by turbulent media, or plasma structures.All these phenomena are also plagued by the problem that the scattering plasma structures must be overpressured with respect to the diffuse ISM, to various degrees.Although it is suspected that turbulent media or plasma structures are responsible for the extreme scattering, it is unknown if they are indeed connected phenomena or not.
Only a handful of ESEs have been discovered even with extensive surveys (Fiedler et al. 1994;Lazio et al. 2001;Pushkarev et al. 2013).Thus, ESEs suffer from small number statistics and their detection rate is inconclusive.This is further troubled by difficulties distinguishing between "small" and "high" amplitude variability.Short time-scale, small-amplitude ESEs are easily confused with interstellar scintillation, while long-timescale small-amplitude ESEs can be contaminated by intrinsic source variability (Stanimirović & Zweibel 2018).Automated methods such as the wavelet transformation of light curves by Lazio et al. (2001) have been able to identify ESEs with moderate accuracy.While Lazio et al. (2001) were able to identify 15 events from 12 sources, five ESEs previously found by eye were not identified using their wavelet transformation tool.Most ESE discoveries have been due to chance occurrence.Bannister et al. (2016), however, were able to achieve a real-time ESE observation by monitoring ∼ 1000 active galactic nuclei for changes in the continuum corresponding to the wavelength-squared dependence of the plasma refractive index indicating the beginning of a possible ESE.Dedicated monitoring for ESEs are needed to gather the necessary statistics and measurements that characterize ESEs.
Since Clegg et al.'s (1998) model, more complicated models have been developed to explain ESEs and the plasma lenses that are theorized to produce them.Whereas Clegg et al. (1998) considered a lens with an electron column density distributed as a two-dimensional Gaussian on the sky, Er & Rogers (2018) considered power law, and exponentially distributed electron volume densities.Er & Rogers (2018) found interesting behaviours in the magnification curves of these models including similar magnification curves between the well-studied Gaussian lens and their softened power-law lens.Rogers & Er (2019) generalized plasma lensing models to dual-component lenses using families of exponential and softened-power laws, finding more complicated magnification curves and caustic structures.Interestingly, Rogers & Er (2019) pointed out a degeneracy between two specific models -a dual-component Gaussian lens, and a single spherically symmetric exponential lens.Since the index of refraction of plasma varies with the wavelength of the signal, this degeneracy may be broken by the active observation at various wavelengths of an ESE (Rogers & Er 2019).Er & Rogers (2019) further expanded the library of analytical plasma lens models to elliptical lenses following exponential and softened power law distributions.
Extending beyond spherical and elliptical lenses, Rogers et al. (2020) explored toy filamentary models confined by helical magnetic fields.The cylindrical geometry, when viewed along the long axis, may create the illusion of a small transverse size without the issue of being over-pressured.Rogers et al. (2020) were able to reproduce the magnification of spherical lenses even with the cylindrical geometry, but determined that these models were physically unrealistic from the resulting dispersion measure.
More recently, Er et al. (2022) investigated plasma lensing as modeled by a double-plane, each with a projected Gaussian-distributed electron density.Such an approach is akin to scattering by a pair of effective planes as opposed to volumetric scattering.Effective scattering planes are often used to describe the effects of gravitational lensing (Schneider et al. 1992), and scintillation (Hill et al. 2005;Reardon et al. 2020;Sprenger et al. 2022;Stinebring et al. 2001, for example).While some image properties from a single Gaussian lens mimicked those produced by a double-lens, there were some properties which distinguished between the two scenarios.In particular, the time delay (and time delay-frequency relation) between the two models were sufficiently different.

Problems with Plasma Lens Models
There are problems concerning the properties of the plasma structures that are predicted to produce ESEs, which speak more broadly to the physics of the ISM.The primary and most accessible problem to be addressed is to better constrain the geometries and physical characteristics of the plasma lenses that produce ESEs.Consequently, additional questions may be addressed by relating the physical properties of plasma lenses to the problem of over-dense, over-pressured models.Given the relatively small scale of these plasma structures and their predicted high densities and pressures, it is unknown how prevalent these objects are, what the origins of these plasma structures are, and how they might influence the physics of the ISM.ESEs provide an opportunity ot study dense structures in the ISM at small length scales that woul be difficult to access otherwise.A plausible model of ESEs would provide new insights in the the physics of the ISM.
One primary concern with plasma lenses is how they are predicted to be highly overpressured relative to the surrounding ISM.Without some confinement mechanism, their high pressures suggest that plasma lenses must be transient structures.Models based on thin sheets viewed (nearly) edge-on could provide a means to reduce the lens pressure.It is possible that other geometries and methods of confinement could provide the necessary conditions to produce the lenses required for ESEs (Romani et al. 1987;Rogers et al. 2020, for example).Between theory and observations, nothing is conclusive at this point concerning the geometric arguments made to address the over-pressure problem.
A big mystery concerns the origin of ESEs and the plasma lensing structures that produce them.Whether these ESEs originate from ionised turbulence or from discrete lensing structures is unknown.Supernova remnant (SNR) shells, bubble walls, cloud edges, and shocks have been postulated as sites of ESEs (Romani et al. 1987).Hamidouche & Lestrade (2007) support the idea that ESEs can be produced by turbulence, however the few instances of multiple-imaging detected during an ESE suggest they are produced by discrete plasma lenses.A turbulent origin of ESEs has also been suggested, but volumetric turbulence, or volumetric ray-tracing more generally, has yet to be explored in investigating ESEs (Stanimirović & Zweibel 2018).
Investigating ESEs and plasma lenses will provide significant insights into the physical processes of the ISM.These small structures are on the scale of turbulent dissipation, and their properties may yield clues informing poorly understood physical processes at these scales.These overpressured plasma structures could be a source of significant heating in the ISM depending on their volume filling factor.A measurement of the volume filling factor, however, has not yet been conducted.Stanimirović & Zweibel (2018) suggest tiny scale ionised structures (TSIS) such as plasma lenses are a ". . .unique probe of stellar feedback efficiency since they are likely produced by local turbulence enhancements in places like SNRs, stellar bubbles, and winds" and ". . .deserve a major spotlight and attention."

Method Outline
Our aim is to construct ESE light curves using numerical volumetric ray-tracing methods.Rays are traced between the source (such as a pulsar) and the observer through a lensing medium.The geometry of the system is set up similarly to a gravitational lensing scenario such as Narayan & Bartelmann (1996, Figure 5).The source is placed a distance D s away from the observer, the lens' origin is placed a distance D d from the observer, and the distance between the lens and the source is given by D ds = D s − D d .In fact, the problem itself draws many analogies with gravitational lensing, which we use as a test of our method (Section 7).
Ray-tracing is calculated by solving the eikonal equation through a medium of spatially varying refractive index (Section 2).The eikonal method requires that a model specifying n(⃗ x) is provided.Multiple-imaging is accounted for by mapping pixels on the observer plane to overlapping or clustered points on the source plane.The mapping of these pixels is done by propagating a wave front formed by Delaunay triangulation over rays traced at the triangle vertices (Section 3).Consequently, magnification maps (Section 4), intensity maps (Section 5), and light curves (Section 6) are constructed from mapping the brightness of a source to calculate the intensity by the observer/detector.Whereas much of the literature has modeled plasma lenses and the corresponding light curves analytically, numerical ray-tracing methods offer an opportunity to explore new lens models, especially those defined by volumetrically extended structures such as regions of plasma turbulence.By approaching the investigation numerically, solutions can be made much more general.Instead of projecting scattering effects onto an effective two-dimensional scattering plane, we use full three-dimensional lens models.

RAY-TRACING
Snell's Law adequately describes refraction at the interface between two uniform media.More realistically, when the medium that light passes through varies continuously in refractive index, Snell's Law becomes inadequate to describe the trajectory of light through the medium.In regions where the wavelength of light is small compared to the size of the optical system, the amplitude and direction of propagation of a wave does not change considerably for several wavelengths and can be approximated as a plane wave.General solutions of the wave equation are described with a frequency ω and wave vector k = ω/c (Römer 2006) (3) where φ k (⃗ x) is a real number quantifying the wave amplitude, and S(⃗ x) is a real function known as the eikonal function.The eikonal function determines the phase of the wave at position ⃗ x such that S(⃗ x) = constant represents surfaces of constant phase.Ray-tracing is performed by solving the eikonal equation where the index of refraction n(⃗ x) spatially varies, ⃗ x are the ray trajectories, and the arc length s serves as the curve parameter.The eikonal equation describes families of rays with trajectories orthogonal to the level surfaces of constant phase and is solved numerically by re-parameterizing Equation 4 in terms of time t.To begin, we define a parameter where ⃗ v is the ray velocity and v is the corresponding speed in the medium is given by With ⃗ q and ds/dt, Equation 4 is rewritten as Using the definition of ⃗ q (Equation 5) and ds/dt (Equation 6a), the velocity of the ray trajectories are d⃗ Equations 6a to 6c are solved numerically for the ray trajectories using the initial conditions for the position ⃗ x 0 and ⃗ q 0 = n⃗ v 0 /∥⃗ v∥ (Equation 5).

Numerically Solving for the Ray Trajectories
Computing ray trajectories requires integrating the parameterized eikonal equation (Equation 6).A custom fifth-order Runge-Kutta (RK5) integrator based on the Dormand-Prince method (Dormand & Prince 1980;Press et al. 2007) is written and tested to numerically solve the parameterized system.A custom RK5 integrator is primarily advantageous for tracing many rays simultaneously and managing stopping conditions.
Rays are initially placed on the image or observer plane (the terms "image" and "observer" planes will be used interchangeably) such that each ray forms the vertices of a pattern of equilateral triangles (Figure 1).From a distant source such as a pulsar or quasar, the rays arrives at the observer (Earth) with a small angle.If gravitational lensing and Einstein rings are indicative of the angle of light rays detected, it suffices to launch the initial rays with an angle of order ∼ θ E (the Einstein angle) primarily in the normal (+z) direction from the observer (x-y) plane toward the source plane.For our volumetric ray-tracing method with plasma lenses, θ E serves as a starting point to launch the initial rays, which is then fine-tuned.Although perhaps it is counter-intuitive to trace rays from the observer plane to the source plane, this is necessary to account for multiple-imaging (see Section 3).Once the rays are launched from the initial positions, they are propogated through a medium with spatially varying continuous index of refraction by solving the eikonal equation with the RK5 integrator.Several models of plasma structures with spatially-varying refractive index are discussed in Section 8.The integrator is halted once the rays reach the source plane.Parameters such as ⃗ x, ⃗ q, and s on the source plane are obtained by interpolation of the ray one time step before and after the ray crosses the source plane.
A map is created from ray-tracing since both ray positions on the observer plane and the corresponding position on the source plane are known.This mapping is the cornerstone of this investigative technique for exploring plasma lenses and ESEs because it links the brightness distribution between the source and the observer regardless of how the medium is structured.It allows for the mapping of pixels of a detector (the observer plane) to the brightness from a distance source (Section 5) and the construction of light curves (Section 6).

WAVEFRONT PROPAGATION WITH DELAUNAY TRIANGULATION
Ray trajectories are directed in the normal direction with respect to the wavefronts.By calculating the ray trajectories through a refracting medium, wavefront elements are approximated by Delaunay triangles.For a two-dimensional set of points, the Delaunay triangulation criterion requires that the circumscribed circle or circumcircle of each triangle contains no other points in its interior.Consequently, the formation of Delaunay triangles is such that triangles with larger minimum angles are formed with preference over those with smaller minimum angles.It is worth noting that the generation of Delaunay triangles can be degenerate for a given set of points, but this is inconsequential as long as the scheme covers the triangulated area.MATLAB's Delaunay functions create a Delaunay object that contains the points or vertices of each Delaunay triangle, and also contains a "connectivity list" establishing the connection between the three vertices forming a Delaunay triangle.Points are listed in an array such that each row is a point whose x, y, z coordinates are listed by column.The connectivity list contains the row numbers of the three vertices of each triangle, which point to the coordinates of each vertex, and is established at the initial positions of the ray trajectories (Figure 1).Once generated on the image plane, the connectivity list provides a static list of connections between triangle vertices defined by rays.Thus, this connectivity list is preserved throughout the rays' trajectories, effectively tracking the propagation of the wavefront.

Accounting for Multiple Imaging
Delaunay triangles provide a significant advantage for handling multiple imaging.Recall that this technique begins with rays on the observer plane which are ray-traced through a medium with spatially varying refractive index to intersect the source plane.Delaunay triangles on the source plane may be distorted from their initial positions on an equilateral grid on he observer plane.Each Delaunay triangle is identified and tracked thanks to the connectivity list.
Accounting for multiple imaging involves creating a mapping between pixels on the observer plane and the positions where these pixels end up as points on the source plane after travelling through the medium.Pixel centres are first created on the observer plane in a grid formation.Each pixel is then identified with the Delaunay triangle it resides within.Assuming barycentric coordinates are preserved as the wavefront of Delaunay triangles propagates through the medium, then the pixels are mapped onto the source plane since the connectivity list of Delaunay triangles between the observer plane and source plane are known (Section 2).Some of the points may be closely clustered together or even overlap on the source plane.Consequently, this mapping accounts for multiple-imaging as the brightness of many pixels on the observer plane map to the same point or cluster of points in close proximity on the source plane (Figure 2).This mapping is used to generate magnification maps (Section 4) and intensity maps (Section 5).

MAGNIFICATION MAPS
Delaunay triangles provide an intuitive way to understand and calculate magnification.Recall that the Delaunay triangles are formed from connecting the vertices or rays that propagate through a medium with spatially varying index of refraction.The rays may scatter and take non-uniform, irregular trajectories through the medium, which deforms the wavefront or surface of Delaunay triangles.The magnification in a triangle is given simply by where A i and A s are the respective areas covered by a mapped triangle on the image plane and on the source plane.The magnification is assigned to the centroid of each triangle on the observer plane and interpolated onto the observer plane pixels to generate a magnification map.As no further calculations are performed from the magnification map, one magnification map is shown in Section 8.1 for brevity.Intensity maps provide a useful way to visualize the images created by this method of ray-tracing and wavefront propagation, and are a precursor to constructing light curves (Section 6).Mapping pixels on the observer plane to points on the source plane (Section 3.1) is key to constructing intensity maps.

INTENSITY MAPS
First the brightness distribution of the source is constructed on the source plane.The brightness profile of the source is typically assumed to be Gaussian distributed where σ 2 is the variance (σ is the standard deviation) and ⃗ x s is the centre of the source, although any brightness distribution may be used.An intensity map on the source plane is produced by first creating a grid of pixel centres on the source plane and then calculating the brightness at each pixel centre based on the source's position.
Having established the mapping of points on the source plane to pixels on the observer plane, the brightness of the points on the source plane is mapped directly onto the pixels of the observer plane.Note that this mapping does not directly use the source plane pixels; rather, the image plane pixels are mapped back to coordinates on the source plane, allowing us to assign intensity values from the source model to the image plane pixels.This procedure generates an intensity map of the image, which automatically accounts for multiple imaging since multiple pixels on the observer plane may sample the same points on the source plane.Instances of multiple imaging are clear especially in the examples and test cases such as gravitational lensing by a point mass presented in Section 7.

Subpixels
Creating intensity maps from mapping pixels to points on the source plane is an effective way to visualize images resulting from lensing systems.However, to decrease noise in the intensity map, each pixel of the observer plane is subdivided into subpixels to calculate an averaged intensity value.In this work, the subdivision of pixels is controlled by an oversampling parameter (OSP ) such that pixels are divided into OSP × OSP subpixels.Typically OSP = 2 is chosen as a good compromise between image fidelity and computation speed.The same mapping of pixels from Section 3.1 is used to map the subpixels' barycentric coordinates from the observer plane to the source plane, the source brightness is calculated for each point (mapped from a subpixel), and the brightness is mapped back to the observer plane to construct the intensity map.The difference is, however, that the intensity of a pixel is calculated as the average over its subpixels.With mapping the subpixels to the source plane, the region where the points are on the source plane is more thoroughly sampled.Consequently, the intensity averaged from the subpixels is a more accurate reflection of the brightness of the region sampled compared to using single pixel centres.
6. LIGHT CURVES Time-dependent changes in the lensing system may cause changes in the total intensity measured, resulting in a light curve.We assume that the lensing medium and observer are stationary while the source object moves, and we construct a light curve from intensity maps generated for a sequence of time steps.Telescope beam effects can be (optionally) introduced in this step by assuming a beam comprised of spatially variable sensitivity on the observer plane.We typically assume a Gaussian beam sensitivity profile or a flat beam sensitivity, which is the default.
To calculate the total intensity from the intensity maps (Section 5), the intensity of each pixel in the map is multiplied by the area of the pixel, and summed for all pixels in the map.This process is repeated for a sequence of times steps to generate a light curve.The light curve is normalized by the total intensity on the image plane when the source object is far off-axis with the lens and the lensing effects are negligible.

TEST CASES
We test our ray-tracing method by comparing to a simple point mass gravitational lens solution, and diverging point lens solution, both with known analytic solutions.
For numerical calculations, it is conventional to compute in dimensionless or scaled quantities, which we denote with an overscript tilde.

Gravitational Point Mass Lens
Gravitational lensing by a point mass, and gravitational lensing more generally, is well-studied (Schramm & Kayser 1987;Kayser & Schramm 1988;Schneider et al. 1992;Narayan 1992;Narayan & Bartelmann 1996, for example) and there exists an especially simple analytical solution for lensing by a point mass.As light propagates from a source, its trajectory follows geodesic paths in the curvature of spacetime, which are observed as curved paths.For small perturbations by a gravitational potential Φ in locally flat Minkowski spacetime, the curvature of light by a gravitational lens is described by an effective index of refraction (Schneider et al. 1992;Narayan & Bartelmann 1996) For a gravitational lens that is a point mass M , the gravitational potential is given by Consequently, the effective index of refraction for a point mass is which can be simplified by collecting the constant terms such that where k = 2GM/c 2 is the corresponding Schwarzschild radius.The Schwarzchild radius provides a natural length scale such that the effective index of refraction (Equation 12) in terms of the dimensionless radius r = r/r 0 , where r 0 is the length scale (r 0 = k in this case), is written as The corresponding dimensionless gradient is With n and ⃗ ∇n modeled for the point mass lens, the eikonal equation (Equation 4) is solved numerically for the ray trajectories (Section 2), followed by the Delaunay triangle wavefront (Section 3), which results in the magnification map (Section 4), intensity map (Section 5), and light curve (Section 6).16) is overlaid on top of the numerical result.The two results match on the order of ∼ 10 −6 , which supports the validity of the numerical method.Parameters were chosen such that k = 1 AU and D s = 1 kpc, which are sizes representative of pulsar ESEs.
Einstein rings are a visually striking consequence of multiple-imaging due to gravitational lensing by a point mass M .When a source lies directly behind the lens along the line of sight, it is imaged as an Einstein ring of characteristic Einstein angle which is a well-known analytic and observed result (Schneider et al. 1992;Narayan & Bartelmann 1996, for example).Scaling θ E by the Schwarzchild radius yields To test the validity of the ray-tracing procedure described previously, the resultant size of the Einstein ring on the observer plane intensity map must agree with the analytic Einstein ring size.Figure 3 demonstrates that the numerically ray-traced pixels on the source plane map to an Einstein ring of the same size as the analytic Einstein ring as calculated by Equation 16, which supports the validity of the numerical ray-tracing method.
Furthermore, the light curves computed numerically by ray-tracing should agree with the analytic result for a point mass gravitational lens (Narayan & Bartelmann 1996) where y = β/θ E and β is the angle on the source plane with respect to the optic axis of the lensing system (Narayan & Bartelmann 1996, Figure 5).Figure 4 shows the normalized light curves for a point mass lens.Red lines of different shades in Figure 4 represent different sizes of the source brightness, which are compared with the analytic solution for a point source (black dashed line).
As the Gaussian source becomes more point-like, the solution converges to the analytic solution, supporting the validity of this method.8) becomes smaller, shades of red lines) to the analytic solution (black dashed line) for a point mass gravitational lens.Total intensity or magnification µ is normalized, and is plotted as a function of angle on the source plane β.The gravitational point mass and source are aligned at β = 0.A slight bulging effect is noticeable in the σ = 0.25 light curve due to the small source size moving moving across regions with a higher density of points on the source plane.

Point Diverging Lens
As a toy model, the point diverging lens serves as a useful test case to compare the numerical ray-tracing with another analytic result (Safonova et al. 2001).In general, the index of refraction for the point diverging lens is described by which is written in dimensionless form as where 0 ≤ k ≤ 1. Comparing the point diverging lens (Equation 18) to the point mass gravitational lens (Equation 12), it is easy to notice that the two equations are different only by a change in sign of the second term resulting in a diverging lens.Similar to the gravitational point mass lens (Section 7.1), the light curves produced numerically by ray-tracing should converge to the analytic light curve for a point source.Safonova et al. (2001) found the toy divergent point lens model produced a light curve described by Figure 5 demonstrates the convergence of the numeric solution as the size of the source becomes point-like (red lines) to the analytic solution (black dashed line).The convergence of the numeric light curve to the analytic light curve further supports the validity of the numerical ray-tracing method.8. PLASMA LENS EXAMPLES Having tested our volumetric numerical ray-tracing method against the analytic cases of Section 7, we now use volumetric plasma lens models.Specifially, we will use a Clegg et al. (1998)-inspired Gaussian lens and an analytic model of a non-self-gravitating, magnetized, rotating filament (Grafton et al. 2023).
In the most general case, the model for a plasma lens describes the particle mass density ρ.Assuming most of particles are hydrogen ρ ≈ ρ H , the hydrogen mass density ρ H is converted to an electron number density via where X is the ionisation fraction, µ is the mean molecular mass, and m H is the mass of a hydrogen atom.For electromagnetic radiation at frequency f propagating through a cold, ionised plasma such as the ISM, the electron number density contributes to the calculation of the plasma frequency f p (Equation 2), which in turn determines the index of refraction n(⃗ x) (Equation 1).The refractive index and its gradient ⃗ ∇n are the main components needed to solve eikonal equation (Equation 6) and calculate ray trajectories (Section 2).The signal frequency f = 10 9 Hz is adopted as ESEs are detected in the radio (∼ GHz).

Gaussian Lens
where the scaling of the spherical radius follows r = r 0 r = σr, σ is the standard deviation, and n e0 = ρ 0 X/µm H is the electron volume density at the core r = 0.The electron volume density is used to calculate the plasma frequency f p (Equation 2), which is used to calculate the refractive index n (Equation 1).Since the Gaussian lens possesses spherical symmetry, the dimensionless gradient of the refractive index is Recall that the overscript tilde notation denotes dimensionless or scaled quantities.
Unlike the gravitational point mass lens (Section 7.1) and point diverging lens (Section 7.2), the Gaussian lens is a plasma lens, so we account for the effects of electromagnetic waves propagating through plasma via f p (Equation 2).Consider a system with a Gaussian lens placed halfway between an observer on Earth and a pulsar ESE event.The Gaussian lens has a length scale (standard deviation) r 0 = σ = 10 AU, and hydrogen mass density ρ H = 10 −20 g cm −3 , typical of plasma lenses based on observations and models of ESEs (Clegg et al. 1998;Stanimirović & Zweibel 2018).Furthermore, pulsars in the Galaxy are located on the order D s = 1 kpc from the Earth, and the lens is placed halfway in between D d = 0.5 kpc.The resultant source and observer planes, and light curve are shown in Figures 6 and 7, respectively, and a magnification map is shown in Figure 8.A snapshot of the source is shown in the left panel of Figure 6 with the corresponding observer image on the right panel of the same figure.White dots in the source plane (which are too small and numerous to be individually resolved in Figure 6, left panel) illustrate the correspondence of the pixels on the observer plane mapped to the source plane after ray-tracing and Delaunay triangulation (Section 3). Figure 7 shows the complete normalized light curve after the source has travelled horizontally from the left edge of the source plane to the right edge.

Rotating Magnetized Filament as a Plasma Lens
Grafton et al. ( 2023) constructed analytic models of rotating magnetized filaments observed in the extreme region of the Galaxy's Central Molecular Zone (CMZ) known as molecular tornadoes.The formation of molecular tornadoes is theorized to arise from the shearing of a magnetic flux tube, which drives a torsional Alfvén wave that propagates along the tube and triggers a cylindrical instability.Both axial (poloidal, B z ) and toroidal (B ϕ ) magnetic fields thread the filament, which are characterized by the poloidal and toroidal flux-to-mass ratios formulated by Fiege & Pudritz (2000) Γ  and taken as positive quantities.From the ideal magnetohydrodynamic (MHD) equations, Grafton et al. (2023) found the radial mass density profile for the non-self-gravitating filament in equilibrium (Grafton et al. 2023, Equation 10) where P ∞ is the total pressure including both gas and magnetic pressures everywhere within the filament and out to radial infinity, and σ is the velocity dispersion.Here, we re-imagine the Grafton et al. ( 2023) molecular tornado as a smaller scale rotating plasma filament.Specifically, the case of constant Γ ϕ and Γ z (Grafton et al. 2023, Section 3) is the filamentary model used here as a plasma lens.The same length scale as Grafton et al. ( 2023) is adopted where with slight changes to the scaling relations between dimensionless quantities and physical quantities The scaling relations (Equation 27) have free parameters r 0 , ρ 0 , and σ.ESE plasma lenses are on the order of r 0 ∼ AU, while typical molecular tornadoes have density ρ 0 = 10 −21 g cm −3 , and velocity dispersion σ = 15 km s −1 (Sofue 2007;Rathborne et al. 2014;Au & Fiege 2017).Au & Fiege (2017) previously explored a parameter space 10 −3 ≤ Γz ≤ 100 for their molecular tornado models.As a rough estimate, assuming B ≈ B z = 0.4 mG (Sofue 2007), ρ 0 = 10 −21 g cm −3 , and σ = 15 km s −1 , Equations 24a with 27 reveal Γz is of order unity.
In dimensionless form, the radial mass density profile of the filament with constant Γ z and Γ ϕ is (Grafton et al. 2023, Equation 21) which may be most conveniently scaled to a physical quantity via Equation 27 to calculate the refractive index (Section 8).The dimensionless gradient of the refractive index is where (Grafton et al. 2023, Equation 24) For example, consider a similar system to the Gaussian lens previously shown in Section 8.1, except now with the filament lens such that the length scale is r 0 = 10 AU, D s = 1 kpc, and D d = 0.5 kpc.The filament is oriented such that its long axis runs vertically from top-to-bottom, perpendicular to the observer's perspective.In this orientation, the filament's inclination angle is defined as i = 0. Rotating the top of the filament down such that the filament axis gradually becomes more aligned with the observer's line-of-sight corresponds to an increase in i.The resulting light curves are shown in Figure 9 for Γz = 10 −3 .

DISCUSSION
To explain the results shown in Section 8, consider Figures 6 and 7 for the Gaussian lens (Section 8.1).Some notable features include multiple-imaging seen in the observer plane in the right panel of Figure 6.The apparent radial stretching of the source on the observer plane due to multipleimaging is similar to other divergent lens models (Safonova et al. 2001;Er & Rogers 2018;Rogers & Er 2019).Notable in Figure 7 are the caustic spikes around x/r 0 = 5 and dramatic drop in flux around x/r 0 = 0 which are signature features of ESEs.Comparing Figure 6 with 7, it is easy to see that many pixels on the observer plane are mapped the same areas on the source plane as evidenced by the dense circular region of white point on the source plane centred on (0,0).When the source crosses these regions of dense points on the source plane, its brightness is mapped to multiple pixels on the observer plane, which results in multiple-imaging and an increase in brightness in the light curve.The opposite is true of the drop in the light curve in Figure 7; the source traverses through a region devoid of points on the source plane, which results in little to no source brightness mapped to the observer plane.This point is discussed further in Section 9.1.2023) filament reimagined as as plasma lens for inclination angles i = [0, π/6, π/4, π/3] (shades of red).The source travels from left to right in the source plane, perpendicular to the long-axis of the filament when i = 0, and takes the same path for each inclination angle.The lensing system is set up with the filament half way between the observer and source plane that D s = 1 kpc, D d = 0.5 kpc, with length scale r 0 = 10 AU, core density ρ 0 = 10 −21 g cm −3 , and Γz = 10 −3 .
A lack of caustic spikes is notable in the light curves of Figure 9.The light curves are mainly flat far away from the filament core and gradually decreasing more dramatically toward the filament core.Interestingly, these light curves resemble some of those produced by the filamentary plasma lens models in Rogers et al. (2020).As the inclination angle increases and the filament axis becomes more aligned with the observer's line-of-sight, the drop in the light curve decreases more.Although these light curves do not possess the characteristic caustic spikes of ESE light curves, the relationship between an increasing inclination angle with a larger drop in the light curve is suggestive that the excessive density and pressure of plasma lenses are a geometrical effect (Romani et al. 1987).When the filament's long axis is more aligned with the observer, the greater deflection in the light rays away from the observer occurs as the rays travel through a region of greater electron column density.
It is also worth noting that the core densities used in the plasma lens models (Section 8) are still a few orders of magnitude in excess of typical ISM values.A wider exploration of the parameter space is underway to investigate possible models without excessive densities and produce ESEs (Section 9.3).

Practical Considerations
As with many computational problems, computational speed is a limitation to the accuracy of the solution and practicality of solving the problem.Although physical hardware limitations do not necessarily affect the overall results and conclusions, it is important to recognize and to be cognizant of the limitations and consequences.For example, the memory limit is particularly important in choosing the number of rays and number of pixels.Generally, the calculation of identifying pix-els on the observer plane with the Delaunay triangles they reside within (Section 3.1) is the most memory-intensive procedure.While this calculation is vectorized to increase the calculation speed, vectorization also increases memory use.There are a few consequences to limiting the number of rays and number of pixels.Qualitatively, the number of rays is most easily observed when the pixels on the observer plane are mapped to and shown on the source plane.When relatively fewer rays are traced, the points on the source plane are not as smoothly distributed.This is because fewer Delaunay triangles are created, resulting in larger Delaunay triangles containing more pixels per triangle, which show up as points distributed in triangular-shaped artifacts on the source plane, and some regions may be noticeably more densely populated with points than others.The resulting light curves may suddenly jump or drop in intensity.
Reducing the number of pixels reduces the number of points sampled by the source on the source plane, meaning less of its brightness mapped to the observer plane and accounted for in the calculation of the total intensity.In severe cases, there may be regions on the source plane which are devoid of points.When a significant portion of the source's total brightness (typically at or around the peak of the Gaussian brightness profile) is located in the gaps between points, the light curve shows a significant drop in intensity.If the source encounters several of these gaps then the light curve will appear "bumpy" or "wavey".The slight ringing effect seen in some curves of Figure 5 demonstrate this effect to a minor degree.This issue is largely mediated by increasing the number of pixels such that the source region is well populated with points, and to a lesser extent by the number of Delaunay triangles.For a fixed number of pixels residing within a Delaunay triangle on the observer plane, a larger triangle on the source plane can spread out the points more compared to a smaller triangle on the source plane.Additionally, increasing the size of the source's Gaussian brightness profile helps mediate this problem, and there are other consequences to consider concerning the source size, which are discussed next.

Source Size
Ideally, the source of the radio signal would be a point source to replicate a pulsar ESE scenario.However, as previously mentioned, practical limits to the number of pixels translates to the distribution and population of points on the source plane.If the source is smaller than the gaps between the points then the same problem arises, as previously discussed, where a "bumpy" light curve results.
In some situations, it is also difficult to tell whether some features in a light curve may be a natural consequence of the lensing model or artifacts due to a combination of the source size and finite number of points on the source plane (or pixels on the observer plane).For example, while generating light curves for the Gaussian plasma lens, some sets of parameters produced what seemed to be outer caustics like those shown in Clegg et al. (1998).However, many of these light curves are riddled with artifacts due to the comparably small Gaussian brightness profile size of the source relative to the size of gaps between source points.It is possible that the overall shape of the light curve including the outer caustics are genuinely from the lens model, but the inability to confidently sample a small source size casts some doubt into the interpretation.The lack of outer caustics in our volumetric Gaussian lens light curve (Figure 7) could also be due to a genuine difference between Clegg et al.'s (1998) electron density model distributed as a two-dimensional Gaussian on the sky.Also of note is that features in the light curve such as caustics may be missed or smoothed out depending on the size of the source relative to the density of points on the source plane.Clegg et al. (1998) also remark and demonstrate for their model that the light curve does change depending on the source size, however.
Despite the fact that the source can never be a perfect point, tests converge to known analytic solutions in the limit where the source becomes smaller and more point-like as previously shown in Section 7.

Future Work
The numerical method presented here lends itself well to ray-tracing through simulated volumetric turbulence, which has not yet been explored.Plasma lenses producing ESEs have so far been modeled as discrete structures.As mentioned previously (Section 1), it is theorized that ESEs due to TSIS may be the result of density fluctuations from turbulence.Often, the (possibly related) phenomenon of interstellar scattering and scintillation arcs are statistically modeled and/or by an effective scattering screen (Cordes et al. 1985(Cordes et al. , 2016;;Pen & King 2012;Simard & Pen 2018;Stinebring et al. 2001Stinebring et al. , 2022;;Sprenger et al. 2022;Baker et al. 2022;Jow & Pen 2022, for example).The eikonal ray-tracing method explored here opens the possibility for investigating three-dimensional volumetric turbulence as a lens.Efforts are currently underway to apply the numerical ray-tracing method presented here to density fluctuations following Kolmogorov turbulence, and simulated turbulence from PythonMHD, a Python-focused astrophysical MHD code under development at the University of Manitoba (Leboe-McGowan 2022).Sheet-like structures are also being considered as their geometry may explain the overpressure problem when viewed edge-on (Section 1).A wider exploration of the parameter space for the Gaussian lens and filament lens inspired by Grafton et al. (2023) is underway to understand the overpressure problem which has plagued plasma lens models so far.

CONCLUSION
Whereas prior studies have often utilised effective two-dimensional planes to study lensing, we have developed a numerical method to explore volumetric plasma lensing configurations that produce ESEs by solving the eikonal equation to trace ray trajectories.Rays start from a grid of pixel centres on the observer plane and are traced to the source plane to account for multiple imaging.Delaunay triangles connect adjacent rays to approximate the wavefront.The pixels are divided into subpixels, and assuming the barycentric coordinates are preserved throughout the trajectory, form a map between the observer and source planes.Gridded subpixels on the observer plane are mapped as points on the source plane.An intensity map is generated by mapping a Gaussian-distributed brightness profile on the source plane back to pixels on the observer plane.A light curve is generated by summing the total intensity contributed by all the pixels on the intensity maps generated as the source brightness profile moves across the source plane.Tests with known analytic solutions demonstrate convergence with our numerical solutions as the Gaussian source brightness distribution becomes more point-like.Examples of lensing models inspired by Clegg et al. (1998) and Grafton et al. (2023) are shown with characteristic ESE light curve features in the former.Further investigations include exploring a larger parameter space for the examples presented here, sheet-like lenses, and volumetric turbulence.support of the University of Manitoba Graduate Fellowship (UMGF).J.D.F.acknowledges the support of a Discovery Grant from NSERC, RGPIN/5930-2017.

Figure 1 .
Figure 1.Illustrative figure of initial ray positions (large red dots) on the observer plane distributed in a grid of equilateral triangles (Section 2.1).Delaunay triangulation (blue lines) are formed from connecting the initial ray positions, which act as the triangle vertices (Section 3).A grid of pixels (small black dots) are identified with the Delaunay triangle they reside within (Section 3.1).

Figure 2 .
Figure 2. Illustrative figure of lensed ray positions (large red dots) on the source plane (Section 2.1) corresponding to Figure 1.The connectivity list tracks the Delaunay triangle wavefront (blue lines) from the observer plane to the source plane (Section 3).Pixels on the observer plane are mapped onto the source plance as points (small black dots) because their barycentric coordinates are preserved (Section 3.1).A cluster of points, Delaunay triangles, and rays in the centre of the image are too densely spaced to be individually distinguished.Multiple-imaging is conducive in the centre of the figure where the points are most tightly clustered because the brightness within that small region will map to many pixels on the observer plane.This specific figure uses a point mass gravitational lens (Section 7.1) for illustrative purposes only.

Figure 3 .
Figure 3. Testing the validity of the ray-tracing method by comparing the numerical results with analytic results of an Einstein ring.r is the angular radius of the Einstein ring and θ is the angle around it.The numerically ray-traced pixels from the source plane are mapped to an Einstein ring (multi-coloured) on the observer plane.A red line representing the analytic Einstein ring (Equation16) is overlaid on top of the numerical result.The two results match on the order of ∼ 10 −6 , which supports the validity of the numerical method.Parameters were chosen such that k = 1 AU and D s = 1 kpc, which are sizes representative of pulsar ESEs.

Figure 4 .
Figure 4. Convergence of the numerical ray-traced solution as the Gaussian source brightness becomes more point-like (the standard deviation σ (Equation8) becomes smaller, shades of red lines) to the analytic solution (black dashed line) for a point mass gravitational lens.Total intensity or magnification µ is normalized, and is plotted as a function of angle on the source plane β.The gravitational point mass and source are aligned at β = 0.A slight bulging effect is noticeable in the σ = 0.25 light curve due to the small source size moving moving across regions with a higher density of points on the source plane.

Figure 5 .
Figure 5. Convergence of the numerical ray-traced solution as the Gaussian source brightness becomes more point-like (the standard deviation σ (Equation8) becomes smaller, shades of red lines) to the analytic solution (black dashed line) for a divergent point lens.Total intensity or magnification µ is normalized, and is plotted as a function of angle on the source plane β.The divergent point lens and source are aligned at β = 0.A slight ringing effect is noticeable in the σ = 0.25 and σ = 0.5 light curves due to the small source size moving in gaps between points on the source plane.
Clegg et al. (1998) presented a classic explanation for the phenomenon of ESEs.In their model of a plasma lens,Clegg et al. (1998) considered a Gaussian-distributed electron column density profile, which produced the notable caustic spikes and flux drop that are characteristic of ESEs.Here, theClegg et al. (1998) plasma lens model serves as inspiration for a full three-dimensional plasma lens with Gaussian-distributed electron density.In dimensionless form the electron volume density for theClegg et al. (1998)-inspired spherical Gaussian plasma lens is

Figure 6 .
Figure 6.Snapshot of a source with Gaussian-distributed brightness traversing horizontally across the source plane (left) and the resulting image on the observer plane (right) for a Gaussian plasma lens.Multipleimaging results in an apparent stretching of the source on the image plane.The parameters are typical of a pulsar ESE with the lens placed halfway between observer and source: D s = 1 kpc, D d = 0.5 kpc, ρ H = 10 −20 g cm −3 , with a lens of length scale r 0 = σ = 10 AU.The colour bar indicates the normalized intensity.

Figure 7 .
Figure 7. Complete light curve as the source traverses horizontally from the left to right edge of the source plane (Figure 6, left panel).

Figure 8 .
Figure 8. Magnification map of the volumetric Gaussian lens example described in the text.Magnification is calculated by the ratio of Delaunay triangle areas on the image and source plane (Section 4).

Figure 9 .
Figure 9.Light curves for the Grafton et al. (2023) filament reimagined as as plasma lens for inclination angles i = [0, π/6, π/4, π/3] (shades of red).The source travels from left to right in the source plane, perpendicular to the long-axis of the filament when i = 0, and takes the same path for each inclination angle.The lensing system is set up with the filament half way between the observer and source plane that D s = 1 kpc, D d = 0.5 kpc, with length scale r 0 = 10 AU, core density ρ 0 = 10 −21 g cm −3 , and Γz = 10 −3 .