# The Thousand-Pulsar-Array programme on MeerKAT - X.
Scintillation arcs of 107 pulsars

R. A. Main1, A. Parthasarathy1, S. Johnston2, A. Karastergiou3,4, A. Basu5, A.
D. Cameron6,7, M. J. Keith5, L. S. Oswald3,8, B. Posselt3,9, D. J. Reardon6,7,
X. Song5, P. Weltevrede5
1 Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn,
Germany
2 CSIRO Astronomy and Space Science, Australia Telescope National Facility, PO
Box 76, Epping NSW 1710, Australia
3 Department of Astrophysics, University of Oxford, Denys Wilkinson Building,
Keble Road, Oxford OX1 3RH, UK
4 Department of Physics and Electronics, Rhodes University, PO Box 94,
Grahamstown 6140, South Africa
5 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Manchester M13 9PL, UK
6 Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, P.O. Box 218, Hawthorn, Victoria 3122, Australia
7 Australian Research Council Centre of Excellence for Gravitational Wave
Discovery (OzGrav)
8 Magdalen College, University of Oxford, Oxford OX1 4AU, UK
9 Department of Astronomy & Astrophysics, Pennsylvania State University, 525
Davey Lab, 16802 University Park, PA, USA Email<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present the detection of 107 pulsars with interstellar scintillation arcs
at 856–1712 MHz, observed with the MeerKAT Thousand Pulsar Array Programme.
Scintillation arcs appear to be ubiquitous in clean, high S/N observations,
their detection mainly limited by short observing durations and coarse
frequency channel resolution. This led the survey to be sensitive to nearby,
lightly scattered pulsars with high effective velocity – from a large proper
motion, a screen nearby the pulsar, or a screen near the Earth. We measure the
arc curvatures in all of our sources, which can be used to give an estimate of
screen distances in pulsars with known proper motion, or an estimate of the
proper motion. The short scintillation timescale in J1731$-$4744 implies a
scattering screen within 12 pc of the source, strongly suggesting the
association between this pulsar and the supernova remnant RCW 114. We measure
multiple parabolic arcs of 5 pulsars, all of which are weakly scintillating
with high proper motion. Additionally, several sources show hints of inverted
arclets suggesting scattering from anisotropic screens. Building on this work,
further targeted MeerKAT observations of many of these pulsars will improve
understanding of our local scattering environment and the origins of
scintillation; annual scintillation curves would lead to robust screen
distance measurements, and the evolution of arclets in time and frequency can
constrain models of scintillation.

###### keywords:

pulsars: general – ISM: general

††pubyear: 2022††pagerange: The Thousand-Pulsar-Array programme on MeerKAT -
X.  Scintillation arcs of 107 pulsars–A2

##  1 Introduction

Observations of the radio emission of pulsars exhibits scintillation, an
interference pattern owing to electron density fluctuations in the ionized
interstellar medium (IISM). Recent studies of several pulsars have revealed
“scintillation arcs” - a parabolic distribution of power in the 2D power
spectrum of scintillation commonly referred to as the secondary spectrum.
Discovered by Stinebring et al. (2001), scintillation arcs have proven a
powerful probe of the distribution and geometry of scattering structures in
the IISM, as well as a probe of pulsar emission regions (e.g. Pen et al.
2014), and binary orbits (e.g. Reardon et al. 2020).

In addition to a main parabola, many sources show ‘arclets’, inverted
parabolae of the same curvature stemming from apices on the main parabola,
which persist over timescales of months. This is difficult to explain with a
picture of turbulent, isotropic, volume filling IISM, instead suggesting
compact scattering structures along highly anisotropic thin ‘screens’ (Walker
et al., 2004; Cordes et al., 2006). Some scattering screens have been seen to
persist over many years (Reardon et al., 2020; McKee et al., 2022; Mall et
al., 2022), suggesting the extent of these screens to be $\gtrsim 1000$ AU.
Moreover, arclets are seen to track between observations and across frequency
(Hill et al., 2005; Sprenger et al., 2022), consistent with compact scattering
structures at fixed angular positions.

It is currently unknown precisely how and where these screens originate,
although screen associations have placed screens in supernova remnants (Cordes
et al., 2004; Yao et al., 2021), HII regions (Mall et al., 2022), potentially
near the edge of the local bubble (Reardon et al., 2020; Sprenger et al.,
2022), or in ionized gas surrounding hot stars (Walker et al., 2017). The
electron densities required to scatter radio light at such large angles are
difficult to reconcile with pressure balance in the IISM, leading to models of
scattering in thin over- or under-dense plasma sheets seen edge on (Pen &
Levin, 2014), or in magnetic filaments (Gwinn, 2019).

While most studies of scintillation arcs to date have focused on specific
sources of interest, wider surveys are useful in determining the behaviour and
occurrence of scintillation arcs, as well as the distribution of screens in
the Milky Way. Recent surveys suggest scintillation arcs are common, if not
prevalent. Stinebring et al. (2022) found scintillation arcs in 19 out of 22
detected pulsars in a survey of bright, low Dispersion Measure (DM) pulsars
using the Green Bank Telescope and the Arecibo Observatory. Wu et al. (2022)
detect scintillation arcs in 9 out of 31 sources with LOFAR at $110-190\,$MHz;
scintillation becomes finer in both time and frequency, requiring very fine
channelization and time sampling, resulting in harsh signal-to-noise
limitations in detecting arcs.

The Thousand Pulsar Array (TPA), is part of the Large Survey Project
‘MeerTime’ on the MeerKAT telescope (Bailes et al., 2020), which has observed
$1270$ pulsars in the southern sky in a wide band from 856$-$1712 MHz (details
of the TPA in Johnston et al. 2020). With the wide band, and the massive
upgrade in sensitivity of MeerKAT compared to Parkes, the TPA can be used as a
representative survey of scintillation properties of pulsars in the southern
sky.

In this paper, we present the first results of scintillation using the TPA,
focusing on detection and measurements of scintillation arcs of 107 pulsars.
Section 2 revisits the necessary theory of scintillation arcs, Section 3
describes the data, the reduction, the sample and curvature measurements,
Section 4 describes our results, and Section 5 includes ramifications of our
results and prospects for future work.

##  2 Background

In this section, we briefly review the relevant theory of scintillation arcs,
and refer the reader to Walker et al. (2004); Cordes et al. (2006) for more
detailed treatment. Additionally, in section 2.2, we outline limiting cases,
and the approximations necessary to estimate the screen distance or pulsar
proper motion from a single measure of a scintillation arc curvature.

###  2.1 Scintillation Arcs

The dynamic spectrum $I(\nu,t)=|E(\nu,t)|^{2}$ expresses the pulse intensity
as a function of time and frequency. The 2D power spectrum, or secondary
spectrum $|I(\tau,f_{D})|^{2}$ expresses the intensity in terms of its
conjugate variables $f_{D}$ and $\tau$ which are related to the Doppler shift
and phase delay. In the limit of strong scattering, the secondary resulting
interference pattern can be described through the stationary phase
approximation as a discrete set of deflected images. The secondary spectrum
then contains the interference of all pairs of images

$\displaystyle f_{D,ij}$
$\displaystyle=\frac{(\boldsymbol{\theta_{i}}-\boldsymbol{\theta_{j}})\cdot\boldsymbol{v}_{\rm
eff}}{\lambda},$ (1) $\displaystyle\tau_{ij}$ $\displaystyle=\frac{d_{\rm
eff}(\theta_{i}^{2}-\theta_{j}^{2})}{2c},$ (2)

where $\lambda$ is the observing wavelength, and the effective distance
$d_{\mathrm{eff}}$ and effective velocity $\boldsymbol{v}_{\mathrm{eff}}$ are
defined as

$\displaystyle d_{\rm eff}$ $\displaystyle=(1/s-1)d_{p},$ (3)
$\displaystyle\boldsymbol{v}_{\rm eff}$
$\displaystyle=(1/s-1)\boldsymbol{v}_{p}+\boldsymbol{v}_{\oplus}-\boldsymbol{v}_{l}/s,$
(4) $\displaystyle s$ $\displaystyle\equiv 1-d_{l}/d_{p},$ (5)

where $d_{\mathrm{p}}$, $d_{\mathrm{l}}$ and $\boldsymbol{v}_{\rm p}$,
$\boldsymbol{v}_{\mathrm{l}}$ are the distances and velocities of the pulsar
and the scattering screen, respectively, and $\boldsymbol{v}_{\oplus}$ is the
Earth’s velocity. In weak scattering, the main pulsar image is much brighter
than subimages, resulting in a sharp forward parabola. With the main pulsar
image at $\boldsymbol{\theta_{j}=0}$, the common dependence of $\tau$ and
$f_{D}$ on $\theta$ results in a relation of

$\tau=\eta f_{D}^{2},\quad\eta=\frac{\lambda^{2}d_{\rm eff}}{2cv_{\rm
eff,||}^{2}}.$ (6)

In strong scattering from a highly anisotropic screen, subimages are of
comparable brightness, and the cross-terms result in inverted parabolae
(‘arclets’) of the same curvature with apices located on the main arc.

Following Mall et al. (2022) and Sprenger et al. (2022), we choose not to work
directly with $\eta$, but rather with

$W\equiv\frac{|v_{\rm eff,||}|}{\sqrt{d_{\rm
eff}}}=\frac{\lambda}{\sqrt{2c\eta}},$ (7)

separating the measurements and constants from physical quantities of
interest, where $W$ has the added benefit of approaching $0$ rather than
diverging when $v_{\rm eff,||}\rightarrow 0$.

###  2.2 Approximations

For an isolated pulsar, $W$ contains the degenerate combination of 6 unknown
parameters; the pulsar’s distance $d_{p}$ and 2D velocity
$\boldsymbol{v_{p}}$, the screen’s distance $d_{l}$, and either the parallel
velocity $v_{l,||}$ and axis of anisotropy $\psi$ for a 1D screen, or the 2D
velocity $\boldsymbol{v_{l}}$ for an isotropic 2D screen. Even in pulsars
where the proper motion and distance are independently well measured, a single
measure of $W$ is insufficient to measure the screen properties, where this
degeneracy is typically either broken by using annual variations of $W$ or
$\eta$ (Reardon et al., 2020; Mall et al., 2022; McKee et al., 2022; Sprenger
et al., 2022), or through VLBI or multi-telescope dynamic spectra (Brisken et
al., 2010; Simard et al., 2019). However, one does not typically have annual
scintillation curves, or scintillation cross-spectra; under certain limits
reviewed here, one can obtain meaningful information from $W$.

Screen is near to the pulsar or the Earth –  A value $W$ much greater than
$v_{p}$ is likely to imply local scattering, either near to the pulsar or the
Earth. When scattering is in the pulsar’s local environment (e.g. supernova
remnants Cordes et al. 2004; Yao et al. 2021, or eclipsing binaries, Main et
al. (2018); Lin et al. (2021) ), then $d_{l}\rightarrow d_{p}$, and
$s\rightarrow 0$. The Earth’s velocity is negligible, and the only relevant
terms are the distance and velocity between the source and the screen.
Assuming the screen velocity is negligible, and the screen is isotropic or
aligned with $\boldsymbol{v}_{\rm eff}$,

$d_{pl}\approx\frac{v_{p}^{2}}{W^{2}},$ (8)

where $d_{pl}\equiv d_{p}-d_{l}$.

Interstellar screens $\lesssim 5\,$pc of the Earth have been observed from
scintillation arcs of B1133+16 (McKee et al., 2022), and scintillation of
active galactic nuclei (Wang et al., 2021). Under the same assumptions as
above, for a local screen,

$d_{l}\approx\frac{v_{\earth}^{2}}{W^{2}}.$ (9)

We caution however that screen velocities or misalignment can greatly
influence the estimates of $d_{l}$ or $d_{pl}$, as their contributions to
either are squared.

Screen is halfway to the pulsar -  when the pulsar proper motion is unknown,
it can be estimated given an assumption of the screen distance. If e.g. the
screen is halfway to the pulsar, then $d_{\rm eff}=d_{p}$, $s=1/2$, and

$\mu\approx\frac{W}{\sqrt{d_{p}}}.$ (10)

While highly uncertain, this can result in a rough estimate of pulsar proper
motions. This argument is equivalent to determining proper motions using
scintillation velocities, e.g. Cordes & Rickett (1998).

##  3 Data

The data and processing for this paper are identical to those used in the
census of Posselt et al. (2022); all of the data were taken under the Thousand
Pulsar Array project (Johnston et al., 2020) as part of MeerTime (Bailes et
al., 2020). In brief the observations were taken with the MeerKAT L-band
receiver (centered at $\approx$ 1283.6 MHz), and are folded with 8 s
integrations and $\approx 0.836\,$MHz channels, corresponding to a total range
of 62.5 mHz in $f_{D}$ and $\approx 0.6\mu$s in $\tau$. For one pulsar (see
Section 4.3) we used the data available at a high time resolution of 38.28
$\mu$s to extend the range of $f_{D}$. As observation duration is one of the
limiting factors to detecting scintillation arcs, we limited this study to the
single longest observation of every pulsar, typically $\sim 5-20$ min,
although in a few cases much longer.

###  3.1 Creating Dynamic Spectra

To create the raw dynamic spectra, radio frequency interference (RFI) is
flagged and masked by the standard deviation in the off-pulse region, using
methods identical to Main et al. (2020); Mall et al. (2022). After RFI
masking, the off-pulse is subtracted in every time and frequency
subintegration. The background-subtracted folded spectrum is then multiplied
by the frequency averaged pulse profile, and summed over phase and
polarization to form $I(t,\nu)$.

Due to the high sensitivity and wide band of MeerKAT, the dominant source of
noise in the dynamic spectrum of many pulsars is not receiver noise, but
rather the pulse-to-pulse variations and the window function caused by the RFI
mask. The measured dynamic spectrum is $I_{\rm meas}(t,\nu)=I_{\rm
s}(t,\nu)M_{\rm RFI}(t,\nu)I_{\rm int}(t,\nu)$. The intrinsic variations are
often removed by dividing the frequency-averaged intensity from every
subintegration, but this has the effect of up-weighting noisy integrations
(and vice-versa), and diverges when the pulsar nulls. We try to solve this
problem using a Wiener Filter, following the same steps as Lin et al. (2021).
The one addition is that we multiply our window function by the frequency-
averaged flux in each time bin, $I_{\rm int}(t,\nu)$, in an attempt to filter
out the intrinsic pulse variations. For computational reasons, the full band
is split into 8 sub-bands, and the time axis is split in longer observations
(> 24 minutes, 180 time bins), and the Wiener Filter computed separately in
each. This leads to slight artefacts at the boundaries, resulting in low
$f_{D}$ noise in the secondary spectrum. Characteristic examples of filtered
dynamic spectra are shown in Figure 1.

Figure 1: Characteristic examples of dynamic spectra at 1.4 GHz (bottom) and
corresponding secondary spectra (top). Dynamic spectra are normalized such
that their average intensity is 1, and secondary spectra are plotted in log-
scaling, spanning the full range in $f_{D}$ and $\tau$. J1136+1551 and
J0922+0638 both show arcs; J1136+1551 is weakly scintillating with the
intensity slightly modulated around 1, while J0922+0638 strongly scintillates
with intensity extending to many times the mean. J2215+1538 shows strong
scintillation resolved in frequency, but unresolved in time due to low $v_{\rm
eff}/\sqrt{d_{\rm eff}}$ and short duration. J0835$-$4510 shows strong
scintillation resolved in time but far unresolved in frequency; due to high
S/N, scintillation is still apparent but with greatly reduced modulation.
J1818$-$1422 is scattered to $\approx 12$ ms at L-band (Oswald et al., 2021),
corresponding to $\nu_{s}\approx 10\,$Hz, entirely unresolved by our channels.

###  3.2 Creating Secondary Spectra

Due to the $\lambda^{2}$ scaling of the arc curvature $\eta$, a regular 2D FFT
over a wide band will result in a smearing of scintillation arcs. We account
for this by using the ‘NuT’ transform applied in Sprenger et al. (2021), where
the Fourier Transform in time is performed on a scaled axis
$t^{{}^{\prime}}=\frac{\nu}{\nu_{\rm ref}}t$, which fixes arcs to constant
curvature. In addition, compact features of fixed angular positions on the sky
will be seen at constant positions across frequency in the secondary spectrum.
We choose $\nu_{\rm ref}=1400\,$MHz for all plots in this paper. Before taking
the Fourier transform, the mean of the dynamic spectrum was subtracted, and a
$10\%$ border in both frequency and time were tapered with a Hanning window.

###  3.3 Arc Detections

For every pulsar, diagnostic plots were created which included the dynamic and
secondary spectra, and 2D autocorrelation functions (ACFs) across the full
band, as well as zoom ins to the two clean bands of 920-1020 MHz, and
1300-1500 MHz. Any sources with observable modulation were deemed sources with
either fully or partially resolved scintillation - for each of these pulsars,
secondary spectra were produced with a range of $f_{D}$ limits in both clean
bands. These plots were inspected manually; scintillation arcs were identified
as sources with significant cross-power in the secondary spectrum; ie. power
away from the $f_{D}=0$, $\tau=0$ axes. Characteristic examples showing arc
detections, as well as sources with partially resolved scintillation (in
either time or frequency) but without an arc, are shown in Figure 1.

###  3.4 Arc Curvature Fitting

To fit for $W$, we used the ‘Normalized Secondary Spectra’, which remaps the
$f_{D}$ axis to $f_{D,\rm norm}=f_{D}\sqrt{\tau_{\rm ref}/\tau}$ (Reardon et
al., 2020). This transformation maps parabolae to vertical lines of $f_{D,\rm
norm}\propto W$, and in this paper we take $\tau_{\rm ref}=0.6\mu$s, the
maximum value of $\tau$ allowed by the channel widths. Summing over the $\tau$
axis results in the power as a function of arc curvature $S(W)$, where arcs
can be identified as peaks. We fit peaks in $S(W)$ with an inverted parabola,
where the starting values and fitting range of $W$ are chosen to vary by
source, using an interactive fitter. Arc asymmetries may originate from local
phase gradients across the scattering screen (Cordes et al., 2006; Rickett et
al., 2014), so we fit $W$ separately for positive and negative values of
$f_{D}$; final values of $W$ are the mean from both sides, while we also
record the difference $\Delta W_{lr}$. The error of $W$ is computed from the
half width at half maximum (HWHM) divided by the S/N of the peak, where the N
is the standard deviation of $S(|W|>2W_{\rm fit})$. For sources with screen
distances discussed in the text (all of which are cases of large $W$
suggesting local scattering), the two screen solutions are estimated from
equations 8 and 9. The screen solution near the pulsar uses previously
measured proper motions, while the solution near Earth uses scintools
(Reardon, 2020) to predict Earth’s velocity vector in the plane of the sky in
the direction of the pulsar.

##  4 Results

Of the 1270 pulsars observed, we found 107 to have observable scintillation
arcs. The measured values of $W$, $\Delta W_{lr}$ and derived estimates of
$V_{\rm ISS}$ are in Tables LABEL:table:arcmeas and LABEL:table:arcmeasnov,
and secondary spectra plots are shown in Figure A1.

The detection of arcs is greatly limited by the short durations of
observations and the coarse channels; many sources showed a clear
scintillation pattern in frequency, but with too short a duration to resolve
the pattern in time, and many more sources have scintles far smaller than the
channel bandwidth, and are averaged out in each channel. The channel
resolution means we are biased towards detecting arcs in mildly scattered,
primarily nearby / low DM sources. The short observing durations mean that
only high $W$ sources will show an arc; ie. sources with large velocities, or
screen near to the pulsar or Earth. This is seen in our results $-$ the
inferred proper motion velocities are $\gtrsim 100\,$km/s in most sources, in
many cases much larger.

In addition to the constraints from channel bandwidth and duration, the
detection of arcs can be limited by S/N, or by unmasked RFI, pulsar moding or
nulling, which lead to convolved noise in the secondary spectrum. While we
tried to account for the pulsar’s variable emission in the dynamic spectra
using a Wiener Filter, this is highly suboptimal when pulsars are nulling for
$\gtrsim 50\%$ of the time, and the RFI flagging and excision may not be
perfect in all observations. We plot the DM and S/N of our arc detections
compared to the full sample in Fig. 2. The sources with arcs are primarily
with DM $\lesssim 100\,$pc cm-3, and with S/N $\gtrsim 100$, where higher DM
sources are typically more scattered, with scintillation unresolved due to our
coarse channels.

Figure 2: Sample DM, S/N, and the detections of arcs. Top: cumulative
distribution of sources vs. DM. Bottom: Scatter plot of source S/N vs. DM in
the longest TPA observation for each source, which was inspected for arcs.
Sources with detected arcs were found with DM$\lesssim 100\,$pc cm-3,
S/N$>100$.

###  4.1 Multiple Screens, and Evolution with Frequency

We detect 5 sources with clear evidence of multiple arcs, shown in Figure 3.
These sources all tend to be very low DM, nearby sources, and exhibit weak
scintillation (ie. with $\sigma_{I}/\langle I\rangle\ll 1$, see J1136+1551 in
Fig 1). The fact that the nearest pulsars show multiple screens suggests that
there are many screens along any given line of sight. However, strongly
scintillating sources tend to show a single, dominant arc (e.g. J0922+0638 in
Fig. 1). This can arise either from physical reasons or from selection - it is
possible that there tends to be one dominant screen along the line-of-sight,
perhaps due to having the highest $\Delta$DM, and being near to the halfway
point which maximises the path length difference (and thus $\tau$) for a given
deflection angle. In addition, 2-screen effects could quench the scintillation
from additional screens once the source appears extended, while in weakly
scattered sources the pulsar line-of-sight image dominates and 2-screen
effects are subdominant. A possible selection effect is that sources of strong
scattering result in inverted arclets, which may drown out additional screens
in the secondary spectra. One additional selection effect to consider is that
all of the multi-screen sources have large proper motions of $140-660$ km/s,
where the slowest of these sources (J1057$-$5226) had a comparatively long
87.8 min observation. This suggests that the detection of multiple arcs is
greatly influenced by the resolution in $f_{D}$, in addition to the effects
described above. We discuss future avenues to break these ambiguities in
Section 5.

The outermost arcs in three of our multi-arc sources imply scattering screens
very near the source, with $d_{pl}\approx 35\,$pc, $d_{pl}\approx 27\,$pc,
$d_{pl}\approx 65\,$pc in J0536$-$7543, J1239+2453, and J2048$-$1616
respectively, or within a parsec of the Earth. If the scattering is near
Earth, there would be a strong annual modulation of the arcs, similar to the
outermost arc in B1133+16 (McKee et al., 2022), and the corresponding screens
would have a large angular extent. These multi-screen sources imply a screen
density along the line of sight of $\lesssim 200-300\,$pc/screen, while the 6
screens in B1133+16 shown by McKee et al. (2022) imply $\lesssim
60\,$pc/screen.

Figure 3: Secondary spectra of multi-arc systems, and a selection of sources
showing interesting evolution across frequency, discussed in sections 4.1 and
4.2.

###  4.2 Secondary Spectrum Evolution with Frequency

In Figure 3, in addition to multi-screen sources, we included sources which
show interesting behaviour across frequency. This includes sources with
discrete features at the same location across frequency, consistent with
scattering from compact regions of fixed angular position. PSR J1614+0737
shows an evolution from weak to strong scattering across the band, showing
evidence of inverted arclets at 1 GHz indicative of a highly anisotropic
screen. PSR J1740$-$1311 looks like a diffuse arc at $1\,$GHz, but appears to
show structured power at $1.4\,$GHz. A quantitative comparison of the feature
positions across frequency (e.g. Brisken et al. 2010; Rickett et al. 2021),
and with time (e.g. Hill et al. 2005; Sprenger et al. 2022) will concretely
determine if scattering occurs from compact regions, and may distinguish
between scattering from over- or under-dense regions.

###  4.3 PSR J1731–4744

PSR J1731$-$4744 (B1727$-$47) was recently determined to be a high proper-
motion pulsar with $\mu=151\pm 19$ mas/year, and tracing back the pulsar
motion suggests an association with supernova remnant RCW 114 (Shternin et
al., 2019). In our sample, J1731$-$4744 stood out as having the highest value
of $W$, or equivalently showing an extremely fast $\sim$several second
scintillation timescale with resolvable scintles in frequency, shown in Fig.
4. For this source, we re-reduced the filterbank data with single pulse time
integrations ($\approx 0.83$ s), to fully resolve the arc in $f_{D}$. The
measured value of $W=3830\pm 440$ km s-1 kpc-0.5indicates either a screen near
to the source, with estimates of $s\approx 0.02$, and $d_{pl}\approx 12$ pc,
or a screen at $d_{l}\approx 0.03\,$pc from the Earth. The SNR shell of RCW
114 has a radius of $\approx 2^{\circ}$, $\approx 20\,$pc at $d_{\rm p}\approx
600$ pc, comparable to $d_{\rm pl}$. Our result strongly suggests the
association of J1731$-$4744 with RCW 114, in agreement with Shternin et al.
(2019).

Figure 4: Dynamic and secondary spectrum of J1731$-$4744, rereduced with
single pulses as described in section 4.3. The large effective velocity
implies a screen distance of $d_{\rm ls}\sim 12\,$pc, comparable to the radius
of SNR RCW114 which the pulsar is coincident with.

###  4.4 Comparison with Literature

Several of the sources in our sample have had previously observed
scintillation arcs. We highlight and compare a few of these results.

#### J0837+0610/B0834+06

is one of the best studied pulsars for scintillation, showing dramatic
inverted arclets, being used for phase retrieval (Walker et al., 2008), and
imaged through VLBI to be a highly anisotropic screen (Brisken et al., 2010)
at a distance of $d_{l}\approx 420$ pc. In our data, we detect only a forward
arc without inverted arclets. Our estimate of $d_{l}\approx 412\,$pc is quite
close to the more robust value from Brisken et al. (2010), suggesting we are
still seeing the same scattering screen, where this agreement is likely aided
by the fact that the screen is (seemingly coincidentally) closely aligned with
the pulsar proper motion. The transition from strong to weak scattering is
consistent with the findings of Smirnova et al. (2020), who observed the same
change in the source’s scintillation between 2012 and 2015.

#### J1057-5226/B1055-52

was studied with Parkes data at 1400 MHz by Kerr et al. (2018), finding a
prolonged period of increased scattering (ie. decreased scintillation
bandwidth) from MJD 54725–55725, interpreted as an extreme scattering event
(ESE). During the ESE, the secondary spectra were dominated by a single arc
with inverted arclets at curvature of $\eta_{\rm outer}\approx 0.02$s3. In
observations at MJD 57506–57509, long past the ESE, a second, interior arc is
visible at $\eta_{\rm inner}\approx 0.06$s3. These arcs correspond to $W_{\rm
outer}\approx 340$ km s-1 kpc-0.5, $W_{\rm inner}\approx 200$ km s-1 kpc-0.5,
in close agreement with the two arcs seen in our data at $W_{\rm outer}=372\pm
79$ km s-1 kpc-0.5, $W_{\rm inner}=166\pm 14$ km s-1 kpc-0.5on MJD 58655. This
suggests that the nature of scattering is similar to the state of ‘quiescent’
scattering observed in the later observations of Kerr et al. (2018); the same
arcs being visible over many years, combined with the pulsar’s transverse
velocity of $\approx 140$ km/s implies the transverse extent of the two
screens is $\gtrsim 100$ AU, $\gtrsim 300$ AU respectively.

#### J1136+1551/B1133+16

was studied by McKee et al. (2022), who measure the variations of 6 arcs
across 40 years. In our single observation, we clearly detect three of these
arcs, which correspond to their arcs B, C, and E. In longer observations with
finer frequency resolution, we could likely further disentangle the multiple
screens.

#### J1239+2453/B1237+25

shows clear evidence of 3 screens in our data, across the full band. In
previous LOFAR observations at 110-190 MHz, the secondary spectra show only a
diffuse blend of power (Wu et al., 2022). Additionally, only the inner arc was
faintly seen by Fadeev et al. (2018), who observed the source with GBT at
$316-332$ MHz.

#### J1921+2153/B1919+21

showed a clear arc at LOFAR in Wu et al. (2022). Their values of $\eta\approx
4.4\,s^{3}$ at observing frequency of $150\,$MHz correspond to $W\approx 216$
km s-1 kpc-0.5, in agreement with our measured value of $W=243\pm 25$ km s-1
kpc-0.5, suggesting we are observing the same screen.

#### J0452$-$1759/B0450$-$18

was studied in detail by Rickett et al. (2021) at 340 and 800 MHz with the
GBT, showing clear inverted arclets and being fit well with a 1D screen
brightness distribution. In our data, we miss this structure, seeing only the
outline of an arc. This shows the limitations of our short observing time (9
min, compared to 60 min observations), and that rich structures can be
revealed using longer integrations. It is also possible that our observation
is at a low point of $v_{\rm eff,||}$ from annual motion.

##  5 Ramifications

In summary, we have presented the discovery of 107 sources with detectable
scintillation arcs with the Meerkat TPA. In this section we discuss
ramifications and future avenues.

Our findings are consistent with recent suggestions that scintillation arcs
are prevalent (Stinebring et al., 2022), their detection primarily limited by
frequency resolution and observation duration, provided the dynamic spectrum
is high S/N and largely devoid of RFI. The nature of this survey manifests in
strong biases, where arcs are easiest to detect in weakly scattered sources
with high effective velocity - either through a high proper motion, or local
screens, as in the case of J1731$-$4744\. The detection of multi-arc sources
include even more selection biases, with multi-arc sources primarily being
nearby, weakly scintillating sources with very high transverse velocity. There
are likely many screens on any given sightline which cannot be seen in
secondary spectra due to insufficient time duration, 2-screen interactions, or
inverted arcs.

This survey leaves many avenues for further research. In the TPA alone, there
are hundreds more sources with resolvable scintillation, but without the
detection of arcs. Analysis of the full sample can provide measures of
scattering time $\tau$ on many more sightlines, which will help improve
galactic electron models (e.g. Cordes & Lazio 2002; Yao et al. 2017). The more
distant, scattered sources will be best studied at higher frequencies where
arcs with be sharper, more easily resolvable. Through wide-band/multi-
frequency observations, or through phase-retrieval techniques such as cyclic
spectroscopy (Demorest, 2011; Walker et al., 2013), one could measure the
correspondence between time-domain scattering and scintillation arcs. This may
help elucidate the morphology of scattering tails, and investigate if
seemingly isotropic scattering could be reproduced through the ensemble of
many anisotropic screens (e.g. Oswald et al. 2021).

With only a single arc curvature measurement, only crude estimates of screen
distances can be made. Additional observations to obtain annual curves
(Reardon et al., 2020; Main et al., 2020; McKee et al., 2022; Mall et al.,
2022), or VLBI (Brisken et al., 2010) can fully solve for screen distances and
geometry. In multi-screen sources, VLBI combined with dynamic cross-spectra
can be used to simultaneously solve for multiple screens (Simard et al.,
2019), and characteristic 2-screen effects may also provide additional
constraints on screen geometry, as seen in B1508+55 (Sprenger et al., 2022).
With a growing number of scattering screen measurements, it may be possible to
map the distribution of screens in the Milky Way, cross-matching with
foreground stars, HII regions, supernova remnants, and with the local bubble
and other local plasma. With high S/N, one could probe magnetism in scattering
screens, and compare scintillation across pulse phase to resolve pulsar
emission. The results of this survey will then serve as a valuable stepping
stone to studying the distribution of free electrons in the Milky Way, probing
the dynamics and emission of pulsars, and elucidating the astrophysical origin
of interstellar scintillation.

## Acknowledgements

We thank the reviewer William Coles for many useful comments which improved
this manuscript. RAM thanks Tim Sprenger and Olaf Wucknitz for their
discussions and insights, and Fang Xi Lin and Daniel Baker for help with the
Wiener Filter. LSO acknowledges the support of Magdalen College, Oxford. The
MeerKAT telescope is operated by SARAO, which is a facility of the National
Research Foundation, an agency of the Department of Science and Innovation.
MeerTime data are housed and processed on the OzSTAR supercomputer at
Swinburne University of Technology, funded by Swinburne University of
Technology and the National Collaborative Research Infrastructure Strategy
(NCRIS).

## Data Availability

The dynamic spectra used in this paper are included on Zenodo at
https://doi.org/10.5281/zenodo.7261413.

## References

  * Bailes et al. (2020) Bailes M., et al., 2020, Publ. Astron. Soc. Australia, 37, e028
  * Brisken et al. (2002) Brisken W. F., Benson J. M., Goss W. M., Thorsett S. E., 2002, ApJ, 571, 906
  * Brisken et al. (2003) Brisken W. F., Thorsett S. E., Golden A., Goss W. M., 2003, ApJ, 593, L89
  * Brisken et al. (2010) Brisken W. F., Macquart J. P., Gao J. J., Rickett B. J., Coles W. A., Deller A. T., Tingay S. J., West C. J., 2010, ApJ, 708, 232
  * Chatterjee et al. (2001) Chatterjee S., Cordes J. M., Lazio T. J. W., Goss W. M., Fomalont E. B., Benson J. M., 2001, ApJ, 550, 287
  * Chatterjee et al. (2009) Chatterjee S., et al., 2009, ApJ, 698, 250
  * Cordes & Lazio (2002) Cordes J. M., Lazio T. J. W., 2002, arXiv e-prints, pp astro–ph/0207156
  * Cordes & Rickett (1998) Cordes J. M., Rickett B. J., 1998, ApJ, 507, 846
  * Cordes et al. (2004) Cordes J. M., Bhat N. D. R., Hankins T. H., McLaughlin M. A., Kern J., 2004, ApJ, 612, 375
  * Cordes et al. (2006) Cordes J. M., Rickett B. J., Stinebring D. R., Coles W. A., 2006, ApJ, 637, 346
  * Deller et al. (2009) Deller A. T., Tingay S. J., Bailes M., Reynolds J. E., 2009, ApJ, 701, 1243
  * Deller et al. (2019) Deller A. T., et al., 2019, ApJ, 875, 100
  * Demorest (2011) Demorest P. B., 2011, MNRAS, 416, 2821
  * Fadeev et al. (2018) Fadeev E. N., Andrianov A. S., Burgin M. S., Popov M. V., Rudnitskiy A. G., Shishov V. I., Smirnova T. V., Zuga V. A., 2018, MNRAS, 480, 4199
  * Gwinn (2019) Gwinn C. R., 2019, MNRAS, 486, 2809
  * Hill et al. (2005) Hill A. S., Stinebring D. R., Asplund C. T., Berwick D. E., Everett W. B., Hinkel N. R., 2005, ApJ, 619, L171
  * Johnston et al. (2020) Johnston S., et al., 2020, MNRAS, 493, 3608
  * Kerr et al. (2018) Kerr M., Coles W. A., Ward C. A., Johnston S., Tuntsov A. V., Shannon R. M., 2018, MNRAS, 474, 4637
  * Lin et al. (2021) Lin F. X., Main R. A., Verbiest J. P. W., Kramer M., Shaifullah G., 2021, MNRAS, 506, 2824
  * Liu et al. (2016) Liu S., Pen U.-L., Macquart J. P., Brisken W., Deller A., 2016, MNRAS, 458, 1289
  * Main et al. (2018) Main R., et al., 2018, Nature, 557, 522
  * Main et al. (2020) Main R. A., et al., 2020, MNRAS, 499, 1468
  * Mall et al. (2022) Mall G., et al., 2022, MNRAS, 511, 1104
  * McKee et al. (2022) McKee J. W., Zhu H., Stinebring D. R., Cordes J. M., 2022, ApJ, 927, 99
  * Oswald et al. (2021) Oswald L. S., et al., 2021, MNRAS, 504, 1115
  * Pen & Levin (2014) Pen U.-L., Levin Y., 2014, MNRAS, 442, 3338
  * Pen et al. (2014) Pen U. L., Macquart J. P., Deller A. T., Brisken W., 2014, MNRAS, 440, L36
  * Posselt et al. (2022) Posselt B., et al., 2022, MNRAS, submitted
  * Reardon (2020) Reardon D. J., 2020, Scintools: Pulsar scintillation data tools, Astrophysics Source Code Library, record ascl:2011.019 (ascl:2011.019)
  * Reardon et al. (2019) Reardon D. J., Coles W. A., Hobbs G., Ord S., Kerr M., Bailes M., Bhat N. D. R., Venkatraman Krishnan V., 2019, MNRAS, 485, 4389
  * Reardon et al. (2020) Reardon D. J., et al., 2020, ApJ, 904, 104
  * Rickett et al. (2014) Rickett B. J., et al., 2014, ApJ, 787, 161
  * Rickett et al. (2021) Rickett B. J., Stinebring D. R., Zhu H., Minter A. H., 2021, ApJ, 907, 49
  * Shternin et al. (2019) Shternin P., Kirichenko A., Zyuzin D., Yu M., Danilenko A., Voronkov M., Shibanov Y., 2019, ApJ, 877, 78
  * Simard et al. (2019) Simard D., Pen U. L., Marthi V. R., Brisken W., 2019, MNRAS, 488, 4963
  * Smirnova et al. (2020) Smirnova T. V., et al., 2020, MNRAS, 496, 5149
  * Sprenger et al. (2021) Sprenger T., Wucknitz O., Main R., Baker D., Brisken W., 2021, MNRAS, 500, 1114
  * Sprenger et al. (2022) Sprenger T., Main R., Wucknitz O., Mall G., Wu J., 2022, MNRAS, 515, 6198
  * Stinebring et al. (2001) Stinebring D. R., McLaughlin M. A., Cordes J. M., Becker K. M., Goodman J. E. E., Kramer M. A., Sheckard J. L., Smith C. T., 2001, ApJ, 549, L97
  * Stinebring et al. (2022) Stinebring D. R., et al., 2022, arXiv e-prints, p. arXiv:2207.08756
  * Walker et al. (2004) Walker M. A., Melrose D. B., Stinebring D. R., Zhang C. M., 2004, MNRAS, 354, 43
  * Walker et al. (2008) Walker M. A., Koopmans L. V. E., Stinebring D. R., van Straten W., 2008, MNRAS, 388, 1214
  * Walker et al. (2013) Walker M. A., Demorest P. B., van Straten W., 2013, ApJ, 779, 99
  * Walker et al. (2017) Walker M. A., Tuntsov A. V., Bignall H., Reynolds C., Bannister K. W., Johnston S., Stevens J., Ravi V., 2017, ApJ, 843, 15
  * Wang et al. (2021) Wang Y., Tuntsov A., Murphy T., Lenc E., Walker M., Bannister K., Kaplan D. L., Mahony E. K., 2021, MNRAS, 502, 3294
  * Wu et al. (2022) Wu Z., et al., 2022, A&A, 663, A116
  * Yao et al. (2017) Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29
  * Yao et al. (2021) Yao J., et al., 2021, Nature Astronomy, 5, 788

Table 1: Arc measurements in sources with known proper motions. $W$ and
$\Delta W_{lr}$ are computed as described in Section 3.4. Values of dpsr with
errorbars come from parallax measurements cited below, otherwise they are DM
distances from NE2001 (Cordes & Lazio, 2002). In $\Delta W_{lr}$, $-$ or $+$
denote cases where the arc could only be measured on the negative or positive
side of $f_{\rm D,norm}$. Distance references: (1) Deller et al. (2009), (2)
Brisken et al. (2003), (3) Chatterjee et al. (2009), (4) Liu et al. (2016),
(5) Chatterjee et al. (2001), (6) Deller et al. (2019), (7) Brisken et al.
(2002)

Source | DM | $d_{p}$ | Vpm | $W$ | $\Delta W_{lr}$
---|---|---|---|---|---
| pc cm-3 | kpc | km s-1 | km s-1 | km s-1
|  |  |  | kpc-0.5 | kpc-0.5
J0151-0635 | 25.66 | 1.2 | 69 | $88\pm 29$ | $-$
J0206-4028 | 12.9 | 0.6 | 76 | $98\pm 36$ | $-$
J0304+1932 | 15.657 | 0.6 | 110 | $51\pm 3$ | $-4\pm 3$
J0452-1759 | 39.903 | 2.4 | 145 | $134\pm 25$ | $-35\pm 13$
J0536-7543 | 18.58 | 0.8 | 255 | $1488\pm 18$ | $-25\pm 17$
.. | .. | .. | .. | $907\pm 83$ | $-117\pm 7$
.. | .. | .. | .. | $574\pm 28$ | $-$
.. | .. | .. | .. | $191\pm 119$ | $-$
J0630-28341 | 34.425 | 0.333${}^{+0.051}_{-0.039}$ | 81 | $183\pm 38$ | $-54\pm 14$
J0659+14142 | 13.94 | 0.288${}^{+0.033}_{-0.027}$ | 60 | $209\pm 84$ | $-118\pm 33$
J0758-1528 | 63.327 | 3.0 | 219 | $127\pm 38$ | $54\pm 28$
J0820-13503 | 40.938 | 1.961${}^{+0.167}_{-0.109}$ | 418 | $322\pm 57$ | $-27\pm 57$
J0837+06104 | 12.864 | 0.613${}^{+0.062}_{-0.052}$ | 148 | $267\pm 6$ | $-7\pm 6$
J0908-1739 | 15.879 | 0.9 | 208 | $345\pm 47$ | $63\pm 47$
J0922+06385 | 27.299 | 1.205${}^{+0.224}_{-0.163}$ | 505 | $338\pm 64$ | $91\pm 9$
J1057-5226 | 29.69 | 0.7 | 143 | $166\pm 14$ | $-20\pm 6$
.. | .. | .. | .. | $372\pm 79$ | $50\pm 80$
J1136+15516 | 4.841 | 0.372${}^{+0.002}_{-0.002}$ | 660 | $1171\pm 226$ | $-31\pm 229$
.. | .. | .. | .. | $764\pm 121$ | $26\pm 122$
.. | .. | .. | .. | $419\pm 20$ | $-$
J1141-6545 | 116.08 | 2.5 | 40 | $193\pm 36$ | $47\pm 38$
J1239+24537 | 9.252 | 0.862${}^{+0.064}_{-0.056}$ | 472 | $2462\pm 152$ | $-318\pm 152$
.. | .. | .. | .. | $1288\pm 59$ | $66\pm 59$
.. | .. | .. | .. | $342\pm 22$ | $15\pm 6$
J1257-1027 | 29.634 | 1.6 | 103 | $306\pm 35$ | $50\pm 13$
J1430-6623 | 65.102 | 1.0 | 176 | $305\pm 40$ | $-57\pm 24$
J1534-5334 | 24.82 | 1.1 | 160 | $119\pm 67$ | $95\pm 17$
J1645-03176 | 35.756 | 3.968${}^{+0.324}_{-0.397}$ | 386 | $287\pm 63$ | $-74\pm 63$
J1650-1654 | 43.25 | 1.5 | 114 | $146\pm 28$ | $39\pm 15$
J1703-18466 | 49.551 | 2.857${}^{+0.476}_{-0.357}$ | 229 | $267\pm 38$ | $54\pm 12$
J1731-4744 | 123.056 | 0.6 | 429 | $3826\pm 440$ | $-489\pm 443$
J1740+1311 | 48.668 | 1.5 | 208 | $365\pm 15$ | $13\pm 11$
J1825-0935 | 19.383 | 0.9 | 66 | $82\pm 22$ | $-31\pm 16$
J1844+1454 | 41.486 | 2.2 | 476 | $126\pm 59$ | $-$
J1900-2600 | 37.994 | 1.2 | 286 | $415\pm 64$ | $48\pm 64$
J1901-09066 | 72.677 | 1.961${}^{+0.167}_{-0.237}$ | 183 | $189\pm 36$ | $-50\pm 26$
J1921+2153 | 12.444 | 1.1 | 186 | $243\pm 25$ | $-32\pm 25$
J1941-2602 | 50.036 | 1.7 | 123 | $131\pm 34$ | $-47\pm 34$
J1946-2913 | 44.309 | 1.5 | 277 | $76\pm 8$ | $11\pm 8$
J2046-0421 | 35.799 | 1.7 | 96 | $50\pm 44$ | $-63\pm 22$
J2048-16163 | 11.456 | 0.952${}^{+0.018}_{-0.026}$ | 511 | $466\pm 43$ | $-66\pm 21$
.. | .. | .. | .. | $1508\pm 35$ | $-78\pm 17$
.. | .. | .. | .. | $1927\pm 109$ | $-225\pm 109$
J2144-39331 | 3.35 | 0.165${}^{+0.017}_{-0.014}$ | 130 | $319\pm 36$ | $-51\pm 5$
J2317+21496 | 20.87 | 1.961${}^{+0.213}_{-0.206}$ | 79 | $56\pm 14$ | $-$
J2330-2005 | 8.456 | 0.4 | 138 | $143\pm 8$ | $-11\pm 5$
J2346-06096 | 22.504 | 3.636${}^{+0.548}_{-0.258}$ | 733 | $534\pm 23$ | $+$

∗: J1141$-$6545 is dominated by its binary motion, resulting in orbital
variations of $W$ (Reardon et al., 2019).

Table 2: Arc measurements and estimates of $V_{\rm ISS}$ in sources without known proper motions. The errors on $V_{\rm ISS}$ include only the measurement error on $W$, not the large unknown uncertainty in the distance, and are subject to the assumptions and biases discussed in section 2.2. Source | DM | dp | $W$ | $V_{\rm ISS}$ | $\Delta W_{lr}$
---|---|---|---|---|---
| pc cm-3 | kpc | km s-1 | km s-1 | km s-1
|  |  | kpc-0.5 |  | kpc-0.5
J0211-8159 | 24.36 | 1.0 | $224\pm 35$ | $224\pm 35$ | $9\pm 35$
J0343-3000 | 20.2 | 0.9 | $170\pm 8$ | $164\pm 8$ | $-$
J0455-6951 | 94.7 | 43.4 | $23\pm 9$ | $152\pm 59$ | $+$
J0459-0210 | 21.02 | 0.9 | $144\pm 17$ | $138\pm 16$ | $-7\pm 17$
J0529-6652 | 103.31 | 46.2 | $378\pm 91$ | $2569\pm 619$ | $+$
J0543+2329 | 77.703 | 2.1 | $402\pm 129$ | $577\pm 185$ | $-182\pm 48$
J0633-2015 | 90.7 | 6.6 | $27\pm 8$ | $69\pm 21$ | $-$
J0636-4549 | 26.31 | 1.3 | $376\pm 78$ | $431\pm 89$ | $-110\pm 76$
J0656-2228 | 32.39 | 1.9 | $117\pm 14$ | $161\pm 19$ | $9\pm 14$
J0909-7212 | 54.3 | 1.9 | $96\pm 74$ | $132\pm 102$ | $105\pm 12$
J0919-6040 | 82.5 | 2.5 | $34\pm 12$ | $54\pm 19$ | $-$
J0924-5814 | 57.4 | 1.8 | $140\pm 17$ | $190\pm 23$ | $14\pm 17$
J1001-5939 | 113.0 | 2.8 | $457\pm 38$ | $758\pm 63$ | $-54\pm 38$
J1018-1642 | 48.82 | 3.2 | $89\pm 45$ | $158\pm 80$ | $-63\pm 16$
J1020-5921 | 80.0 | 2.1 | $154\pm 28$ | $222\pm 40$ | $21\pm 28$
J1032-5206 | 139.0 | 3.5 | $161\pm 15$ | $299\pm 28$ | $-$
J1055-6905 | 142.8 | 3.9 | $118\pm 51$ | $234\pm 101$ | $31\pm 51$
J1112-6926 | 148.4 | 4.0 | $106\pm 29$ | $212\pm 58$ | $-41\pm 19$
J1119-7936 | 27.4 | 1.0 | $116\pm 17$ | $117\pm 17$ | $-1\pm 17$
J1132-4700 | 123.0 | 5.2 | $262\pm 81$ | $596\pm 184$ | $-114\pm 48$
J1141-3107 | 30.77 | 1.2 | $257\pm 47$ | $283\pm 52$ | $-$
J1210-6550 | 37.0 | 1.1 | $125\pm 38$ | $134\pm 41$ | $-54\pm 5$
J1232-4742 | 26.0 | 0.9 | $129\pm 12$ | $125\pm 12$ | $4\pm 15$
J1239-6832 | 94.3 | 2.1 | $128\pm 11$ | $186\pm 16$ | $-2\pm 11$
J1240-4124 | 44.1 | 1.5 | $150\pm 59$ | $185\pm 73$ | $84\pm 22$
J1312-6400 | 93.0 | 1.9 | $78\pm 8$ | $109\pm 11$ | $+$
J1340-6456 | 76.99 | 1.7 | $150\pm 54$ | $196\pm 70$ | $-$
J1404+1159 | 18.499 | 1.4 | $97\pm 4$ | $115\pm 5$ | $12\pm 3$
J1418-3921 | 60.49 | 1.2 | $27\pm 8$ | $29\pm 9$ | $11\pm 7$
J1440-6344 | 124.2 | 3.1 | $141\pm 28$ | $247\pm 49$ | $-40\pm 20$
J1443-5122 | 87.0 | 1.9 | $94\pm 42$ | $131\pm 59$ | $60\pm 11$
J1514-4834 | 51.5 | 1.3 | $176\pm 30$ | $201\pm 34$ | $-38\pm 30$
J1527-3931 | 48.8 | 1.5 | $120\pm 29$ | $147\pm 36$ | $-$
J1530-5327 | 49.6 | 1.2 | $1064\pm 218$ | $1181\pm 242$ | $308\pm 59$
J1535-4114 | 66.28 | 2.0 | $160\pm 20$ | $224\pm 28$ | $3\pm 21$
J1537-4912 | 69.7 | 2.5 | $107\pm 20$ | $168\pm 31$ | $-$
J1549+2113 | 24.055 | 2.0 | $480\pm 25$ | $673\pm 35$ | $-36\pm 17$
J1603-2712 | 46.201 | 1.6 | $218\pm 45$ | $272\pm 56$ | $-64\pm 25$
J1605-5257 | 34.9 | 2.2 | $203\pm 57$ | $301\pm 84$ | $-80\pm 26$
J1614+0737 | 21.395 | 1.0 | $524\pm 8$ | $515\pm 8$ | $1\pm 8$
J1646-6831 | 43.0 | 1.3 | $23\pm 7$ | $26\pm 8$ | $0\pm 7$
J1720+2150 | 40.719 | 2.3 | $49\pm 15$ | $74\pm 23$ | $-74\pm 8$
J1720-1633 | 44.83 | 1.3 | $337\pm 69$ | $390\pm 80$ | $-98\pm 11$
J1728-0007 | 41.09 | 1.6 | $298\pm 99$ | $374\pm 124$ | $141\pm 67$
J1740+1000 | 23.897 | 1.2 | $456\pm 191$ | $508\pm 213$ | $271\pm 26$
J1745-0129 | 89.3 | 3.1 | $137\pm 60$ | $242\pm 106$ | $85\pm 30$
J1755-0903 | 63.67 | 1.8 | $257\pm 23$ | $344\pm 31$ | $-32\pm 21$
J1800-0125 | 51.0 | 1.8 | $641\pm 37$ | $850\pm 49$ | $40\pm 41$
J1810-5338 | 45.0 | 1.3 | $72\pm 13$ | $82\pm 15$ | $16\pm 13$
J1811-4930 | 44.0 | 1.2 | $118\pm 6$ | $132\pm 7$ | $+$
J1813+1822 | 60.8 | 2.9 | $96\pm 17$ | $164\pm 29$ | $24\pm 13$
J1819+1305 | 64.808 | 2.8 | $107\pm 16$ | $178\pm 27$ | $22\pm 12$
J1825+0004 | 56.618 | 1.9 | $107\pm 32$ | $148\pm 44$ | $+$
J1840-1419 | 19.4 | 0.8 | $26\pm 5$ | $24\pm 5$ | $7\pm 3$
J1849+0409 | 63.97 | 2.6 | $150\pm 12$ | $241\pm 19$ | $+$
J1900-7951 | 39.0 | 1.4 | $179\pm 40$ | $214\pm 48$ | $57\pm 13$
J1917+0834 | 29.18 | 1.9 | $46\pm 24$ | $64\pm 33$ | $+$
J1929+2121 | 66.0 | 3.4 | $42\pm 7$ | $77\pm 13$ | $-$
J1943-1237 | 28.918 | 1.2 | $180\pm 24$ | $198\pm 26$ | $-33\pm 28$
Table 3: Table 2 continued. Source | DM | dp | $W$ | $V_{\rm ISS}$ | $\Delta W_{lr}$
---|---|---|---|---|---
| pc cm-3 | kpc | km s-1 | km s-1 | km s-1
|  |  | kpc-0.5 |  | kpc-0.5
J1945-0040 | 59.71 | 2.5 | $76\pm 8$ | $119\pm 13$ | $0\pm 8$
J1949-2524 | 23.07 | 0.7 | $337\pm 11$ | $281\pm 9$ | $-$
J1951+1123 | 31.29 | 2.1 | $100\pm 5$ | $145\pm 7$ | $1\pm 5$
J2005-0020 | 35.93 | 1.8 | $587\pm 173$ | $781\pm 230$ | $-$
J2038-3816 | 33.96 | 1.4 | $163\pm 43$ | $190\pm 50$ | $-$
J2043+2740 | 21.021 | 1.8 | $39\pm 11$ | $52\pm 15$ | $-15\pm 3$
J2108-3429 | 30.22 | 1.2 | $180\pm 12$ | $198\pm 13$ | $-$
J2139+2242 | 44.16 | 2.8 | $101\pm 14$ | $168\pm 23$ | $-6\pm 14$
J2155-3118 | 14.85 | 0.6 | $55\pm 12$ | $43\pm 9$ | $-$
J2253+1516 | 29.204 | 1.7 | $54\pm 13$ | $70\pm 17$ | $-$

Figure A1: Panorama of Secondary Spectra. The colourbars are logarithmic,
extending from the median to the maximum intensity. The $\tau$ axes are all
identical, covering the full observable range, while the $f_{D}$ range is
chosen separately for each source to best show the arc.
Figure A2: Continued..


# Multidimensional deconvolution with shared bases

Daria Sushnikova, Matteo Ravasi and David Keyes
Extreme Computing Research Center
King Abdullah University of Science and Technology  0138-B1, KAUST, Thuwal,
Saudi Arabia, zip 23955<EMAIL_ADDRESS>
###### Abstract

We address the estimation of seismic wavefields by means of Multidimensional
Deconvolution (MDD) for various redatuming applications. While offering more
accuracy than conventional correlation-based redatuming methods, MDD faces
challenges due to the ill-posed nature of the underlying inverse problem and
the requirement to handle large, dense, complex-valued matrices. These
obstacles have long limited the adoption of MDD in the geophysical community.
Recent interest in this technology has spurred the development of new
strategies to enhance the robustness of the inversion process and reduce its
computational overhead. We present a novel approach that extends the concept
of block low-rank approximations, usually applied to linear operators, to
simultaneously compress the operator, right-hand side, and unknowns. This
technique greatly alleviates the data-heavy nature of MDD. Moreover, since in
3d applications the matrices do not lend themselves to global low rank
approximations, we introduce a novel $\mathcal{H}^{2}$ -like approximation. We
aim to streamline MDD implementations, fostering efficiency and controlling
accuracy in wavefield reconstruction. This innovation holds potential for
broader applications in the geophysical domain, possibly revolutionizing the
analysis of multi-dimensional seismic datasets.

  * April 1, 2024

##  1 Introduction

Seismic recordings provide crucial insights into the interaction of waves with
an otherwise unknown physical medium. A so-called wave state (defined by the
wave nature, medium parameters, source type, and boundary conditions) can be
associated with any such recording [8]. While, in theory, one has the
flexibility to choose any combination of such conditions, real-world seismic
exploration faces constraints stemming from physical, economical, and
environmental factors.

Several data-driven approaches have been developed under the framework of
seismic interferometry (e.g., [33]) to create ideal settings, also known as
virtual surveys, from physically acquired data. By defining interaction
quantities and employing reciprocity theorems [8], these approaches can
leverage wavefield measurements from distinct receivers to emulate responses
as if a source was placed in one of the receivers’ locations. Ideally, cross-
correlating the wavefield recorded by two receivers surrounded by a closed
boundary of sources in a lossless medium yields the Green’s function between
such receivers. This principle has been extensively used in seismic
applications for a variety of tasks [14, 22, 7, 24, 2, 23].

However, closed boundaries are rarely available in real-life applications.
This, combined with the inevitable presence of losses in the propagating
medium, can hinder accurate reconstruction, yielding a blurred Green’s
function from cross-correlation interferometry [35]. An alternative approach
that can accommodate both dissipative media and open boundaries relies on a
convolution-type reciprocity theorem. [35] showed that this leads to an
implicit Green’s function representation; deconvolution can then be employed
to retrieve such Green’s function mitigating the aforementioned correlation-
induced, spatio-temporal blurring effects. This method, known as seismic
interferometry by a Multidimensional Deconvolution (MDD) has been initially
proposed in the acoustic setting and later extended to more general vector
fields, including elastic [30] and electromagnetic waves [34].

While the adoption of MDD constitutes a significant leap towards achieving
more accurate redatuming compared to traditional correlation methods, it comes
with numerical and computational challenges. The former stems from the ill-
posed nature of the inverse problem that one has to solve to estimate the
sought after Green’s function, whilst the latter is related to the arithmetic
and storage complexity of the problem. MDD algorithms must in fact manipulate
dense matrices for the operator, unknowns, and right-hand side, which
significantly escalates the computational demands as the problem size grows,
especially in 3d applications. This has resulted in a somewhat constrained
uptake of MDD within the geophysical community.

However, the superior capabilities of MDD have led the geophysical community
to identify strategies aimed at bolstering inversion stability and alleviating
computational overhead. For instance, to enhance stability, researchers in
[17] apply singular-value decomposition (SVD) to compute the pseudoinverse
matrix and truncate small eigenvalues as a way to naturally regularize the
inverse process. More recently, [3] proposed to introduce a regularization
term to enforce solution similarity among neighboring sources and/or
receivers.

Alternatively, a time-domain formulation of MDD was proposed in [28]: by
operating in the time domain, one can more naturally regularize the inverse
problem, alongside parametrizing the solution in a suitable transformed domain
(e.g., curvelets) where sparsity-promoting inversion can be employed. Along
similar lines, [29] and [16] proposed to introduce a preconditioner to enforce
reciprocity in the sought after Green’s function. Continuing, [31] suggested
to combine three physics-based preconditioners, namely causality, reciprocity,
and frequency–wavenumber locality within the same time-domain formulation.
This resulted in improved stability also for scenarios with highly complex
propagating overburdens (e.g., salt bodies). Finally, low-rank regularization
has been recently proposed by [13, 5, 6] as an alternative way to parametrize
the sought after solution. Whilst primarily focused on tackling the first
challenge (i.e., stability), this approach can also reduce the size of the
solution when a factorization-based form of low-rankness is considered.

From a computational standpoint, two main approaches have been investigated to
date to reduce the cost of applying MDD. The former leverages the fact that
the MDD equations can be formulated as a finite-sum functional [21]; as such,
the associated inverse problem can be solved by means of stochastic gradient
descent algorithms, where the gradient at each iteration is computed using a
small subset of randomly selected sources. This can reduce the overhead of
having to access the entire dataset at each iteration and has been shown to
lead to faster and more robust convergence. The latter involves the
application of block low-rank approximation techniques, such as the Tile-Low
Rank (TLR) algorithm [11, 15, 12, 20, 19], to the operator matrix. By doing
so, one can first reduce the size of the operator, and later accelerate
matrix-vector multiplications to performing this operation in two steps where
the factors, instead of the full matrix, are utilized.

Our paper extends the concept of block low-rank approximation to encompass
compression of both the operator, itself, and the right-hand side and unknown
matrix. This is particularly advantageous given the data-heavy nature of the
right-hand side and the unknowns in MDD. Furthermore, by embracing a more
sophisticated $\mathcal{H}^{2}$ [10, 9] approximation of the dense matrices,
we show that further compressibility of the wavefield can be achieved compared
to the TLR approach. Through this multifaceted approach, we endeavor to pave
the way for more efficient MDD implementations, ultimately expanding its
potential applications within the geophysical domain.

##  2 Multi-dimensional deconvolution

Let us consider a medium composed of two parts: the first one, usually
referred to as overburden, lies in between the surface of the Earth and the
receiver datum; the second one, below the receiver datum, represents the
target area we are ultimately interested to image. One or more sources,
denoted here as $x_{s}$, are placed anywhere in the overburden, whilst
receivers $x_{0}^{\prime}$ and virtual sources $x_{0}$ are located at the
datum that separates the overburden from the target area. We assume that the
wavefields $w^{-}$ from $x_{s}$ to $x_{0}$ and $w^{+}$ from $x_{s}$ to
$x_{0}^{\prime}$ are known (i.e., physically recorded or obtained by a pre-
processing step). In this problem, we are interested to recover the wavefield
$r$ that travels only below the receiver datum. Figure 1 shows a schematic
representation of the acquisition geometry associated with the MDD problem.

ground or water surface 1d or
2dsourcesreceivers$x_{s}$overburden$x^{\prime}_{0}$$x_{0}$target of
interestvirtual receivers$r$? Figure 1: Visual representation of the MDD
algorithm.

Following [36], the aforementioned wavefields are related to each other via a
convolution-type representation theorem,

$w^{-}(x_{s},x_{0},t)=\mathcal{F}_{t}^{-1}\int
w^{+}(x_{s},x^{\prime}_{0},f)\mathcal{F}_{t}\left(r(x^{\prime}_{0},x_{0},t)\right)dx^{\prime}_{0},$
(1)

where $\mathcal{F}_{t}$ and $\mathcal{F}_{t}^{-1}$ represent the forward and
inversion Fourier transforms, whilst $t$ and $f$ refer to the time and
frequency axes. In order to find a numerical solution of equation (1) for the
wavefield of interest $r$, both space and time are usually discretized,
leading to the following expression

$w^{-}(x_{s,i},x_{0,j},t_{k})=F_{t}^{-1}\left(\sum_{\alpha=1}^{N_{0}}w^{+}(x_{s,i},x^{\prime}_{0,\alpha},f_{p})F_{t}(r(x^{\prime}_{0,\alpha},x_{0,j},t_{k}))\right),$
(2)

where $F_{t}$ and $F_{t}^{-1}$ are the discrete forward and inversion Fourier
transforms, $i=1\dots N_{s}$, $j=1,\dots,N_{0}$, $\alpha=1,\dots,N_{0}$,
$k=1,\dots,N_{t}$, and $p=1,\dots,N_{f}$, with $N_{s}$, $N_{0}$, $N_{t}$,
$N_{f}$ representing the number of sources, receivers, time samples, and
frequency samples, respectively.

The linear system in equation (2) can be reformulated in a tensor form if we
assume that

$B_{t}[i,j,k]=w^{-}(x_{s,i},x_{0,j},t_{k})$
$W[i,\alpha,p]=W^{+}(x_{s,i},x^{\prime}_{0,\alpha},f_{p})$
$X_{t}[\alpha,j,k]=r(x^{\prime}_{0,\alpha},x_{0,j},t_{k})$

yielding

$F^{H}WFX_{t}=B_{t},$

where $F$ is a Discrete Fourier transform operator applied over the time
index, $W\in\mathbb{R}^{N_{s}\times N_{0}\times
N_{f}},X_{t}\in\mathbb{R}^{N_{0}\times N_{0}\times
N_{t}},B_{t}\in\mathbb{R}^{N_{s}\times N_{0}\times N_{t}}$.

By applying the Fourier transform to the system from right and left over the
time index, we obtain an equivalent system of equations in the frequency
domain:

$WX_{f}=B_{f},$

where $X_{f}\in\mathbb{R}^{N_{0}\times N_{0}\times N_{f}}$ and
$B_{f}\in\mathbb{R}^{N_{s}\times N_{0}\times N_{f}}$. This tensor system is
equivalent to a series of linear systems for different frequencies:

$W_{p}X_{p}=B_{p},\quad p=1\dots N_{f},$ (3)

where $W_{p}=W[:,:,p]$, $X_{p}=X_{f}[:,:,p]$, and $B_{p}=B_{f}[:,:,p]$.

$N_{s}\times N_{0}$$W_{p}$$N_{0}\times N_{0}$$=$$X_{p}$$N_{s}\times
N_{0}$$B_{p}$

The three matrices involved in equation 3 are usually dense and have
comparable sizes. In order to reduce the cost of solving these $N_{f}$ systems
of equations, one must find a suitable compressed representation of such
matrices.

##  3 Compression of the operator, right-hand side, and unknowns

In this section, we introduce several approaches that benefit from the global
or block-wise low-rank nature of the matrices in the system of equations (3).
These approaches span from the simplest low-rank approximation to more
advanced techniques based on the $\mathcal{H}^{2}$ matrix approximation.

###  3.1 Global low-rank parametrization

To begin, let us examine the scenario where all the matrices in equation (3)
possess the property of being globally low rank.

####  3.1.1 Approximation based on the non-symmetric low-rank parametrization
of $X_{p}$

Since we expect the solution $X_{p}$ to be globally low rank, let us write
$X_{p}$ in this factorized form:

$X_{p}=L_{p}R_{p},$

where $X_{p}\in\mathbb{R}^{N_{0}\times N_{0}}$,
$L_{p}\in\mathbb{R}^{N_{0}\times r}$, $R_{p}\in\mathbb{R}^{r\times N_{0}}$,
$r\ll N_{0}$. We obtain the system:

$W_{p}L_{p}R_{p}=B_{p}.$ (4)

One obvious option for solving this system of equations is represented by the
alternating least squares method (ALS [26, 27]), where the factors $L_{p}$ and
$R_{p}$ are updated in an alternating fashion. This method has been shown to
have poor convergence unless additional regularization terms are imposed on
the factors (see, for example, [5]). Instead, we decide here to compute also a
low-rank approximation for the matrix $B_{p}$:

$B_{p}=\widehat{L}_{p}\widehat{R}_{p},$

where $B_{p}\in\mathbb{R}^{N_{s}\times N_{0}}$,
$\widehat{L}_{p}\in\mathbb{R}^{N_{s}\times r}$,
$\widehat{R}_{p}\in\mathbb{R}^{r\times N_{0}}$, $r\ll N_{0},N_{s}$. We obtain
the system:

$W_{p}L_{p}R_{p}=\widehat{L}_{p}\widehat{R}_{p}.$ (5)

$N_{s}\times N_{0}$$W_{p}$$N_{0}\times r$$L_{p}$$r\times
N_{0}$$=$$R_{p}$$N_{s}\times r$$\widehat{L}_{p}$$r\times
N_{0}$$\widehat{R}_{p}$

We then assume that

$R_{p}=\widehat{R}_{p}.$ (6)

Such an assumption is based on the intuition that the solution and right-hand
side are wavefields with similar features and therefore can be represented by
a shared basis. The experiments in sections 4.1 and 4.2 confirm the validity
of this assumption as satisfactory solutions are produced in both scenarios.
Limitations and benefits of such an assumption are a subject of future
research.

Using the assumption in equation (6), equation (5) becomes:

$W_{p}L_{p}\widehat{R}_{p}=\widehat{L}_{p}\widehat{R}_{p}.$ (7)

Multiplying both sides of equation (7) by the pseudoinverse of
$\widehat{R}_{p}$ (i.e., $\widehat{R}_{p}^{\dagger}$) we obtain:

$W_{p}L_{p}=\widehat{L}_{p}$ (8)

$N_{s}\times N_{0}$$W_{p}$$N_{0}\times r$$L_{p}$$=$$N_{s}\times
r$$\widehat{L}_{p}$

We can now solve this system for $L_{p}$ and ultimately reconstruct our
solution of interest as $X_{p}=L_{p}\widehat{R}_{p}$. By choosing this
parametrization, we have reduced the size of the problem (and therefore the
number of computations required to obtain the solution) by a factor of
$\frac{N_{0}}{r}$. We refer to this technique as the USV method and as USV-rsp
when a reciprocity preconditioner is added in the solution of equation 8.

However, the proposed approach presents several downsides. Notably, it does
not exploit the fact that also $W_{p}$ may be a low-rank matrix. Moreover, the
chosen parametrization for the solution (i.e., $X_{p}=L_{p}\widehat{R}_{p}$)
ignores the fact that this matrix should be symmetric, owing to the
reciprocity of acoustic wavefields.

####  3.1.2 Approximation based on the symmetric low-rank parametrization of
$X_{p}$.

Let us now define the truncated SVD of the matrix $B_{p}$ as:

$B_{p}=U_{p}\Sigma_{p}V_{p},$

where $B_{p}\in\mathbb{R}^{N_{s}\times N_{0}}$,
$U_{p}\in\mathbb{R}^{N_{s}\times r}$, $\Sigma_{p}\in\mathbb{R}^{r\times r}$,
$V_{p}\in\mathbb{R}^{r\times N_{0}}$, $r\ll N_{0},N_{s}.$

Given that we want i) $X_{p}$ to be symmetric (but not Hermitian), and ii) the
matrices $B_{p}$ and $X_{p}$ to share their right factor, we write $X_{p}$ in
the following format:

$X_{p}=V_{p}^{\top}S_{p}V_{p},$

where $S_{p}\in\mathbb{R}^{r\times r}$ is an unknown matrix. We obtain the
system:

$W_{p}V^{\top}_{p}S_{p}V_{p}=U_{p}\Sigma_{p}V_{p}.$ (9)

$N_{s}\times N_{0}$$W_{p}$$N_{0}\times r$$V_{p}^{\top}$$r\times
r$$S_{p}$$r\times N_{0}$$=$$V_{p}$$N_{s}\times r$$U_{p}$$r\times
r$$\Sigma_{p}$$r\times N_{0}$$V_{p}$

Since $V_{p}$ and $U_{p}$ are unitary matrices (i.e.,
$V_{p}V_{p}^{{}^{\mathsf{H}}}=I_{r}$ and
$U_{p}^{{}^{\mathsf{H}}}U_{p}=I_{r}$), multiplying the system (9) by the
matrix $V_{p}^{{}^{\mathsf{H}}}$ from the right and with the matrix
$U_{p}^{{}^{\mathsf{H}}}$ from the left leads to

$\left(U_{p}^{{}^{\mathsf{H}}}W_{p}V_{p}^{\top}\right)S_{p}=\Sigma_{p},$ (10)

where all three matrices
$\left(U_{p}^{{}^{\mathsf{H}}}W_{p}V_{p}^{\top}\right),S_{p},\Sigma_{p}\in\mathbb{R}^{r\times
r}$.

$r\times r$$U_{p}^{{}^{\mathsf{H}}}W_{p}V_{p}^{\top}$$r\times
r$$S_{p}$$=$$r\times r$$\Sigma_{p}$

We can now solve this system for $S_{p}$ and ultimately reconstruct our
solution of interest as $X_{p}=V_{p}^{\top}S_{p}V_{p}$ and we refer to this
technique as the LR method. Similar to the previous section, we have here made
the assumption that the operator, right-hand side, and unknowns share the
bases $V_{p}$ and $U_{p}$. For the numerical problems in Sections 4.1 and 4.2,
this assumption is validated empirically. A theoretical analysis of the
limitations of such an assumption is the subject of future research.

To summarize, we observe that the parametrization
$X_{p}=V_{p}^{\top}S_{p}V_{p}$ can dramatically reduce the size of the
solution (i.e., from inverting a matrix of size $N_{0}\times N_{0}$ to
inverting a matrix of size $r\times r$) and, thus, the computational
complexity of the associated inverse problem. Impressively, this reduction in
complexity does not incur a significant compromise in the solution’s quality.
Furthermore, this approach empowers us to benefit from the symmetry of the
solution $X_{p}$, without requiring the introduction of any reciprocity
preconditioner in the solution of the inverse problem [31].

###  3.2 $\mathcal{H}^{2}$ approximation

Unfortunately, the matrices $W_{p}$, $X_{p}$, $B_{p}$ are usually not globally
low-rank in 3-dimensional problems. However, on the bright side, those
matrices can be well approximated using block low-rank structures.
Particularly, they can be represented in a TLR structure [11] as well as
$\mathcal{H}^{2}$ structure [10, 9], as we verify it in our numerical
experiments in Sections 4.3.

In this section, we introduce the fundamental concepts of the
$\mathcal{H}^{2}$ matrix approximation. For more comprehensive and formal
definitions, please refer to [9, 10, 4]. It is important to note that our
version for $\mathcal{H}^{2}$ decomposition differs from the standard one, yet
it remains valid. A thorough explanation of our version of $\mathcal{H}^{2}$
is available in [25].

Let us, for example, consider the dense right-hand side matrix $B_{p}$ in the
system (3) and build $\mathcal{H}^{2}$ factorization of it. Noting that rows
and columns of $B_{p}$ are associated with sources and receivers in the
physical acquisition geometry, we can start by hierarchically splitting the
grids of sources and receivers and storing these divisions within tree
structures. Utilizing this grid division, we rearrange the rows and columns of
the matrix $B_{p}$: such a rearrangement creates a block structure in the new
matrix, which possesses now block low-rank properties. Let the matrix $B_{p}$
have $N_{r}$ block rows and $N_{c}$ block columns. Thus, the matrix $B_{p}$ is
divided into blocks

$B_{ij}\in\mathbb{R}^{N_{i}\times N_{j}},\quad i\in 1\dots N_{r},j\in 1\dots
N_{c},$

where $N_{i}$ is a size of $i$-th block row, and $N_{j}$ is a size of $j$-th
block column. We also define block rows and block columns

$B_{i}\in\mathbb{R}^{N_{i}\times N_{0}},$ $B_{j}\in\mathbb{R}^{N_{s}\times
N_{j}}.$

All blocks in matrix $B_{p}$ can be separated into two classes: blocks that
correspond to the sources and receivers that are located closer than some
predefined distance (called ’close blocks‘) and all others (called ’far
blocks‘). Let us also define the block row $\widehat{B}_{i}$, which is $B_{i}$
with zeroed close blocks, and accordingly, block column $\widehat{B}_{j}$ is
$B_{j}$ with zeroed close blocks:

$B_{ij}\in\widehat{B}_{i}=\begin{cases}B_{ij}\in
B_{i}&\quad\mathbf{if~{}}B_{ij}\mathbf{~{}far}\\\
0&\quad\mathbf{if~{}}B_{ij}\mathbf{~{}close}\end{cases}$

$B_{ij}\in\widehat{B}_{j}=\begin{cases}B_{ij}\in
B_{j}&\quad\mathbf{if~{}}B_{ij}\mathbf{~{}far}\\\
0&\quad\mathbf{if~{}}B_{ij}\mathbf{~{}close}\end{cases}$

See illustration of $B_{i}$ and $\widehat{B}_{i}$ in Figure 2.

(a) Illustration of $B_{i}$

(b) Illustration of $\widehat{B}_{i}$

Figure 2: Example of $B_{i}$ and $\widehat{B}_{i}$. Dark green blocks
represent close blocks and light green blocks represent far blocks. White
blocks are zeros.

The key property of the $\mathcal{H}^{2}$ matrix formulation used in this work
is that block row $\widehat{B}_{i}$ and block column $\widehat{B}_{j}$
matrices are globally low rank. In other words, if we consider the left factor
$U_{i}\in\mathbb{R}^{N_{i}\times N_{i}}$ of the SVD decomposition of the
$\widehat{B}_{i}$ matrix, and multiply the block row $\widehat{B}_{i}$ from
the left with $U_{i}^{{}^{\mathsf{H}}}$, the resulting matrix will contain
non-zero values only in the first $r$ rows (whilst all other rows will have
values close to zero), where $r$ is rank of $\widehat{B}_{i}$. Similarly, if
we multiply the block row $B_{i}$ from the left with
$U_{i}^{{}^{\mathsf{H}}}$, the resulting matrix will have the same compression
for the far blocks and some changes in the close blocks. See the illustration
of the block row $U_{i}^{{}^{\mathsf{H}}}B_{i}$ in Figure 3.

Figure 3: Example of block row $U_{i}^{{}^{\mathsf{H}}}B_{i}$.

Analogously, we can compress the block columns by computing the right SVD
factor $V_{j}\in\mathbb{R}^{N_{j}\times N_{j}}$ of the block columns
$\widehat{B}_{j}$ and multiplying $V_{j}$ from the right to the block column
$B_{j}$.

Combining all the obtained $U_{i}$ and $V_{j}$ blocks into two big block-
diagonal unitary matrices, obtain:

$U=\begin{bmatrix}U_{1}&&\\\ &\ddots&\\\ &&U_{N_{b}}\\\ \end{bmatrix},\quad
V=\begin{bmatrix}V_{1}&&\\\ &\ddots&\\\ &&V_{N_{b}}\\\ \end{bmatrix},$

where $U\in\mathbb{R}^{N_{s}\times N_{s}}$, $V\in\mathbb{R}^{N_{0}\times
N_{0}}.$ Note, that in practice, it is more efficient to obtain $V_{j}$ factor
computing SVD not of the original $\widehat{B}_{j}$ block columns, but of the
block column $U^{{}^{\mathsf{H}}}B_{j}$ with excluded zeros.

Multiplying the the original matrix $B_{p}$ from left and right with the
matrices $U^{{}^{\mathsf{H}}}$ and $V^{{}^{\mathsf{H}}}$ we obtain a sparse
matrix $S_{B}\in\mathbb{R}^{N_{s}\times N_{0}}$:

$S_{B}=U^{{}^{\mathsf{H}}}B_{p}V^{{}^{\mathsf{H}}}.$ (11)

More specifically, $S_{B}$ is a block matrix with the following block
structure: it has a dense block
$S_{ij}=U_{i}^{{}^{\mathsf{H}}}B_{ij}V^{{}^{\mathsf{H}}}_{i}$, if the
corresponding block of $B_{p}$ is of close type and it has zero block with
non-zero $r\times r$ upper-left corner, if the corresponding block of $B_{p}$
is of far type. Note that in the actual implementation of our algorithm, we
store the matrix $S_{B}$ in a block-sparse format to avoid storage of the zero
elements. Multiplying the Equation (11) by the matrices $U$ and $V$, obtain:

$B_{p}=US_{B}V$

See Figure 4 for an illustration of the structure of our matrix’s
parametrization:

$N_{s}\times N_{0}$$B_{p}$$=$$N_{s}\times N_{s}$$U$$N_{s}\times
N_{0}$$S_{B}$$N_{0}\times N_{0}$$V$ Figure 4: Illustration of
$\mathcal{H}^{2}$ structure of matrix $B_{p}$.

Note that this is the one-level version of the $\mathcal{H}^{2}$ matrix. The
multi-level version is more advanced, and its application is the subject of
future research.

###  3.3 $\mathcal{H}^{2}$-like parametrization of the MDD

In this section, to simplify our notation, we drop the $p$ index from every
matrix. From now on, we therefore assume that the equations hold true for all
frequencies $p\in\mathbb{R}^{N_{f}}$.

As in Section 3.1.2, we know a priori that $X$ is symmetric (but not
hermitian). We also want the matrices $B$ and $X$ to share the right factor.
Thus, we will search $X$ in the following format:

$X=V^{\top}S_{X}V,$

where $S_{X}$ has the same format as $S_{B}$. We obtain the system:

$WV^{\top}S_{X}V=US_{B}V.$

Then we apply the matrices $U^{\mathsf{H}}$ to the $W$ from the left and
$V^{\mathsf{H}}$ from the right. We obtain the $U^{\mathsf{H}}WV^{\mathsf{H}}$
operator. Thanks to the physical properties of the problem, that operator has
the same structure as $S_{B}$:

$S_{W}=U^{\mathsf{H}}WV^{\mathsf{H}}.$

We obtain the following system:

$US_{W}\overline{V}V^{\top}S_{X}V=US_{B}V.$ (12)

Let us then multiply the system (12) with $U^{\mathsf{H}}$ from the left and
$V^{\mathsf{H}}$ from the right:

$U^{\mathsf{H}}US_{W}\overline{V}V^{\top}S_{X}VV^{\mathsf{H}}=U^{\mathsf{H}}US_{B}VV^{\mathsf{H}}.$

Considering the fact that $U$ and $V$ are hermitian matrices, and thus
$U^{\mathsf{H}}U=I$, $VV^{\mathsf{H}}=I$,
$\overline{V}V^{\top}=(VV^{\mathsf{H}})^{\top}=I^{\top}=I$, we receive:

$S_{W}S_{X}=S_{B}.$ (13)

We can now solve this system with the iterative method LSQR [18] for $S_{X}$.
After computing $S_{X}$, we reconstruct the solution using the formula

$X=V^{\top}S_{X}V.$

After that, using the Fourier transform, we obtain the final solution in the
time domain.

######  Remark 1.

While solving the system (13), the LSQR method employs matrix-by-matrix
multiplication. If the process does not preserve sparsity, the resulting dense
matrix negates the benefits of compression. The solution is to put to zero all
the elements of the resulting matrix that were zero in the initial matrices.
This approach is illustrated in Figure 5, which visually explains how we
preserve the sparse structure during multiplication.

$\times$$=$$\approx$$\approx$ Figure 5: The fill-in during the matrix
multiplication.

######  Remark 2.

In the numerical implementation, matrix-by-matrix multiplication is carried
out in a specific way: the right matrix is converted into a vector that
contains only its nonzero elements. The left matrix is implemented as a linear
operator, which is then applied to the vectorized right matrix. After applying
the linear operator, the resulting vector is converted back into a sparse
matrix format. This approach efficiently handles the multiplication by using
only the nonzero elements. Another significant benefit of this method is its
compatibility with the standard LSQR algorithm, which naturally operates with
vectors and linear operators. This compatibility makes integrating the LSQR
method with this approach easier.

This approach empowers us to benefit from the symmetry of the solution $X$,
enabling, for example, the utilization of the reciprocity precondition. We
refer to this technique as the $\mathcal{H}^{2}$ method.

##  4 Numerical Examples

In this section, we present three numerical experiments to validate the
proposed methods; these involve applying global low-rank (for the
2-dimensional case) and block low-rank (for the 3-dimensional case)
compression techniques to all elements of the linear system we want to solve.
The different experiments are chosen to mirror real-world seismic exploration
scenarios. Among the datasets, two are two-dimensional, offering standard
complexity scenarios, while the third is a significantly larger three-
dimensional dataset, providing a more challenging and comprehensive test of
the capabilities of our approaches in handling large-scale seismic data. All
algorithms are implemented in Python; our focus is to characterize the time
and memory savings of our methods in return for a controllable compromise in
the quality of the results.

###  4.1 Ocean-bottom cable (2d)

In this subsection, we consider a synthetic ocean-bottom cable (OBC) dataset.
In our specific case study, we use a two-dimensional segment of the SEG/EAGE
Overthrust model [1], which includes an additional water layer on top. 401
receivers are arranged on the sea floor. Each is spaced 10 meters apart and
292 meters below the ocean’s surface. 401 sources are placed at a depth of 20
m with a lateral spacing of 10 m.

Figure 6 illustrates the inversion results achieved using MDD, using the LR,
USV, USV-rsp (USV with reciprocity precondition), and $\mathcal{H}^{2}$
compression strategies described in the paper. These results are further
compared to the reference response directly modeled via finite difference
without the water column and free surface and to the "full" solution where we
solve MDD in the time or frequency domain with uncompressed matrices.

(a) Modeled solution

(b) Full, time dom.

(c) LR method

(d) USV method

(e) USV-rsp method

(f) $\mathcal{H}^{2}$ method

Figure 6: MDD inversion results for different solution methods.

The LR method outperforms the USV and USV-rsp methods, using a significantly
larger amount of data in the solution process. In this context, applying
reciprocity preconditioning has minimal impact on the outcome. In the 2D
scenario, the $\mathcal{H}^{2}$ method appears to be the least effective. Its
complex data structure necessitates large computational time, and the
additional approximations result in a visible accuracy loss. However, in 3D,
the large size of the matrices involved makes this format advantageous,
offsetting the initial drawbacks with its benefits, as we will show later.

Table 1 compares the time and memory requirements for the MDD methods
described above. We consider the uncompressed time and frequency domain MDD,
LR, USV, and $\mathcal{H}^{2}$ methods. For timing, we present integral
solution timings, including the approximation to the LR, USV, or
$\mathcal{H}^{2}$ format. For memory, we sum up the sizes of all matrices
involved in the MDD process (operator, right-hand side, and unknown matrices)
in the frequency domain.

Below, we present the formulas for the memory requirements. For this 2d ocean-
bottom cable dataset, the sizes of the arrays are: number of sources
$N_{s}=201$, number of receivers and virtual receivers $N_{0}=201$, number of
time samples $N_{t}=1150$, and number of frequency samples $N_{f}=200$, rank
is set to $r=50$.

The full-matrix right-hand side data in the time domain has a size of

$N_{s}\times N_{0}\times N_{t}\times 4\text{Bytes}=185.8\text{MB},$ (14)

whilst the size of the full-matrix right-hand side data in the frequency
domain is

$N_{s}\times N_{0}\times N_{f}\times 8\text{Bytes}=64.6\text{MB},$ (15)

the size of the LR right-hand side data is

$N_{s}\times r\times N_{f}\times 8\text{Bytes}=16.1\text{MB},$ (16)

the size of the USV right-hand side data is

$r\times r\times N_{f}\times 8\text{Bytes}=4.0\text{MB},$ (17)

and the size of the $\mathcal{H}^{2}$ right-hand side data is computed
empirically as $13.5\text{MB}$.

The full-matrix unknowns in the time domain have a size of

$N_{0}\times N_{0}\times N_{t}\times 4\text{Bytes}=185.8\text{MB},$ (18)

whilst the size of the full-matrix unknowns in the frequency domain is

$N_{0}\times N_{0}\times N_{f}\times 8\text{Bytes}=64.6\text{MB},$ (19)

the size of the LR unknowns is

$N_{0}\times r\times N_{f}\times 8\text{Bytes}=16.1\text{MB},$ (20)

the size of the USV unknowns is

$r\times r\times N_{f}\times 8\text{Bytes}=4.0\text{MB},$ (21)

and the size of the $\mathcal{H}^{2}$ unknowns side data is again computed
empirically as $13.5\text{MB}$.

The operator is always stored in the frequency domain. The operators have the
same size as the frequency domain right-hand side of the same method for all
methods except for LR, where operators have the same size as the full-matrix
right-hand side.

| full, time dom. | full, freq. dom | LR | USV | $\mathcal{H}^{2}$
---|---|---|---|---|---
Time (s) | 21.3 | 12.8 | 6.0 | 0.71 | 33.8
Mem (MB) | 436.3 | 193.9 | 96.8 | 12.0 | 40.4
Table 1: Time and memory comparison of the candidate methods.

In this dataset, all three system matrices across every frequency possess a
globally low rank. This characteristic underscores the effectiveness of the
USV method, which emerges as the optimal choice in terms of both memory and
time efficiency.

###  4.2 Target level below a complex salt body (2d)

For our second experiment, we use a dataset created by redatuming surface
reflection data at a target level below a complex salt body [31]. These
wavefields are obtained by applying the scattering Rayleigh-Marchenko method
as detailed in [32]. This method allows one to estimate wavefields traveling
upwards and downwards as if they were emitted from surface sources and then
captured by virtual receivers located under a complex layer of earth, like
beneath a salt body; see Figure 7.

The experimental area extends over 16.26 km horizontally and reaches a depth
of 8 km. It is designed to imitate a challenging environment with an uneven
salt layer and other complex features like sharp changes in rock layers,
uneven sediment layers, and inconsistencies between different rock layers.
Redatuming in such a complex area is difficult because of the stark
differences in rock properties, which can cause multiple internal reflections.
These reflections can lead to misleading signals in the data, often seen as
artifacts or distortions. Additionally, other challenges like limited
frequency range and insufficient source illumination make the redatuming
process more challenging.

Figure 7: Sub-salt velocity model supporting decomposed wavefields at the
receiver array. A magenta line runs horizontally at a depth of 4.4 kilometers
below the layer above the salt. The sources of the wavefields, shown as a
solid red line, are positioned at the surface. A dashed black line also marks
the specific target area where we aim to capture and analyze the reflected
wave signals.

Figure 8 presents a reference solution, and Figure 9 presents different
results obtained by applying MDD inversion to the considered dataset. These
results were achieved with the various methods detailed in the paper. This
particular dataset benefits from having a globally low-rank structure across
all three matrices in systems for all frequencies. However, despite this
advantage, the problem is inherently unstable. The results are only reliable
when the reciprocity preconditioner is applied, as clearly demonstrated in
Figure 9. (Compare same methods with and without preconditioning.)

(a) Modeled solution

(b) Full, time dom

(c) Full-rcp, time dom

Figure 8: MDD inversion results for reference solution.

(a) LR method

(b) USV method

(c) USV-rsp method

(d) $\mathcal{H}^{2}$ method

(e) $\mathcal{H}^{2}$-rsp method

Figure 9: MDD inversion results for different solution methods.

Table 2 compares the time and memory requirements for the MDD methods
described above. We consider the uncompressed time and frequency domain MDD,
LR, USV, and $\mathcal{H}^{2}$ methods. For timing, we present integral
solution timings, including the approximation to the LR, USV, or
$\mathcal{H}^{2}$ format. For memory, we sum up the sizes of all matrices
involved in the MDD process (operator, right-hand side, and unknown matrices)
in the frequency domain.

We compute the memory requirements using the formulas (14 – 21) presented in
the previous subsection. The sizes of the $\mathcal{H}^{2}$ unknowns side data
are computed empirically. In this 2d target level below a complex salt body
dataset, the number of sources is $N_{s}=201$, the number of receivers and
virtual receivers is $N_{0}=151$, the number of time samples is $N_{t}=2001$,
the number of frequency samples is $N_{f}=400$, and rank is set to $r=50$.

| full, time dom. | full, freq. dom. | LR | USV | $\mathcal{H}^{2}$
---|---|---|---|---|---
Time (s) | 26.5 | 19.3 | 10.8 | 2.1 | 19.5
Mem (MB) | 583.0 | 267.2 | 153.4 | 24.0 | 80.7
Table 2: Time and memory comparison of the candidate methods.

For this dataset, the reciprocity preconditioning is crucial to obtaining
valuable results. Thus, the USV-rsp method outperforms other methods not only
in terms of memory and time but also in terms of solution quality. The LR
method’s inability to use reciprocity becomes critical in this example. As in
the previous example, the $\mathcal{H}^{2}$ method is the least effective for
the above reasons.

###  4.3 Ocean-bottom cable (3d)

In this subsection, we consider the synthetic ocean-bottom cable model. In
this example, we use a full 3d SEG/EAGE Overthrust model [1], which includes
an additional water layer on top.

The grid of $177\times 90$ receivers is arranged on the sea floor. Each spaced
20 meters apart and 292 meters below the ocean’s surface. Additionally, we
create a $120\times 217$ grid of shot gathers with a source spacing of 20 m
and a depth of 10 m.

Figure 10 presents a reference MDD solutions, and Figure 11 illustrates the
visualization of inversion results achieved using frequency MDD, applied
through method $\mathcal{H}^{2}$ with various rank parameters, compared to the
reference results: the modeled solution and the frequency domain full-matrices
solution. This dataset, a substantial 3D one, does not exhibit a globally low
rank across its three matrices at all frequencies. Consequently, the only
feasible method for this dataset is the $\mathcal{H}^{2}$ method. Our analysis
of the previous two datasets shows that the $\mathcal{H}^{2}$ method incurs
significant overhead compared to the LR and USV methods. However, it remains
preferable over using full matrices.

(a) Modeled local reflectivity

(b) Solution $\mathbf{R_{inv}}$, computed in frequency domain, for full
matrices

Figure 10: Reference solutions

(a) $\mathbf{R_{inv}}$, $\mathcal{H}^{2}$ , $r=200$

(b) $\mathbf{R_{inv}}$, $\mathcal{H}^{2}$ , $r=100$

(c) $\mathbf{R_{inv}}$, $\mathcal{H}^{2}$ , $r=1$

Figure 11: MDD results obtained with $\mathcal{H}^{2}$ method for various
ranks.

In Figure 12, we demonstrate a comparison between the solutions obtained
through the full-matrix method and the $\mathcal{H}^{2}$ method within the
frequency domain. The frequency slice number is $p=90$.

(a) Full freq. dom.

(b) $\mathcal{H}^{2}$ , $r=200$

(c) $\mathcal{H}^{2}$ , $r=100$

(d) Full freq. dom., fragment

(e) $\mathcal{H}^{2}$ , $r=200$, fragment

(f) $\mathcal{H}^{2}$ , $r=100$, fragment

Figure 12: MDD solution $X_{p}$ in the frequency domain for frequency slice
number $p=90$.

This comparison reveals that although the $\mathcal{H}^{2}$ method may
represent a compromise in terms of outright solution fidelity, it notably
succeeds in maintaining all essential elements of the solution.

Table 3 compares the time and memory requirements for the MDD methods
described above. We consider the uncompressed time and frequency domain MDD
and $\mathcal{H}^{2}$ method for various ranks. Given the computational
limitations of fully resolving the time domain MDD on even the most advanced
supercomputers, we are confined to presenting only theoretically derived
memory requirements. Conversely, in the frequency domain, leveraging a state-
of-the-art supercomputer equipped with parallel GPU cores has enabled us to
calculate and visually present the results. However, it is important to note
that comparisons of computation time may not be equitable, as the
$\mathcal{H}^{2}$ method calculations were performed on a CPU-core-based
supercomputer.

For timing, we present integral solution timings, including the approximation
to the $\mathcal{H}^{2}$ format. For memory, we sum up the sizes of all
matrices involved in the MDD process (operator, right-hand side, and unknown
matrices) in the frequency domain. We compute the memory requirements using
the formulas (14 – 21). The sizes of the $\mathcal{H}^{2}$ unknowns side data
are computed empirically. In this 3d Ocean-bottom cable dataset, the number of
sources $N_{s}=26040$, number of receivers and virtual receivers
$N_{0}=15930$, number of time samples $N_{t}=1126$, and number of frequency
samples $N_{f}=200$.

| full, time dom. | full, freq. dom. |  $\mathcal{H}^{2}$ , $r=200$ |  $\mathcal{H}^{2}$ , $r=100$ |  $\mathcal{H}^{2}$ , $r=1$
---|---|---|---|---|---
Time (s) | - | - | 4189 | 1530 | 350
Mem (TB) | 3.85 | 2.2 | 0.46 | 0.11 | 0.017
Table 3: Time and memory comparison of the presented methods.

The trade-off becomes apparent: reducing the rank compromises solution quality
(as illustrated in Figures 11 and 10), but it significantly saves time and
memory resources.

##  5 Conclusion

To advance the solution of linear systems that involve massive right-hand
sides and unknowns, this paper addresses the challenge of solving systems with
a data-heavy right-hand side and solution, using the example of
Multidimensional Deconvolution (MDD) inversion. We introduce and compare
several approaches to compress all matrices within the linear system,
significantly improving the solution efficiency, including algorithms based on
global low-rank (USV, LR) and block low-rank ($\mathcal{H}^{2}$).

##  6 Acknowledgements

Authors are thankful to the Supercomputing Laboratory and the Extreme
Computing Research Center at King Abdullah University of Science and
Technology for their support and computing resources.

## References

## References

  * [1] F. Aminzadeh, J. Brac, and T. Kunz.  3d salt and overthrust models.  In SEG/EAGE Modeling Series, Tulsa, Oklahoma, 1997. Society of Exploration Geophysicists.
  * [2] A. Bakulin and R. Calvert.  The virtual source method: Theory and case study.  Geophysics, 71(4):SI139–SI150, 2006.
  * [3] D. Boiero and C. Bagaini.  Up and down deconvolution in complex geological scenarios.  Geophysics, 86(5):WC55–WC65, 2021.
  * [4] S. Börm.  Efficient numerical methods for non-local operators: H2-matrix compression, algorithms and analysis, volume 14.  European Mathematical Society, 2010.
  * [5] F. Chen, M. Ravasi, and D. Keyes.  A physics-aware, low-rank regularization for multidimensional deconvolution.  arXiv preprint arXiv:2312.11004, 2023.
  * [6] F. Chen, M. Ravasi, and D. Keyes.  Solving multi-dimensional deconvolution via a nuclear-norm regularized least-squares approach.  In 84th EAGE Annual Conference & Exhibition, volume 2023, pages 1–5. European Association of Geoscientists & Engineers, 2023.
  * [7] D. Draganov, X. Campman, J. Thorbecke, A. Verdel, and K. Wapenaar.  Reflection images from ambient seismic noise.  Geophysics, 74(5):A63–A67, 2009.
  * [8] J. T. Fokkema and P. M. van den Berg.  Seismic applications of acoustic reciprocity.  Elsevier, 2013.
  * [9] L. Greengard and V. Rokhlin.  A fast algorithm for particle simulations.  J. Comput. Phys., 73(2):325–348, Dec. 1987.
  * [10] W. Hackbusch, B. Khoromskij, and S. Sauter.  On $\mathcal{H}^{2}$-matrices.  In H.-J. Bungartz, et al. (eds.), Lectures on Applied Mathematics, pages 9–30. Springer-Verlag, Berlin Heidelberg, 2000.
  * [11] Y. Hong, H. Ltaief, M. Ravasi, L. Gatineau, and D. E. Keyes.  Accelerating seismic redatuming using tile low-rank approximations on NEC SX-aurora TSUBASA.  2021\.
  * [12] Y. Hong, H. Ltaief, M. Ravasi, and D. E. Keyes.  HPC seismic redatuming by inversion with algebraic compression and multiple precisions.  2022\.
  * [13] R. Kumar, D. Boiero, C. Bagaini, and M. Vassallo.  Transform-domain multidimensional deconvolution–sparsity v/s low-rank.  In 83rd EAGE Annual Conference & Exhibition, volume 2022, pages 1–5. European Association of Geoscientists & Engineers, 2022.
  * [14] X. Li, T. Becker, M. Ravasi, J. Robertsson, and D.-J. van Manen.  Closed-aperture unbounded acoustics experimentation using multidimensional deconvolution.  The Journal of the Acoustical Society of America, 149(3):1813–1828, 2021.
  * [15] H. Ltaief, Y. Hong, L. Wilson, M. Jacquelin, M. Ravasi, and D. E. Keyes.  Scaling the “memory wall” for multi-dimensional seismic processing with algebraic compression on Cerebras CS-2 systems.  In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–12, 2023.
  * [16] N. Luiken and T. Van Leeuwen.  Seismic wavefield redatuming with regularized multi-dimensional deconvolution.  Inverse Problems, 36(9):095010, 2020.
  * [17] S. Minato, T. Matsuoka, and T. Tsuji.  Singular-value decomposition analysis for seismic interferometry by multidimensional deconvolution.  In SEG Technical Program Expanded Abstracts 2011, pages 2694–2699. Society of Exploration Geophysicists, 2011.
  * [18] C. C. Paige and M. A. Saunders.  LSQR: An algorithm for sparse linear equations and sparse least squares.  ACM Transactions on Mathematical Software (TOMS), 8(1):43–71, 1982.
  * [19] M. Ravasi, Y. Hong, H. Ltaief, and D. Keyes.  Tile-low rank compressed multi-dimensional convolution and its application to seismic redatuming problems.  In 83rd EAGE Annual Conference & Exhibition, volume 2022, pages 1–5. European Association of Geoscientists & Engineers, 2022.
  * [20] M. Ravasi, Y. Hong, H. Ltaief, D. Keyes, and D. Vargas.  Large-scale Marchenko imaging with distance-aware matrix reordering, tile low-rank compression, and mixed-precision computations.  In Second International Meeting for Applied Geoscience & Energy, pages 2606–2610. Society of Exploration Geophysicists and American Association of Petroleum Geologists, 2022.
  * [21] M. Ravasi, T. Selvan, and N. Luiken.  Stochastic multi-dimensional deconvolution.  IEEE Transactions on Geoscience and Remote Sensing, 60:1–14, 2022.
  * [22] E. Ruigrok, X. Campman, D. Draganov, and K. Wapenaar.  High-resolution lithospheric imaging with seismic interferometry.  Geophysical Journal International, 183(1):339–357, 2010.
  * [23] G. Schuster, J. Yu, J. Sheng, and J. Rickett.  Interferometric/daylight seismic imaging.  Geophysical Journal International, 157(2):838–852, 2004.
  * [24] G. T. Schuster and M. Zhou.  A theoretical overview of model-based and correlation-based redatuming methods.  Geophysics, 71(4):SI103–SI110, 2006.
  * [25] D. A. Sushnikova and I. V. Oseledets.  Simple non-extensive sparsification of the hierarchical matrices.  Advances in Computational Mathematics, 46(4):52, 2020.
  * [26] R. Tauler and D. Barceló.  Multivariate curve resolution applied to liquid chromatography—diode array detection.  TrAC Trends in Analytical Chemistry, 12(8):319–327, 1993.
  * [27] R. Tauler, A. Smilde, and B. Kowalski.  Selectivity, local rank, three-way data analysis and ambiguity in multivariate curve resolution.  Journal of Chemometrics, 9(1):31–58, 1995.
  * [28] J. van der Neut and F. J. Herrmann.  Interferometric redatuming by sparse inversion.  Geophysical Journal International, 192(2):666–670, 2013.
  * [29] J. van der Neut, M. Ravasi, Y. Liu, and I. Vasconcelos.  Target-enclosed seismic imaging.  Geophysics, 82(6):Q53–Q66, 2017.
  * [30] J. van der Neut, J. Thorbecke, K. Mehta, E. Slob, and K. Wapenaar.  Controlled-source interferometric redatuming by crosscorrelation and multidimensional deconvolution in elastic media.  Geophysics, 76(4):SA63–SA76, 2011.
  * [31] D. Vargas, I. Vasconcelos, M. Ravasi, and N. Luiken.  Time-domain multidimensional deconvolution: A physically reliable and stable preconditioned implementation.  Remote Sensing, 13(18):3683, 2021.
  * [32] D. Vargas, I. Vasconcelos, Y. Sripanich, and M. Ravasi.  Scattering-based focusing for imaging in highly complex media from band-limited, multicomponent data.  Geophysics, 86(5):WC141–WC157, 2021.
  * [33] K. Wapenaar and J. Fokkema.  Green’s function representations for seismic interferometry.  Geophysics, 71(4):SI33–SI46, 2006.
  * [34] K. Wapenaar, E. Slob, and R. Snieder.  Seismic and electromagnetic controlled-source interferometry in dissipative media.  Geophysical prospecting, 56(3):419–434, 2008.
  * [35] K. Wapenaar, J. Van Der Neut, E. Ruigrok, D. Draganov, J. Hunziker, E. Slob, J. Thorbecke, and R. Snieder.  Seismic interferometry by crosscorrelation and by multidimensional deconvolution: A systematic comparison.  Geophysical Journal International, 185(3):1335–1364, 2011.
  * [36] K. Wapenaar, J. Van Der Neut, E. Ruigrok, D. Draganov, J. Hunziker, E. Slob, J. Thorbecke, and R. Snieder.  Seismic interferometry by crosscorrelation and by multidimensional deconvolution: A systematic comparison.  Geophysical Journal International, 185(3):1335–1364, 2011.