Scattering model of scintillation arcs in pulsar secondary spectra

The dynamic spectra of pulsars frequently exhibit diverse interference patterns, often associated with parabolic arcs in the Fourier-transformed (secondary) spectra. Our approach differs from previous ones in two ways: first, we extend beyond the traditional Fresnel-Kirchhoff method by using the Green's function of the Helmholtz equation, i.e.\ we consider spherical waves originating from three dimensional space, not from a two dimensional screen. Secondly, the discrete structures observed in the secondary spectrum result from discrete scatterer configurations, namely plasma concentrations in the interstellar medium, and not from the selection of points by the stationary phase approximation. Through advanced numerical techniques, we model both the dynamic and secondary spectra, providing a comprehensive framework that describes all components of the latter spectra in terms of physical quantities. Additionally, we provide a thorough analytical explanation of the secondary spectrum.


INTRODUCTION
The first observation of radio pulsars goes back to Hewish et al. (1969).The electromagnetic signals of pulsars encounter on their way a varying electron density of the interstellar medium (ISM), resulting in deflection and scattering of the travelling electromagnetic waves.In addition, the relative motion of pulsar, ISM, and the observer's antenna leads to a Doppler shift and time-varying phenomena.The dynamic spectra consist of the frequency resolved pulse sequences observed over a time-span of up to several hours.A two-dimensional (2D) Fourier transform of the dynamic spectra gives the secondary spectra.Stinebring et al. (2001) discovered parabolic arc structures in the secondary spectra, which result from the interference of multiple signal pathways at a specific distance between the pulsar and the observer.The recent catalogue of scintillation arcs compiled by Stinebring et al. (2022) of 22 pulsars shows various structures and contains the basic physical parameters of the pulsar, such as distance and velocity.Walker et al. (2004) developed theoretical descriptions of the parabolic arcs starting from the Fresnel-Kirchhoff integral.Using Monte Carlo methods, Walker et al. (2004) then computed the locations of points in the secondary spectra.Similarly, Cordes et al. (2006) used a thin phase-changing screen approach to study the dynamic and secondary spectra.Since then there have been a variety of arc studies based on observations by different groups, e.g., Hill et al. (2003Hill et al. ( , 2005)); Wang et al. (2005); Bhat et al. (2016); Safutdinov et al. (2017); Wang et al. (2018); Stinebring et al. (2019); Reardon et al. (2020); Rickett et al. (2021); Yao et al. (2021); Chen et al. (2022); McKee et al. (2022).
Here we put forward a different theoretical approach to treat scattering by the ISM using Green's functions.This method is commonly applied to scattering problems in quantum mechanics; see Kramer & Rodríguez (2006) for an application to matter waves originating from a compact source.By solving Helmholtz's equation in Cartesian coordinates using Green's functions we determine the pulsar spectra received after scattering at the interstellar medium (dynamic spectrum) and its 2D Fourier transform with respect to time and frequency by high-precision numerics (secondary spectrum) for a given scattering configuration.In contrast to the Fresnel-Kirchhoff approach, our method enables the determination of the entire spectrum and relates the strengths of the individual components to physical quantities such as the refractive index and the wavenumber.Furthermore, we give a complete analytical description of the secondary spectra.Walker et al. (2004) obtained point-like peaks in the snapshot regime and determined their positions.We considerably extend this analysis by analytically determining also the peak extensions and the intensities.
In the second section we introduce our scattering approach and present our findings for an analytical description of the spectra in the third section.We conclude in the fourth section and relegate to the appendices some detailed explanations and technical details.

SOLUTION OF THE HELMHOLTZ EQUATION
In this section we propose a Green's function method to describe the scattering of pulsar radiation in the ISM.The differences to the Fresnel-Kirchhoff approach are summarised in Appendix A. We consider scattering from an extended plasma cloud (see Fig. 1), described by a region with the scattering potential  (r ′ ) The electron number density   of the plasma cloud determines the plasma frequency   , the refractive index , and : Pulsar signals travel long distances and undergo a dispersion due to the average electron density in the galaxy.This results in pulses where different frequency components arrive at different times.Here, we are interested in the ISM properties affecting the signal on shorter scales compared to the pulsar distance.This allows us to set the dielectric constant of the background to unity, but alternatively, a uniform background could be introduced.Maxwell's equation for the electric field becomes (Schwinger 1998, Eq. (50.35)) The free Green's function reads Within the Born approximation we replace in the integral the electric field by the incoming electric field E inc (r ′ ) and also neglect any change in the polarisation direction by dropping the Hessian in the second term in Eq. ( 3).This can be easily verified considering e.g. a linearly polarised wave possessing solely a non-vanishing  component depending only on the spatial  coordinate.We take the incoming electric field to be a spherical wave emitted from the pulsar where U 0 determines the polarisation direction of the electric field and has units of a voltage.We obtain the electric field from the Green's function where we consider only one interaction with the scattering potential , corresponding to the Born approximation.The electric field in the Born approximation consists of two contributions, first the unscattered component in the absence of any medium, and second the volume integral over the distribution of plasma clouds in the interstellar medium.For dense or compact objects multiple scattering could be included by summing the Born series in terms of the transition matrix.Such processes are neglected in Eq. ( 6).This equation is also used in quantum mechanics to describe the scattering of coherent electron waves at obstacles (Heller 2018, ch. 26).Since we are only interested in the relative contributions of the electromagnetic waves we set |U 0 | = 1 with direction orthogonal to the direction of propagation.The intensity is given by the absolute value squared of the electric field By taking the absolute value, interference terms appear in the exponents related to the free Green's function.The argument of the exponent contains the differences in distances measured in multiples of the wavelength.For the typical pulsar geometry shown in Fig. 1 all path differences are slowly varying functions across the interstellar medium.We introduce coarse-grained integration regions, which result in a collection of three-dimensional clouds of scattering sources across an entire region of the ISM.For simplicity, we consider a Gaussian electron density profile of the th cloud centred around r  of the form where  ,peak denotes the peak electron density.The first order Born approximation requires to evaluate In the last step in the equation above we evaluated the Gaussian integral via a series expansion for  (see appendix C) and introduced the parameter A uniform background density  , along the signal propagation requires to replace the electron density with the local change in density   (r ′ ) → (  (r ′ ) −  , ).Within the Fresnel-Kirchhoff approach in Walker et al. ( 2004), all the contributions in ( 7) and ( 10) are treated on equal footing rendering it impossible to distinguish their individual prefactors.
The coarse-graining of the ISM complements other approaches which focus on the correlation function of scattering from extended sources described in terms of statistical spatial correlation functions, see Tatarski (1961); Coles et al. (2010).

Dynamic spectra
Eq. ( 10) contains the complete description of the electric field and is next evaluated for specific conditions.We consider a coordinate system where the ISM is considered to be at rest and distributed around the origin of the coordinate system, while the pulsar and the observer are moving as shown in Fig. 1.The pulsar is moving with velocity v  and changes position as function of time likewise the observer moves with velocity v  and is located at position The line of sight vector lies along e  and is given by the direct connection of r  (0) and r  (0), the effect of the velocity components along this axis is negligible due to large values of the distance |  | between the pulsar and the ISM plane and the distance |  | between the ISM plane and the observer.Therefore the velocity component of v  along e  is the only relevant one, while for v  the components along e  and e  need to be considered.A dynamic spectrum is obtained by recording  (r Computed spectra are shown in Fig. 2, left panel, obtained from numerically evaluating Eq. ( 7) and the analytical formulae for 25 scattering clouds.Each scattering cloud is assigned the same value of   = 4.5 × 10 18 m ( = 1, . . ., 25).One possible set of parameters for the Gaussian cloud model is  = 10 9 m and  ,peak = 0.1 cm −3 .Further parameters are given in the caption of Fig. 2. We used the extended precision mathematical functions available in Mathematica, Wolfram Research, Inc (2024), and the GCC Quad-Precision Math Library, Free Software Foundation (2024) as we needed to determine trigonometric functions of large arguments.

Secondary spectra
The secondary spectrum is given by the two-dimensional Fourier transform of the preceding expression We note that the secondary spectrum is conventionally defined as the square of the latter quantity.Our usage here coincides with the "conjugate spectrum" used by other authors (e.g.Simard et al. (2019)), although we are considering only the magnitude of that quantity.Often this quantity is also shown on a logarithmic scale (see Fig. F1), in that case H (  ,   ) and H (  ,   ) 2 differ only by a factor 2. The secondary spectrum is shown in Fig. 2. It is obtained by the discrete Fourier transform of (512, 512) points of the dynamic spectra which are zero padded to size (1536, 1536).The secondary spectra show a wealth of sharply delineated features, caused by the presence of the ISM and the specific Fourier integration domain.All results shown in the figures result from a numerical evaluation of Eq. ( 10) and Eq. ( 14).For the interpretation of the numerical results, we discuss different levels of approximations of the integrals in the following sections.

Main parabolic arc features
The Fourier transform of the first term in Eq. ( 10) describes interference between the direct path of the electric field from the pulsar to the observer and the path going through the ISM at a cloud centred at position (0,   ,   ).The Fourier integral comprises terms in the form To evaluate the integrals in the last equation, we set r  () = (  ,   , 0), and r  () = (  , v , , v , ).The expressions derived in this section include the possibility of arbitrary movements of pulsar, observer, and ISM.The argument of the exponential functions is expanded around   = −∞ and   = ∞ to first order.In addition, we consider only the first order in  and  around zero and   , respectively.The denominator is taken to be constant.Using the relation we obtain for the first exponential function The second exponential function in Eq. ( 15  For   = 0 and v , = 0 the last expression reduces to the parabolic arc expression crossing the origin of the   ,   coordinate system.Eq. ( 19) shows that   induces a shift of the parabolic structures in the   -direction and v , a shift in   -direction, moving the parabola away from the the origin of the   ,   coordinate system.The corresponding expressions within the Fresnel-Kirchhoff approach (Walker et al. ( 2004)) for arbitrary two-dimensional scatterer positions and velocities have been given in (Xu et al. 2018, Fig. A1) and in Shi (2021).The   coordinates of the points allow one to construct a projected spatial distribution of the scatterers along the e  axis from the secondary spectra, or, for points clearly offset from the main parabola, to determine their e  coordinate.Note that Cordes et al. (2006) introduce a 1/(2)-scaled variant of the conjugate quantities: conventionally referred to as  and   in the literature, respectively.A detailed quantitative description of the features requires to move beyond the linearised exponents and to evaluate the integrals in terms of special functions, see Appendix D. For simplicity of presentation, we concentrate on the special case v  = 0.For the main parabolic arc we obtain a trapezoid around the centre point ( Here, the ±, ∓ in  in the secondary spectrum with a magnitude proportional to the inverse distance 1/  .The extension of the Fourier window (and in general shape of the chosen window) changes the secondary spectra by affecting the size of the rectangular areas in the secondary spectra.Each of these trapezoids comes with a complex phase leading to interference effects in the case of overlaps of different trapezoids.The right panel of Fig. 2 shows a close-up of the patterns with solid lines drawn according to Eq. ( 21).

Inverted parabolic arcs
The Fourier transform of the second term in Eq. ( 10) arises from interference between two waves travelling through the ISM at positions (  ,   ) and at (  ,   ).
This expression is evaluated as in the last subsection and yields maxima at points By the same replacements as described after Eq. ( 17) the expressions Eqs.(24,25) can be shown to be identical to the ones in Hill et al. (2003Hill et al. ( , 2005)); Cordes et al. (2006).These points lie on inverted parabolic arcs and illuminate rectangular areas in the secondary spectra with vertices and magnitude given by The right panel of Fig. 2 shows a close-up of the patterns with dashed rectangles drawn according to Eq. ( 26).Thus, our analytic expressions (21,26) are in excellent agreement with our numerics.In the "noodle model of scintillation arcs" proposed by Gwinn a smearing of these structures (Gwinn 2019, Eqs. (62,63)) leads to partly similar expressions as our Eqs.( 17), ( 26).The differences in the results originate from the differing approaches: whereas Gwinn analyses how in the expressions for the spot positions in (17,18) and (??,25) change during the integrations in ( 15) and ( 23), respectively, we directly analyse the analytical results for ( 15) and ( 23).Furthermore, the respective magnitudes (22,27) and consequences with respect to the visibility of main vs. inverted arcs are not discussed in Gwinn (2019).The issue of resolution was described in Walker et al. (2004).The authors obtain sinc-functions which lead to decreasing/increasing spot sizes in   (  )-direction with increasing/decreasing Δ  (Δ  ).Existing techniques to combat this are performing the Fourier transform with respect to wavelength instead of frequency Fallows et al. (2014);Reardon et al. (2020) or with respect to time times frequency instead of time Sprenger et al. (2020).We emphasise that our expressions derived in Appendix D describe both effects, the smearing and the resolution.

Magnitudes of main and inverted parabolic arcs
The ratio of the magnitude of an inverted parabolic arc structure compared to a main parabolic structure is given by the expression Eq. ( 28) implies an increased visibility of the inverted arc structures for higher plasma densities (corresponding to a larger value of ).Some pulsars (see pulsar B1508+55 discussed by Sprenger et al. (2022) show a transient evolution of the secondary spectra at different epochs with less and more pronounced inverted arc structure.For a pulsar located 1313 pc from the ISM and observed at a distance of 788 pc from the ISM, | ≈ 2 gives a value of  = 7.5 × 10 18 m to distribute equal intensities to the direct and inverted arc features.This estimate is in general agreement with the transition from a single main arc ( = 10 18 m) to inverted arcs ( = 5 × 10 18 m) seen in Fig. 3.The Helmholtz equation describes the underlying physics in terms of a local change in the refractive index, and thus cannot distinguish between an increase or decrease in matter relative to an average background density.If we take ionised gas (plasma) as the origin of the refractive index change, we can estimate the electron density of the ISM structures, since the  parameter is directly related to the refractive index change.
For the model in Fig. 3 we chose  = 317 scattering clouds in order to avoid considerable overlap of the corresponding panels in the secondary spectrum.For simplicity, all clouds are assigned the same  value.The parameter  = 1 × 10 18 m could be realised with a Gaussian cloud with  = 5 × 10 8 m and  ,peak = 0.15 cm −3 , which is about ten times the average electron density  0 = 0.015 cm −3 in our galaxy (Ocker et al. (2020)).In the large -regime areas in the secondary spectra overlap and additional interference of the individual structures in the secondary spectrum occurs.Refractive index changes might also be coming from neutral gas clouds.The finite extension of all interference structures leads to further interference between overlapping rectangles and trapezoids as seen in the right panel of Fig. 2. The interference causes smaller scale structures compared to the extension of the rectangular or trapezoidal areas.
Fig. 4 displays the dynamic and secondary spectra of a ISM region with a split-off part (see black arrow in the left panel), causing an offset feature in the secondary spectra (white arrow in the right panel).The split-off part produces shifted inverted arclets by including a clump of scatterers offset in the -direction, but still in the same scattering screen (see Fig. 4, left panel).Similar structures have been observed by Brisken et al. (2010) for pulsar B0834+06, but are attributed there to a lens-like concentration of plasma due to their different movement with wavelengths (see also Simard & Pen (2018)), or to multiple screens (Simard et al. (2019), Zhu et al. (2023)).

CONCLUSION
We show that the Born approximation to Green's function is a suitable method for computing dynamic and secondary spectra of pulsar signals.The theoretical description does not use the Fresnel-Kirchhoff diffraction integral.Our method paves the way for the effective extraction of physical parameters such as the refractive index change and the spatial structure of the ISM from secondary spectra.The main and inverted arcs seen in secondary spectra are obtained without assuming a quasi one-dimensional structure of the ISM.Furthermore, within this approach, we are able to compute and analyse the spectra with high precision numerically, as well as explain them analytically.The method can be generalised in several directions, e.g., it is straightforward to describe configurations containing multiple screens or screens extended in three dimensions.The Green's function method could also be applied to plasma lens structures, as e.g. the Gaussian plasma lens (Clegg et al. (1998)), not considered here.
at the observer at (  , 0, 0) and another path starting at the pulsar, going to the ISM offset from the line-of-sight at (0,  0 , 0) and from there to the observer is given by where the last relation holds in the limit of | 0 | ≫  0 and |  | ≫  0 .Shifting the scatterer along the line-of-sight from (0,  0 , 0) to (Δ,  0 , 0) changes the result in Eq. (B1) to As long as |Δ| ≪ |  | and |Δ| ≪ |  |, the impact of Δ on the interference is quite small, i.e. for the setup in Fig. 3 and a scatterer at a distance of 1 a.u.away from the line of sight, a 170,000 a.u.shift along  gives a 0.3 m change in path difference (radio wavelength  = 0.3 m at 1 GHz).In the case |  | = |  |, there is no first order dependence on Δ.To produce a change of Δ of 1 meter at |  | = 500 pc, Δ can extend up to 13 pc, corresponding to 3 percent of the LOS.
Due to these observations we do not consider explicitly the extension of the ISM in -direction.It would be straightforward to include this effect in our calculation and extend our result to that case.
In contrast to a change of the  position of the scatterer, a change in  direction by Δ leads to Comparing the effects of Δ in Eq. (B2) and of Δ in Eq. ( B3), respectively, we see that under the assumption |  | ≈ |  | the effect of Δ is by a factor  0 /|  | smaller than the effect of Δ.

APPENDIX C: BORN APPROXIMATION FOR A GAUSSIAN DISTRIBUTION OF ELECTRONS
Consider the electron density of the th cloud To evaluate the integral of the first order Born approximation ) we use the propagator representation of the free Green's function, which has a Gaussian kernel: For definiteness we set r  = (0,   , 0), r = (, 0, 0), and r ′ = ( ′ , 0, 0).This allows us to perform the spatial integration r ′′ analytically and the remaining integral reads Expanding the integrand in a power series around  = 0 and integrating term by term yields where we introduced  2 =  2 +  2  ,  ′2 =  ′2 +  2  and used the definition of   (Eq.11).The first term is the direct path from the pulsar to the observer, the second term is identical to the interaction of the pulsar pulse with a point scattering source obtained by contracting the Gaussian cloud, and the third term leads to a direction dependent scattering amplitude.We conclude that the contraction of the Gaussian cloud to a point is a valid approximation if the last term in Eq. (C5) can be neglected; otherwise, it should be included and leads to a diminishing effect of scattering clouds off the line of sight.

APPENDIX D: SADDLE POINT EVALUATION
For deriving the finite extensions of the interference regions in the secondary spectra, it is convenient to introduce the effective perpendicular velocity of the interstellar medium (D1)

Figure 1 .
Figure1.Sketch of the scattering setup, including the pulsar, the ISM, and the observer on Earth.The drawing is not to scale, the extension of the ISM along the vertical axis is about 10 9 times exaggerated.We take the ISM to be at rest, while the pulsar and the observer are possibly moving in orthogonal directions with respect to each other and the connecting line observer -pulsar.The coordinate origin is taken to be at the intersection of the ISM and the line of sight.The cones mark the scattering disk.Within the Born approximation structures within the ISM are contracted to point scatterers (black dots).

Figure 2 .
Figure 2. Theoretical dynamic spectra  ( , ), upper left panel, and secondary spectra (| H (  = 2    ,   = 2    ) |, upper right panel.The lower row panels show the secondary spectra in detail (left panel: numerical evaluation of Eq. (14), right panel: analytic result from Eqs. (D3), (D9).The polygons indicate the analytical boundaries of the main arc features (solid lines) and inverted arcs (dashed lines).The centre frequency is set at   = 0.1 GHz to produce large features in the secondary spectra.All scatterers are located on a line perpendicular to the line of sight, intersecting the line of sight.( = 4.5 × 10 18 m,   = −214 pc,   = 429 pc,   = 640 km/s).
) leads to the expressions in Eqs.(17,18)  with the replacement  → −.These expressions agree with the ones inHill et al. (2003Hill et al. ( , 2005) ) andCordes et al. (2006), which can be shown by introducing for    and v ⊥ the corresponding components.

Figure 4 .
Figure 4. Left panel: example of an extended scatterer set with one separated scatterer region (black arrow), leading to the formation of an inverted arc offset from the main parabolic arc in the secondary spectrum | H (  = 2    ,   = 2    ) | (right panel), marked by the white arrow.The dynamic spectrum  ( , ) is shown in the centre panel.Conjugate time   = 2    , conjugate frequency   = 2    .Parameters:   = 0.3 GHz,  = 10 18 m (for all scatterers),   = −214 pc,   = 429 pc,   = 160 km/s.

Figure F1 .
Figure F1.Left panels: Secondary spectra | H (  = 2    ,   = 2    ) | shown in Figs.3,4 on a linear scale, right panels: logarithmic scale of the same data after subtracting the mean of each row.