Mapping Plasma Lenses

Dynamic spectral mapping of interstellar plasma lenses


Compact radio sources sometimes exhibit intervals of large, rapid changes in their flux-density, due to lensing by interstellar plasma crossing the line-of-sight. A novel survey program has made it possible to discover these “Extreme Scattering Events” (ESEs) in real time, resulting in a high-quality dynamic spectrum of an ESE observed in PKS 1939–315. Here we present a method for determining the column-density profile of a plasma lens, given only the dynamic radio spectrum of the lensed source, under the assumption that the lens is either axisymmetric or totally anisotropic. Our technique relies on the known, strong frequency dependence of the plasma refractive index in order to determine how points in the dynamic spectrum map to positions on the lens. We apply our method to high-frequency (4.2-10.8 GHz) data from the Australia Telescope Compact Array of the PKS 1939–315 ESE. The derived electron column-density profiles are very similar for the two geometries we consider, and both yield a good visual match to the data. However, the fit residuals are substantially above the noise level, and deficiencies are evident when we compare the predictions of our model to lower-frequency (1.6-3.1 GHz) data on the same ESE, thus motivating future development of more sophisticated inversion techniques.

Subject headings:
ISM: general – ISM: structure – scattering – gravitational lensing – methods: data analysis

1. Introduction

Extreme Scattering Events (ESEs) were discovered almost thirty years ago (Fiedler et al., 1987), yet their cause remains obscure. That is in part because it is difficult to construct acceptable physical models for the phenomenon. Refraction by ionised interstellar gas appears to be the mechanism responsible for the variations (Fiedler et al., 1987; Romani et al., 1987), but the inferred properties of the plasma lenses are problematic. In particular the very high gas pressures, confined in very small regions, are challenging to interpret. And the challenge is heightened by the fact that the lenses are commonplace in the Galaxy, with local to the Sun.

Despite their large numbers, the small size of the lenses (a few AU) means that the event rate is actually quite low; Fiedler et al. (1994) estimated roughly one in 200 sources undergoing an ESE at any instant, but using more restrictive criteria other authors have suggested a fraction as small as one in 2000 (e.g. Karastergiou and Walker 2010). That is the other reason why ESEs are still not understood: we have had very few examples to study. However, that problem is now being addressed with a new survey program which recognises events in progress, and is thus able to follow-up each event with intensive multi-wavelength monitoring (Bannister et al., 2015). As a result of that program, we now have a high-quality dynamic radio spectrum, from the Australia Telescope Compact Array (ATCA), for an ESE towards the source PKS 1939–315. Our data span a factor in wavelength, corresponding to a factor in the focal length of the lens, and thus provide strong constraints on viable lens models. This paper details our analysis of that dynamic spectrum.

Previous analyses of ESE light-curves have proceeded by fitting to models based on analytic lens profiles, notably the Gaussian form (e.g. Clegg et al. 1998). Here, however, we attempt to invert the dynamic spectrum of the lensed source and thus discover the column-density profile of the plasma lens without prejudice. To do so we restrict attention to two possible geometries: axisymmetric lenses, and extremely anisotropic lenses (i.e. lenses which refract in only one direction). Most of the physical scenarios which have been proposed for ESEs can be approximated with models which fall into one of these two categories (Henriksen and Widrow, 1995; Walker and Wardle, 1998; Walker, 2007; Pen and King, 2012). And for these geometries the inversion is strongly over-constrained because the profile we seek is one-dimensional, whereas the data are two-dimensional.

This paper is organised as follows. In section 2 we recap the basics of lensing theory, using geometric optics and the thin lens approximation, and apply the theory to the case of a plasma lens. Section 3 then presents our approach to reconstructing the plasma density in the lens, using the well-defined scalings of the plasma refraction angle and image magnification with wavelength. Section 4 applies these ideas to our data on PKS 1939–315 and presents our main results, and section 5 discusses their interpretation.

2. General lens model

In the geometric optics approximation (see e.g., Schneider, Ehlers and Falco 1992 where most of the material in this section can be found), the optical effects of a plasma lens can be encapsulated in the “lens equation”


which relates the angular position of the source, , to the angular position of its image, , via the “deflection angle,” , and the distances from lens-to-source (), and observer-to-source (). The deflection angle is actually the negative of the refraction angle introduced by the lens, with the sign chosen for consistency with the majority of the gravitational lensing literature. In general there may be more than one image of the source, corresponding to the case where there are multiple solutions to equation (1) for a given value of .

For a thin lens, the refraction can be characterised in terms of a dimensionless scalar quantity, , referred to as the “lens potential”, such that


where is the gradient with respect to angular position in the lens plane (Schneider, 1984). Our lens equation then takes on the simple form


and provides a complete description of the lensing behaviour. We will see, below, that in the case of plasma lensing the potential function is proportional to the electron column-density profile of the lens, .

For sufficiently small sources, and we assume a point source for the rest of this paper, the flux of an image is determined by the magnification of the lens at the location , because intensity is conserved along rays if there is no absorption or gain introduced by the lens. Here the term “magnification” means the ratio of solid-angular size of the image to the solid-angular size of the source. Magnification is thus calculated from the differential properties of the lens mapping, i.e. the relationship of to , where the latter is the locus that describes the perimeter of the source. In general the image is distorted, with different amounts of stretching in different directions, and one evaluates the image magnification from the Jacobian matrix, (Schneider, 1985):




where the convergence, , and shear, , are invariants of ,


familiar from optics. The former describes the isotropic deformation of a small light beam due to the lens, whereas the latter characterises stretching of the beam along one axis relative to the other. Equations (1-6) are quite general and can be used to describe various types of lensing in the geometric optics approximation.

2.1. Potential function of a plasma lens

The additional phase, , impressed by the lens on a wave which has propagated through it, is given by


where is the refractive index of the lens, and we have assumed propagation in the direction. Here we are using the thin lens approximation, in which one neglects the effect of deflection of the wave during passage through the lens, and , independent of the angle of incidence.

As noted above, the angle is the negative of the refraction angle introduced by the lens. By definition the wavevector is normal to the wavefront, i.e. normal to surfaces of constant phase, and for small deflections the refraction angle is thus given by the ratio of the transverse phase gradient to the longitudinal one, so


Here the operator is a two-dimensional vector quantity, and the subscript emphasises that this gradient is taken with respect to the position , where is the distance of the lens from the observer. We can relate this operator to the angular gradient via .

For an electromagnetic wave of angular frequency , the refractive index of a cold plasma is (e.g., Stix 1992)


where the plasma frequency, , is assumed to be such that . The plasma frequency is given in terms of the electron density, , by


Here and are the charge and mass of the electron, respectively, and is the permittivity of free-space.

Using equations (7-9) we thus obtain


where is the classical radius of the electron. Combining this result with equation (7) then yields the potential function




and is an arbitrarily chosen reference wavelength.

Note that the deflection angle and its derivatives are all proportional to :


Thus, a definite scaling should hold for the inverse magnification2 of an image measured at a fixed position :


We now discuss the possibility of such measurements.

3. Interpreting data

3.1. Characteristics

In the limit where , the lens mapping of equation (3) reduces to the identity , so there is only one image and the magnification of that image is . As increases there is still only one image, formed close to , but the deflection angle becomes greater, and the magnification deviates further from unity. The single-image/high-frequency regime is convenient for two reasons. First, if the lens is so weak that there is only one image present then there is no ambiguity of interpretation: the current source position corresponds to a unique image position, and the measured flux of the source tells us about the magnification – and via equation (4) the curvatures of the lens potential – at that point. Furthermore the image has positive parity, i.e. (Burke, 1981), and it is not necessary to take the absolute value in equation (5). Secondly, if the image magnification is not far from unity then the influence of source structure on the observed fluxes is less likely to be important, and a point-source model can be used in the first instance.

In the absence of imaging data – i.e. data which reveal the angular structure of the lensed source – the argument  of the potential is not observationally accessible and needs to be reconstructed along with the potential function from equations (3, 5). It is therefore convenient to consider itself as the primary dependent variable when reconstructing the lens. The deflection angle can then be obtained by interpolating between the pairs and the locus of clearly defines the portion of the lens plane constrained by the available data. The potential function and electron column density are subsequently found by integrating the deflection angle. Relegating to the status of a dependent variable has the considerable practical advantage of bypassing the computationally costly solution of the lens equation. It can also provide an insight into the qualitative structure of the data, and a simple visualisation of the predictive capabilities of the model.

The latter is made possible by identifying the locus of points in the data that correspond to the same position in the lens plane. From the lens equation (3) with the deflection angle (14), all along the line


in the domain of our independent variable project to the same dependent variable in the image plane for all . In order to compute , we need to know at a single point on the line . Substituting the corresponding value of into equation (16), we obtain


Thus, knowing the mapping at one particular position and wavelength determines the solution along the entire line (17) in its domain, and by construction. By analogy with the theory of partial differential equations this line can be referred to as a “characteristic”. Later, in Section 4 (Figure 2), we will see examples of systems of characteristic lines appropriate to particular lens models.

The reference wavelengths in equations (12-16, 17) need not be equal and it is convenient to use , which we assume in the following.

3.2. One-dimensional problems

Two-dimensional (epoch and wavelength) data of the dynamic spectrum over-constrain models of electron density that depend on only one spatial coordinate, such as the two cases we consider – extremely anisotropic and axially symmetric lens models. Characteristics are useful in defining which multiple data points constrain the same location in the lens.

Extremely anisotropic case

When the electron density depends on only one of the two Cartesian coordinates in the lens plane, , the deflection angle is always parallel to it and the source position is unaffected in the other direction:


we will only consider the component, dropping the subscript. Using equation (3), the image magnification is


and therefore, on a characteristic, where the argument of the derivative in the denominator is constant, the derivatives, which are measurable as , should give the same result when rescaled to the reference wavelength:


Comparing equation (20) to equation (15) we also establish the relationship between the optical scalars in the extremely anisotropic case:


which can be used as a consistency check if can be obtained independently – e.g., by fitting a biquadratic polynomial (15) to the magnifications along the characteristic defined by .

Axially symmetric case

In complete axial symmetry, , the deflection is always in the radial direction:


The relation between the observable and is slightly different as one has to account for the curvature of the polar coordinate system when computing magnification:


But as the lens equation (22) is formally equivalent to equation (18), the scaling (20) remains valid along the characteristic in the axisymmetric case. A slight complication of the to conversion by the factor presents no difficulty in practice because, in order to consider variation of the derivative along a characteristic, the characteristic’s needs to be specified in the first place. For the axisymmetric case, the optical scalars are given by


Formalisation of one-dimensional problems

By definition of a characteristic line, is constant along the line. We therefore require that the right-hand side of equation (20) – with computed as or in the extremely anisotropic and axisymmetric cases, respectively – is constant on characteristics:


This relation expresses the over-constrainedness of the model by the data in the one-dimensional case.

3.3. Charting the characteristic set

The actual characteristic passing through a given point is not known in advance and one can envisage a number of approaches to correctly drawing the family of these lines through the data. The latter task is equivalent to the lens reconstruction problem since each characteristic corresponds to a measurement and can be converted into by interpolating on pairs. There are two constraints to guide the solution process: equation (25) and equations (19, 23); the former expressing the expected behaviour of the data along the characteristics and the latter across them. The relative weights given to the two in a particular combination define the various ways in which the solution could be attempted.

Condition (25) can be explicitly resolved with respect to to define the characteristic suggested by the data in the neighbourhood of a given point:


but such would not generally be constant along characteristics, rendering this relation of little practical use. The reason is that it implies a consistency condition on the data themselves which, even if met by the underlying signal, will be violated by the random and systematic errors which are inevitable in the measurement process.

A milder version of this method would be to select, among all the possible lines (17) passing through a given point , the one on which the gradient (20) is least variable. However, this approach makes no attempt to satisfy equations (19, 23) and the resulting does not have to be smooth, continuous or even self-consistent. We have found that both of these “morphological methods” performed poorly on our data.

The approach that has so far proven most useful, and the one which we use to analyse our data in Section 4, is as follows. We first select a pivot point, , through which we would like to draw a characteristic. With an assumed value of the corresponding position in the lens plane, , the trajectory of the characteristic in the data domain, i.e. , is then fixed by equation (17). Data points which lie on or near this trajectory give us information on the magnification along the characteristic. Scaling to the reference wavelength, we then use these magnifications to form the average as our estimate of the gradient of the lens mapping at the position . We then select a new pivot point, , lying close to, but not on the first characteristic. The corresponding position in the lens plane, , is given by


where is the gradient of the first characteristic evaluated at , and is the distance along the axis between the new pivot point and the initial characteristic. Knowing fixes the entire trajectory of the new characteristic, and so we can use our data to form an estimate of the gradient, , for the new characteristic. This process is then iterated, allowing us to cover the data domain with an entire family of characteristics. Once we fix the set of pivot points, the foregoing procedure furnishes us with a unique lens model for each assumed value of .

As we have no prior knowledge of the appropriate value of , we must choose a suitable figure by determining which of the resulting lens models we prefer. This is readily achieved, e.g. by introducing a figure of demerit that measures the root-mean-square difference between the model predictions and the data. The procedures just described are summarised in pseudocode as Algorithm 1.

Because the characteristic lines are never parallel to the axis, it is convenient to choose the same wavelength for all pivot points. We chose all pivots to be at 8.4 GHz. And because deflection angles are small at high frequencies, and our data are not uniformly sampled in time, we placed characteristics in one-to-one correspondence with the sampled epochs, . Over much of the 4 cm observing band, this choice gives us a close correspondence between the locations of the data samples and the locations of the characteristics.

Data: Magnification
Result: Deflection curve
Construct interpolating function ;
Select a set of characteristic pivot points ;
Select a set of initial values to try;
for  do
       Set first characteristic parameter ;
       Compute average scaled derivative from
           magnification data interpolated onto the
           characteristic ;
       for k=2..K do
             Compute distance to next characteristic;
             Rescale to ;
             Compute parameter ;
             Compute average scaled derivative along
                 characteristic ;
       end for
      Compute model along all characteristics;
       Interpolate onto data sampling points ;
       Compute demerit ;
       Store demerit, characteristic parameter set ;
end for
Identify initial value that results in optimal demerit;
interpolates optimal , .
Algorithm 1 Averaging along trial characteristic set. The derivative is obtained from the magnification using equation (19) or equation (23), for the extremely anisotropic or axisymmetric geometry, respectively, and scales with wavelength according to equation (20).

A clear advantage of this method is that there is just a single variable with respect to which the model is to be optimised, and yet it takes account of both equations (25) and (19, 23), even if only in an average sense for the latter. It also performs well in practice, producing reasonable models and doing so stably with respect to variations in the spectral and kinematic models or readjustment of weighting schemes. A serious drawback, however, is that errors in determining the properties of one characteristic propagate into all subsequent characteristics in the inner loop of Algorithm 1. Of particular concern with the real data, which inevitably contain gaps in the coverage, is the possibility of missing episodes of strong (de-)magnification, which strongly deforms the paths on which the derivatives are averaged.

A simple way to mitigate problems associated with gaps in the coverage is to break the dataset into a few separate pieces at suspect positions – particularly in the longer gaps in the data coverage, where a missed high-maginfication event might lurk. The algorithm is then run on each of the pieces independently. If the results align well, it is unlikely that an important event has been missed and the data can be restitched together for a more reliable reconstruction. If the algorithm returns considerable offset between the end of one reconstructed piece and the beginning of the next, it is a good indication an important event in the light curve has been missed.

4. Application to ATCA data on PKS 1939–315

4.1. Kinematic and spectral model

The algorithm described above deals with the angular position and magnification of the source whereas the available data come in the form of a dynamic spectrum . To convert the latter to the former, models of the instrinsic spectrum of the source, , and of its motion relative to the Earth-lens system, , have to be assumed. We aimed to keep both models as simple as practical, restricting ourselves to a power-law spectrum and rectilinear motion:




Here the effective transverse velocity, , is given in terms of the source, lens and Earth velocities (, , and , respectively) by


(Cordes and Rickett, 1998). In the case of a quasar seen through a plasma lens in our own Galaxy, , so . It is convenient to choose the epoch so that is the impact parameter in the axisymmetric case: .

In the extremely anisotropic case, only the component of equation (29) is relevant and the photometry is insensitive to changes in or . Further, only affects the overall scale of the angular coordinates by setting the time-to-angular-position conversion factor. That too is irrelevant for the photometry, which is determined by the dimensionless derivatives. As a result, only the normalisation, , and slope, , of the intrinsic spectrum need to be specified, and those choices can be optimised, while the angular positions are effectively measured in time units – i.e., .

In the case of axial symmetry the observables depend on the magnitude of ,


and it is now only the overall angular scale of the problem, that remains free. Thus, the parameter inventory of this case includes both the spectral and kinematic portions. Optimisation in this four-dimensional space is slow, so we only optimised with respect to the impact parameter and epoch, , adopting the best-fit spectral model of the extremely anisotropic case. The impact parameter and epoch of closest approach, as determined by our optimisation, are given in Table 1 along with the parameters of the spectral model.

, mJy , days , MJD , %
Table 1Best fit spectral and kinematic parameters for PKS 1939–315.

4.2. Implementation details

We applied the methods presented in section 3.3 to our Australia Telescope Compact Array data on the ESE towards PKS 1939–315. This event was discovered via a novel method utilising the wide bandwidth of the Compact Array Broadband Backend (CABB, Wilson et al. 2011). The method of discovery and the follow-up data are described by Bannister et al. (2015). The data consist of flux density measurements taken in 16 cm (), and 4 cm (), bands with -wide spectral channels.3 The source was observed every few days after the ESE was identified in June 2014 until March 2015. Gaps in this coverage lasting up to a month are present, particularly in the 16 cm band. The data are presented in Figure 1.

Figure 1.— Observed flux density of PKS 1939-315 (in , as per the colour bar) as a function of Modified Julian Date and frequency. The horizontal, white strip is the gap between the and bands in our observing configuration. The upper (lower) set of tickmarks in this strip denote the epochs of observation in the () band. Between observing epochs, we use linear interpolation to determine appropriate flux densities. The time axis in this figure is related to the apparent source position, , via the adopted kinematic model.

We tried the different methods described in section 3.3 when reconstructing the characteristic system of the lens mapping , for both extremely anisotropic and axisymmetric cases, but only the method of Algorithm 1 proved useful in practice. We tested for missing high- or de-magnification events in gaps in temporal coverage by breaking the data into separate pieces, but we found no conclusive evidence for any. The test instead revealed a small but consistent lag of change in compared to the predictions of , as if the latter was biased low. We interpret the bias as due to the presence of a steady component that is not subject to lensing, either because of its large size, or because it is substantially offset from the lensed component of the source. To counter the bias we introduced into our models a parameter, , to describe the unlensed fraction of the source flux. To keep the spectral model simple we adopted the same unlensed fraction at all frequencies; the optimal values found were for the extremely anisotropic and for axisymmetric geometries. We thus employ a total of three model parameters to describe the source.

In the vicinity of caustics, a point-source model yields arbitrarily large fluxes and the root-mean-square flux difference is overwhelmingly influenced by these regions. Root-mean-square flux difference is therefore not a good statistic for assessing the goodness-of-fit of a model. A better choice of demerit is to use the root-mean-square difference in the inverse-flux:


which we used in our optimisation. We also quote the r.m.s. of the flux density residuals, to permit comparison with measurement noise levels. Although only a single number, , is freely adjusted during the optimisation in Algorithm 1, we include in our parameter count the total number of characteristics in the model, which is 86. We do so because for each characteristic the gradient of the lens mapping is chosen so as to match, in an average sense, the data close to that characteristic. This also allows us to compare models with varying resolution and regardless of how they were obtained.

The refractive index of a plasma increases strongly with wavelength, as per equation (14), leading to multiple refracted images, and to diffractive scattering at sufficiently low frequencies. Moreover, the angular size of compact radio quasars tends to increase with wavelength. These considerations suggest that the methods discussed above are best suited to short radio wavelengths. We therefore focused our modelling effort exclusively on the 4 cm band, with the 16 cm band data being used as an independent check on the model. We present our model predictions in the 16 cm band, below, for illustration only; we compare them to the data and discuss the discrepancies in Section 5.

Figure 2.— Characteristic curves (upper panels: a, b), and dynamic spectra (lower panels: c, d), for the optimal extremely anisotropic model (left panels: a, c), and the optimal axisymmetric model (right panels: b, d). The characteristic curves are plotted in blue for positive parity and red for negative; the red-blue boundary exhibits the caustic structure of these models. In order to permit a direct visual comparison with figure 1, the model dynamic spectra are constructed with identical time and frequency sampling as our data, and use the same colour scale as figure 1. At a frequency of GHz, the locations of the characteristic curves coincide with the cm band observing epochs.

The characteristic lines used in modelling cannot coincide with the data at all frequencies. To estimate flux densities on characteristic points that fall in between the data sampling locations (and vice versa), we used linear interpolation in coordinates4. When computing the average of equation (20), which depends non-linearly on the observable flux densities, , we weighted the estimates at each wavelength in proportion to the (squared) derivative of with respect to , so as to avoid biasing the average by a few estimates with a low denominator. This presumes uniform uncertainty in along the spectrum, which should be a good approximation for the thermal noise in our data.

For computational reasons, we averaged the data over spectral channels, yielding a regular spectral grid sampled at 100 MHz intervals. We chose all the characteristic pivot points at , and placed them at positions corresponding to the 86 epochs available in the 4 cm band between June 2014 and March 2015. This choice minimises the influence of interpolation, because bending angles are small at high radio-frequencies — as can be seen from the characteristic curves plotted in Figure 2.

4.3. Modelling results

The results of fitting for both the extremely anisotropic and axisymmetric geometries are quite similar, both qualitatively and in terms of the demerit figures, and we cannot decide between these geometries. Table 2 summarises the performance of the resulting models and Figure 2 shows the predicted dynamic spectra. Both models produce a reasonable fit to the data in the 4 cm band. Nevertheless, the uncertainties on our data are significantly lower than the r.m.s. flux residual between model and data, 5 and therefore there is considerable scope for improvement in the models. At lower frequencies the residuals are much worse and exceed the average flux density of the source itself – a result of the very high magnifications in the vicinity of the model caustics.

Figure 3.— Deflection angles of the two optimal models shown for the reference frequency . The horizontal axis is position in the lens plane, , with . The natural modelling units are shown on the left and bottom axes, whereas the right and top axes assume the effective velocity and distance of , for conversion to physical units. The latter scale as and on the horizontal and vertical axes, respectively. The origin coincides with the position of the source at the reference epoch MJD56800 for the anisotropic model, and the symmetry centre for the axisymmetric model.
Anisotropic Axisymmetric
#model parameters 86+3 86+3+2
demerit (), 4 cm
r.m.s. residual flux (mJy), 4 cm
demerit (), 16 cm
r.m.s. residual flux (mJy), 16 cm
Table 2Parameters of the optimal models
Figure 4.— Panel (a) shows the potential function, , obtained by integrating the deflection curves in Figure 3. Panel (b) shows the same curves with the mean gradient subtracted, in order to show more clearly the lumps and bumps which are present. The horizontal axis is position in the lens plane, as in Figure 3. The potential is presented in the natural modelling units on the left-hand vertical axis. The right-hand axis shows the physical units corresponding to and and scales as .

Figure 3 presents the deflection curves of the two characteristic sets. The curves are presented both in the “modelling units” of days (i.e., with and ), and for illustration we have also converted to physical units by assuming specific values for the effective velocity and distance: and , respectively. Integrating the deflection curves via equations (2, 12) we obtain the potential function, , which is proportional to the electron column density (13).

Figure 4 presents the inferred potential on top, and the bottom panel shows the same estimates with the mean gradient of removed (so as to reveal the lumps and bumps which are superposed on the trend). The gradient does not affect magnification and would be impossible to measure at a single frequency. It does affect the alignment between light-curve features at different wavelengths, making it detectable with spectral data. The mean gradients that we infer through averaging deflection angles in Figure 3 are


where the upper and lower figures are for the extremely anisotropic and axisymmetric geometries, respectively.

Our method is completely insensitive to the density zero point so readers are free to add any constant vertical offset they like to the points in figure 4. On the other hand the curvature and gradient of the density directly affect the light curves and their alignment across the frequency dimension, so those quantities are fixed by the modelling at every plotted point in figure 4.

The uncertainties in mean-gradient and mean-curvature may be of interest. By varying the parameter we obtain models which have a different mean gradient from our optimum model, but which necessarily have a greater demerit. Because the uncertainties on our data are primarily systematic (), we can only give rough estimates of the uncertainty associated with model parameters; we adopt uncertainty intervals defined by the mean-square residual being no larger than its minimum value plus . With this criterion we find that the uncertainty on the mean-gradient in the lens is % around the value given in equation (4.3).

Even a small change in the mean curvature away from our best lens model would, by itself, greatly increase the residuals of the fit, because of the change in the time-averaged spectrum of the lensed source. However, the inferred mean curvature of the lens is, to some extent, degenerate with the shape of the intrinsic source spectrum, and the intrinsic spectrum is not precisely known. We therefore modified our source spectrum model in such a way as to precisely mimic the effect of magnification by a plasma lens, and then determined the mean convergence of the resulting lens model. As expected, the lens model tries to compensate for the change in the spectral model in such a way that the lensed source spectrum is largely unchanged. By varying the degree to which the intrinsic spectrum is modified, and evaluating the mean-square residuals of the resulting lens models we were thus able to estimate the uncertainty in the mean-convergence of the lens. We determined that the data permit an additive constant convergence of up to . In other words, the zero-point in figure 5 is uncertain by only . The sensitivity of our data to a constant offset in the lens curvature stems, in part, from the large region of the lens-plane which we have sampled: a range in of order 250 days. Over this interval a constant curvature of only changes the gradient by about 5 days, which is roughly twice as large as the mean gradient in our best models — see figures 3 and 4.

We emphasise that these constraints are strictly local and do not depend on the behaviour of the column-density at large distances from the line-of-sight — there is no analogue of the “external shear” which is encountered in gravitational lensing. Nor is there an analogue of gravitational lensing’s “mass-sheet degeneracy” — that degeneracy is broken by the strong frequency-dependence of the plasma refractive index. For example: the large negative curvature in (hence ) around day 60 in Figure 4 is responsible for the period of strong demagnification around MJD 56850 in Figure 1.

5. Discussion

Even with an ideal inversion algorithm, the dynamic spectrum alone cannot provide a complete characterisation of a plasma lens. Such an inversion returns , and the scaling factors required to obtain plasma column-density as a function of position (see section 2) must be determined by other means. In common with previous authors we can, of course, simply assume values for the effective transverse velocity and lens distance, in order to fix those conversion factors. And, naturally, if we choose values similar to those used by previous authors then we obtain similar properties for the lenses themselves, because in all cases the basic observational requirement is to obtain magnification variations for a few months. Thus if we adopt and , then we infer lensing by a structure of size with electron column-density gradient . For an axisymmetric lens, which presumably arises from a structure that has spherical symmetry, the column-density gradient is expected to be comparable to the volume density of electrons in the lens, implying tiny lenses which are highly over-pressured. For highly anistropic lenses it is not obvious what the appropriate three-dimensional structure is, so the relationship between and remains unspecified. Many authors have argued in favour of sheet-like structures seen edge-on, as a means of moderating the required electron pressure (Romani et al., 1987; Goldreich and Sridhar, 2006; Pen and King, 2012; Pen and Levin, 2014).

In addition to ESEs, two other radio-wave scattering phenomena appear to require a large population of tiny structures incorporating ionised gas with large gradients in . These phenomena are the Intra-Day Variability (IDV) of some compact radio quasars (e.g. Kedziora-Chudczer et al. 1997; Dennett-Thorpe & de Bruyn 2000; Bignall et al. 2003; Lovell et al. 2008), and pulsar parabolic arcs (Stinebring et al., 2001; Cordes et al., 2006). These two phenomena can be plausibly attributed to a single type of scattering structure, differing only in the nature of the background radio source (Tuntsov et al., 2013). There is good evidence that these structures are often very highly anisotropic in their scattering (e.g. Walker et al. 2004; Cordes et al. 2006; Walker et al. 2009; Brisken et al. 2010), which suggests that strong magnetic fields may be present. In the modelling of PKS 1939–315 reported in this paper we were unable to distinguish between axisymmetric and highly-anistropic lenses. However, it is notable that our best-fit axisymmetric model demanded a very small impact parameter for the event (consistent with zero), so that the relative orientation of the effective velocity vector and the plasma density contours is essentially constant during our monitoring interval. Put another way: when we ask for an axisymmetric model, the data prefer us to choose one which looks very much like an extremely anisotropic model. This deduction, made from fitting the dynamic spectrum alone, is consistent with the astrometric shifts of PKS 1939–315, measured with the Very Long Baseline Array (Bannister et al., 2015), which are all consistent with refraction in the same direction on the sky. We note, however, that the interpretation of the observed astrometric shifts is unclear at present (see Bannister et al. 2015).

One further aspect of our data is worth drawing attention to. In the dynamic spectrum of PKS 1939–315, there appear to be flux variations at all times. Thus, although there is an interval of outstandingly strong variability in our dynamic spectrum of PKS 1939–315 – i.e. the two months either side of MJD56850, which triggered our follow-up observations (Bannister et al., 2015) – that interval should perhaps not be thought of as an isolated event. Rather, it might be better thought of as a large example of the ongoing fluctuations which are currently manifest in this particular souce.

Although the foregoing points are suggestive, not conclusive, they are all consistent with the Extreme Scattering “Event” in PKS 1939–315 being closely related to the episodic phenomenon of Intra-Day Variability in radio quasars, and thus also to the pulsar parabolic arcs. It has previously been noted that in pulsar parabolic arcs one sometimes sees deflection angles consistent with a lens that is strong enough to explain the large magnification fluctuations observed in ESEs (Hill et al, 2005).

A valuable, independent check on the quality of our modelling is provided by our ATCA 16 cm band data on PKS 1939–315, which were not utilised in constructing our model plasma lens profiles. The lensing effects of a plasma increase with wavelength, and caustics are expected to form at some point. Comparison of Figure 2 with Figure 1 immediately reveals the deficiencies of our models, particularly in the poor correspondence between the locations of model caustics and local peaks in observed flux density. Nor do the peaks in the data appear strong and sharp, as expected for caustics. These qualitative deficiencies are quantitatively reinforced by the figures of demerit and r.m.s. flux residuals given in table 2. The poor performance of our models at low-frequencies is not surprising, for two reasons. First, because the optical-depth to synchrotron self-absorption decreases with frequency, radio quasars are most compact at high-frequencies, and thus our point-source approximation is a poorer approximation at low-frequencies. Indeed, VLBA observations reveal (Bannister et al., 2015) that the resolved jet component contributes only a few per cent of the total flux density at but at and the resolved component dominates the compact core. Although the foregoing point is not necessarily responsible for errors in the location of the caustics, the figures of demerit and r.m.s. flux residuals are both strongly influenced by the presence of caustics when the moderating influence of source-size vanishes. Second, our treatment is based on geometric optics and thus excludes the effects of diffraction by small-scale column-density structures in the lens. Diffraction is also expected to be more important at low frequencies because the phase profile of the lens scales in proportion to the wavelength (11), and because the Fresnel scale, which separates refractive and diffractive regimes, increases with wavelength. Unfortunately, non-zero source size and diffraction both appear to be incompatible with the core ideas of this paper – i.e. characteristics, and a well-defined scaling of magnification with wavelength – so the exclusion of these physical effects is a fundamental limitation of our method. On the other hand forward modelling that includes these phenomena is feasible, which suggests that parameterisation of the lens and source structures, followed by optimisation, should be a profitable avenue to explore.

Our method of determining the lens profile provides built-in consistency checks (21, 24) that can be applied to the reconstructed characteristic sets, by estimating the optical scalars from the fit to equation (15) along each characteristic. We performed these checks, fitting for the coefficients of and in equation (15), and the models largely satisfy their consistency relations. But the uncertainty and level of degeneracy between the two coefficients preclude us from drawing any further conclusions. The uncertainties are in fact so large that the optical scalars inferred for the axisymmetric model easily satisfy the consistency check for the anistropic model.

Figure 5.— Beam convergence (at the reference frequency ) of our optimal extremely anisotropic model computed as the average of equation (21) along each characteristic. The horizontal axis is the position in the lens plane, as in Figure 3. The convergence is clearly negative in the course of the strong demagnification event. The error bars show the r.m.s. variation of the individual estimates along the characteristic.

The results of our modelling are nevertheless very definite regarding the nature of the lens responsible for the period of strongest demagnification. As shown in Figure 5, it was caused by a diverging lens, with . Pen and King (2012) argued that the intervals of strong demagnification observed during ESEs are naturally explained by a converging lens that is over-focused. Our lens inversion is predicated on the assumption that only a single image is present at high frequencies, and does not directly address the circumstance of an over-focused lens (which is necessarily in the caustic regime). However, towards the upper end of the frequency range of our data (Figure 1) the fractional variations in flux are both small and decreasing, making it difficult to sustain the idea that this lens is over-focused.

Attempting to extend our approach, which identifies characteristics by minimising the deviation from equation (15), to the case of a general two-dimensional lens incurs the penalty of a broadened interpretation of both the solution and the characteristic lines. It can be shown that the general reconstruction problem in 2D does not have a unique solution, and therefore one is obliged to look for particular classes of solution or their representatives. Likewise, a one-dimensional line (16) in the three-dimensional model space would not generally lie on the surface defined by the two-dimensional data but instead only intersects this surface. However, the lines on this surface that minimise the deviation from equation (15) do exist, and by construction they project to intervals with little variation in the optical scalars. This makes it possible, in principle, to constrain the physical properties of the screen in both transverse coordinates.

6. Conclusions

The advent of new, high quality data on plasma lensing events motivates the development of novel techniques which exploit the information in those data to constrain the lens properties. We have presented one such method, which takes a dynamic spectrum as input and yields the plasma column-density profile as output. In so far as possible, our approach avoids the proliferation of model freedoms: we use geometric optics, the thin-screen approximation, a point-source model, and a lens with reduced dimensionality. These restrictions are all of questionable validity and one or more may need to be relaxed in future. Nevertheless, restricting the solution space in this way permits a straightforward and rapid exploration of the lens properties, providing immediate insights into lens structure and a firm foundation for more sophisticated modelling.

Applying our method to ATCA data on the ESE in PKS 1939–315 we recovered profiles of the lens responsible for this event in the case of either axisymmetric or extremely anisotropic lens geometries. The two profiles share many features, including a large mean gradient and a column-density peak associated with the region of strongest lensing. Both our models provide a good visual match to the high frequency data (4.2-10.8 GHz), which were used as input, but neither did a good job of predicting the lower-frequency data (1.6-3.1 GHz), which were used only as a test of the models. There is, therefore, considerable scope for improvement in the inversion procedure.

The Australia Telescope Compact Array is part of the Australia Telescope National Facility, which is funded by the Commonwealth of Australia for operation as a National Facility, and managed by CSIRO. We especially thank the scheduler, Phil Edwards, for making the regular ATCA observations possible. We have benefitted greatly from the guidance of Ron Ekers on a variety of issues relating to this work.


  1. slugcomment: Submitted to ApJ 10th December 2015
  2. The minus sign is for images of negative parity.
  3. The original data covered bandwidths slightly larger than the values quoted here. We trimmed the band edges, as the extrema are difficult to calibrate accurately.
  4. Specifically, we interpolated in first – at for data, or along the characteristics for models – and then interpolated in .
  5. With channels of 100 MHz, the thermal noise is only mJy and the uncertainties in the data are dominated by residual calibration errors and source confusion, which we expect to be of order 1% (i.e. mJy).


  1. Bannister K.W., Stevens J., Tuntsov A.V., et al. 2015, submitted
  2. Bignall H. E., Jauncey D. L., Lovell J. E. J., et al. 2003, ApJ, 585, 653
  3. Brisken, W.F., Macquart J.-P., Gao J.J., Rickett B.J., Coles W.A., Deller A.T., Tingay, S.J., and West C.J., 2010, ApJ, 708, 232
  4. Burke W.L., 1981, ApJ,244, L1
  5. Clegg A.W., Fey A.L., Lazio T.J.W., 1998, ApJ, 496, 253
  6. Coles W.A., Kerr M., Shannon R.M., et al., 2015, ApJ, 808, 113
  7. Cordes J.M., Rickett B.J., 1998, ApJ, 507, 846
  8. Cordes J. M., Rickett B. J., Stinebring D. R., Coles W. A., 2006, ApJ, 637, 346
  9. Dennett-Thorpe J., de Bruyn A. G., 2000, ApJ, 529, L65
  10. Fiedler R.L., Dennison B., Johnston K.J., Hewish A., 1987, Nature, 326, 675
  11. Fiedler R.L., Dennison B., Johnston K.J., Waltman E.B., Simon R.S., 1994, ApJ, 430, 581
  12. Goldreich P., Sridhar S., 2006, ApJ, 640, L159
  13. Henriksen R.N., Widrow L.M., 1995, ApJ, 441, 70
  14. Hill A.S., Stinebring D.R., Asplund C.T., et al. 2005, ApJ, 619, L171
  15. Karastergiou A.,Walker M.A., 2010, in Astronomy with Megastructures. Joint Science with the E-ELT and SKA, ed. I. Hook et al. (Crete, Greece: Crete University Press), 127
  16. Kedziora-Chudczer L., Jauncey D. L., Wieringa M. H., et al. 1997, ApJ, 490, L9
  17. Lovell J. E. J., Rickett B. J., Macquart J.-P., et al. 2008, ApJ, 689, 108
  18. Pen U.-L., King L., 2012, MNRASL, 421, 132
  19. Pen U.-L., Levin Y., 2014, MNRASL, 442, 3338
  20. Petroff E., Keith M.J., Johnston S., van Straten W., Shannon R.M., 2013, MNRAS, 435, 1610
  21. Romani R.W., Blandford R.D., Cordes J.M., 1987, Nature, 328, 324
  22. Schneider P., 1984, A&A, 140, 119
  23. Schneider P., 1985, A&A, 143, 413
  24. Schneider P., Ehlers J., Falco E.E., 1992, Gravitational lenses (Springer-Verlag Berlin Heidelberg New York)
  25. Stinebring D. R., McLaughlin M. A., Cordes J. M., et al. 2001, ApJ, 549, L97
  26. Stix T.H., 1992, Waves in plasma (American Institute of Physics, New York)
  27. Tuntsov A.V., Bignall H.E., Walker M.A., 2013, MNRAS, 429, 2562
  28. Walker M.A., 2007, in ASP Conf. Ser. 365, SINS - Small Ionized and Neutral Structures in the Diffuse Interstellar Medium, ed. M. Haverkorn & W. M. Goss (San Francisco, CA: ASP), 229
  29. Walker M. A., de Bruyn A. G., Bignall H. E., 2009, MNRAS, 397, 447
  30. Walker M. A., Melrose D. B., Stinebring D. R., Zhang C. M., 2004, MNRAS, 354, 43
  31. Walker M.A., Wardle M.J., 1998, ApJL, 498, L125
  32. Wilson W.E., Ferris R.H., Axtens P., et al., 2011, MNRAS, 416, 832
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description