#  Stressed GANs snag desserts
a.k.a
Spotting Symmetry Violation
with Symmetric Functions

Christopher G. Lester  Cavendish Laboratory, University of Cambridge, CB3 0HE,
United Kingdom  Rupert Tombs  Cavendish Laboratory, University of Cambridge,
CB3 0HE, United Kingdom

(1st November 2021)

###### Abstract

Is parity violated at the LHC in some way that no-one has anticipated? This
paper is intended to illustrate a simple way one could try to answer that
question without relying on theoretical predictions from specific mechanisms
of parity violation. Alternatively the method could be used as a calibration
tool to remove parity-violating misalignments in imperfect detectors. A longer
sister paper extends the concepts to general symmetries (beyond parity) and
addresses issues connected to acceptance and reliability of real detectors.
The work here deliberately ignores those issues and works in an idealised
environment to allow greater simplicity of exposition.

##  1 Introduction

The core principles underlying the whole of this paper are:

  1. 1.

if a probability distribution $p(x)$ assigns equal probability densities to
objects $x$ and their mirror images $Px$ then the expectation of a parity-odd
function $f(x)$ on such a distribution will be zero; and therefore

  2. 2.

it is only possible for a parity-odd function’s expectation to be non-zero (to
a statistically significant degree) if the sample of events on which it was
evaluated violated parity – i.e. had some objects more or less frequently then
their mirror images.

Therefore, if one wishes to look for parity violation one can attempt to train
a parity-odd neural net (or any other sufficiently general function which is
unavoidably parity-odd by construction)111Note that it is absolutely trivial
to force a neural net (or any function) by-construction to be parity-odd. You
only need to evaluate it twice, once on $x$ and once on the mirror image of
$x$, and then subtract one answer from the other. I.e. if $g(x)$ is a general
net (or function) then the function $f(x)$ defined by $f(x)\equiv g(x)-g(Px)$
is automatically parity-odd. to give a non-zero mean on one portion of a real
dataset. The evidence for parity violation would then be provided if the
trained model (when evaluated on the remaining test portion of the dataset)
gives a statistically significant non-zero mean.

###  1.1 Discussion

It cannot be stressed how different this pattern of net use is to those most
frequently seen in particle physics. There are many ways one could illustrate
this, but perhaps the most easy one to highlight is that the training step is
(technically) optional!

If the training step were entirely omitted, the net would be in a random state
whose only property was its (unavoidable) characteristic of being parity-odd.
Nonetheless, if such a net were able to generate a statistically non-zero mean
for data in the test sample, then that non-zero mean would be evidence for
parity violation! In practice, of course, the training step very useful.
Without training there is little chance that the net will actually be
attending to any useful features of the data. But the point we wish to stress
is that the training process exists not to make the discovery of parity-
violation possible, but only to increase the probability of being able to
claim discovery.

For related reasons ‘overtraining’ does not have the same negative overtones
it would normally have. The possibility of overtraining still exists but it
cannot lead to an accidental mis-discovery of parity-violation as there simply
is not a concept of ‘right’ or ‘wrong’ training, only of ‘more efficient’ or
‘less efficient’ training. All that matters is that the net (whether trained
or not) manages to show a statistically significant non-zero mean on its test
data’ – if it does, it’s a good net, however it got there. The worst that can
happen is that in training the net focussed on entirely the wrong thing,
overtrains, and so does not know how to bias the test-set away from zero.

We note that it is also unusual that the test and training datasets are
identical, consist entirely of data (no Monte Carlo) and are unlabelled – and
so are immediately usable in high luminosity and intense enviroments where
simulation can be both complex and costly.

Finally we note that the discussion, though couched in terms of neural nets
could in principle use any function approximators.222Whether or not the method
could benefit from using Generative Adversarial Networks (GANs) has not been
seriously considered by the authors, but the inclusion of GANs in the title
seemed too good an opportunity to pass on, given the importance of baked-in
symmetries.

###  1.2 No free lunch – but win-win outcome

No detector is perfectly aligned or calibrated or without regions where
acceptance can vary. Consequently evidence of what we have called ‘parity
violation’ cannot be pinned immediately at nature’s door. On the contrary, it
is quite possible that the biggest sources of parity violation could be mis-
alignments and mis-calibrations. This is not really a disadvantage of the
method, though. Instead it is a sign that if you can get a non-zero mean on
the test-sample, then you know immediately that either you don’t understand
something about your detector, or you have an interesting discovery. Both are
equally important, so a signal is a win-win outcome even if the source of the
parity-violation is not immediately obvious.

###  1.3 Working in the real world

Real detectors and experiments are, of course, imperfect. Not all parts work
with perfect efficiency, and acceptance does not extend uniformly into all
directions. A sister paper [7] addresses those absolutely critical questions,
and also explains how to generalise the symmetry-violation discovery process
to general symmetries (including continuous symmetries) rather than the narrow
case of parity discussed herein.

The current paper, however, deliberately eschews those details as the
intention here is to provide a simple and pedagogically useful illustration of
the key message, free from other distractions.

##  2 The MNIST Calorimeter

###  2.1 Visualising calorimeter deposits

Calorimetric deposits in idealised versions of detectors like ATLAS or CMS are
commonly represented as functions on a cylinder about the beam axis as
depicted in the upper left corner of Figure 1. This cylinder may be opened out
into a rectangle with pseudorapidity $\eta$ on one axis and azimuthal angle
$\phi$ on the other as depicted in the lower left corner of the same figure.
In this flattened view the upper and lower edges are identified with each
other.

Figure 1: This figure aims to show how transformations of a ‘real’ three-
dimensional of a calorimeter (top row) correspond to induced transformations
in the unrolled two-dimensional $(\eta,\phi)$-space in which those deposits
are often represented. The first column shows deposits that would be typical
of a five jet event. The second column shows an easier-to-visualise set of
deposits resembling a letter ‘R’. The other three columns show how the
deposits of the first ‘R’ would move around if affected by (i) an arbitrarty
rotation about the $z$-asis, (ii) a rotation of 180 degrees about the $x$-axis
and (iii) a parity inversion.

Real calorimetric deposits are blotchy and messy, and this makes it hard for
humans to picture what they might look like when rotated or mirrored. Since
the goals of this paper are pedagocial, we choose to imagine a strange world
in which calorimetric deposits look just like hand-written letters of the
alphabet. For example, the second column of Figure 1 imagines what would
happen if calorimeter deposits happened to have landed in locations which
(once unfolded) would resemble the capital letter ‘R’. Making use of this
visual simplification, the other columns of Figure 1 illustrate the effect
that different transformations would have on the first ‘R’.

###  2.2 The LHC datasets

(a) A mode-1-chiral LHC dataset

(b) A mode-1-achiral LHC dataset

(c) A mode-2-chiral LHC dataset

(d) A mode-2-achiral LHC dataset

Figure 2: Four types of LHC dataset are shown. Further descriptions can be
found in the body-text.

A real analysis would use real data not Monte Carlo simulations.
Unfortunately, as there is no real detector that produces such energy deposits
that look like letters, we must generate some toy data of our own for our
pedagogical purposes. There is a ready supply of hand written numerals in the
MNIST dataset [2] so we implement a generator that draws its calorimeter
deposits from there. We call these Large datasets of Handwritten Calorimeter
digits – or ‘LHC datasets’ for short. We produce LHC datasets in four modes as
shown in Figure 2 and defined as follows:

  * •

mode-1-chiral : These are normal MNIST digits. These samples are chiral
because humans write hand-written digits with handed-appendages called hands.
The word chiral even come from the ancient Greek word for hand:
$\chi\epsilon\iota\rho$;

  * •

mode-1-achiral : These are MNIST digits which have been reversed left-to-right
with probability 50%; The symmetrisation process makes these samples achiral.

  * •

mode-2-chiral : These are random pairs of normal MNIST digits superimposed on
each other. They inherit their chirality from their constituent digits.

  * •

mode-2-achiral : These are random pairs of normal MNIST superimposed on each
other with one member of the pair forwards and one reversed. The
symmetrisation process makes these samples achiral.

Mode-1-achiral acts as a control for mode-1-chiral, and mode-2-achiral acts as
a control for mode-2-chiral. The only purpose of mode-2 is to illustrate that
the analysis is not dependent on there being only a small number of ‘types’ of
deposit — in this case ten as there are ten sorts of digits. Both types of
Mode-2 therefore act as a form of control against the possibility that mode-1
is too simple to be a good proxy for real and highly variable calorimeter
data.

###  2.3 Removing remaining differences between real calorimeter data and
MNIST derived calorimeter analogues

MNIST digits are ‘oriented’ because the humans who originally generated them
were instructed to write them in boxes the conventional way up.

Real jet data, in contrast, is not oriented: if energy deposits were to to
appear in a configuration resembling a capital ‘R’ then the deposits would be
as likely to arrive in any of the configurations shown in columns two, three
and four of Figure 1, provided that the collider has symmetric beams and no
reason to choose any $\phi$-orientation more than any other.

To avoid the MNIST-based dataset having an unfair handle (orientability) not
possessed by real jet data, this discrepancy needs to be removed. One way of
removing it would be to augment our dataset by adding $y$-translations of
uniformly random magnitude to every image in an LHC dataset (being sure to
wrap pixels lost from the top back onto the bottom) before then giving every
image a 180 degree rotation with probability 50%.

Although one could augment the dataset in that way, it is more efficient to
simply make the neural net invariant with respect to those transformations.
Doing so is as easy as making the net parity-odd, the only difference being
that one seeks symmetry rather than antisymmetry and so may combine sub-nets
with symmetric operations (like maxpooling or summing) rather than
antisymmetric operations (like subtractions).

##  3 Our parity-odd axisymmetric network

Recall that the necessary axisymmetry comes from establishing (i) an
invariance with respect to $y$-translations of 28x28-pixel images
corresponding to $\phi$-invariance in the real world, and (ii) invariance with
respect to 180 degree rotations of the 28x28 images corresponding to
interchange of beams and isotropy.

Our final network $f(x)$ gains its parity-oddness and one-half of its
axisymmetry (the part requiring invariance with respect to 180 degree
rotations) by being defined in terms of a sub-network $g(x)$ as follows:

$\displaystyle f(x)$ $\displaystyle=g(x)+g(R_{180}x)-g(Px)-g(R_{180}Px).$ (1)

The sub-network $g(x)$ (implemented in pytorch [6]) therefore only needs to
respect the $\phi$-rotation (or $y$-translation) invariance. The $g(x)$ sub-
network starts by having two 2D-convolutional network layers (one with kernel
size $4\times 4$ in $\eta\times\phi$, and the other with size $14\times 6$ in
$\eta\times\phi$, always stride 1) to provided basic image processing of the
28x28 pixel images fed to it. These two layers bring the number of features up
to 32 and then to 64. They are then followed by a max-pooling over the $\phi$
direction to finally achieve the desired $\phi$-invariance. $g(x)$ then
concludes with two fully connected layers reducing 640 features to 128 and
thence to a single feature so that $f(x)$ can then make used of $g(x)$ in the
manner outlined in (1). The package used (pytorch, [6]) self-describes the
structure of our net as follows:

    Net(
      (conv1): Conv2d(1, 32, kernel_size=(4, 4), stride=(1, 1))
      (conv2): Conv2d(32, 64, kernel_size=(14, 6), stride=(1, 1))
      (dropout1): Dropout2d(p=0.25, inplace=False)
      (dropout2): Dropout2d(p=0.5, inplace=False)
      (fc1): Linear(in_features=640, out_features=128, bias=True)
      (fc2): Linear(in_features=128, out_features=1, bias=True)
    )

The network is trained on images provided in batches. Each batch

$B=\\{x_{1},x_{2},\ldots,x_{|B|}\\}$

contains $|B|=64$ images. For each batch of the mean net score $\mu_{B}$ is
calculated:

$\mu_{B}=\frac{1}{|B|}\sum_{x\in B}f(x).$

The standard deviation of the scores of the images in the batch is also
recorded as $\sigma_{B}$. The net is trained to maximise $\mu_{B}/\sigma_{B}$.
This process encourages to the net to do its best to give every image a common
constant positive value, and to eschew negative values if possible. The
complete code for the network may be found in [4].

##  4 Results

(a) A mode-1-chiral LHC dataset

(b) A mode-1-achiral LHC dataset

(c) A mode-2-chiral LHC dataset

(d) A mode-2-achiral LHC dataset

Figure 3: Four types of LHC dataset from the test dataset are shown after
having been coloured according to the output of the net. Bright green
represents increasingly positive numbers. Shades of increasingly vibrant red
represent increasingly negative numbers. Yellow represents zero.

The preponderance of green in plots (a) and (c) within Figure 3 shows that the
networks trained on the chiral modes are able to label the majority of their
inputs with positive scores, meaning that they have identified that these
datasets are chiral (i.e. violate parity). In contrast the nets have not
managed to achieve this for the achiral controls ((b) and (d) in the same
figure), despite having been trained on achiral control data. It is notable
that the networks are not flummoxed by the differences between mode-1 and
mode-2.

Mode | $\rho^{+}_{\mathrm{chiral}}$ | $\rho^{+}_{\mathrm{achiral}}$
---|---|---
mode–1 | 90$\pm$1 % | 51$\pm$2 %
mode–2 | 89$\pm$4 % | 42$\pm$6 %
Table 1: When tested on groups of 1024 images in the relevant test dataset,
the above numbers show the typical fraction $\rho^{+}$ of images per batch
which were assigned a positive value by trained networks. As desired, the
chiral LHC datasets are flagged as chiral at a statistically significant
level, while those LHC datasets which are achiral are not, with their
fractions remaining consistent with the 50% expected by chance alone.

Quantitative results are presented in Table 1 and confirm that parity
violation is discovered where it should be, and is not discovered where it
should not be.

Figure 4: Here the same 64 letters from a mode-1-chiral LHC dataset are
evaluated by the network in six different ways by the same trained network. In
the first pair of rows they are evaluated without any changes at all, and most
are coloured green as expected. In the next pair of rows they are re-evaluated
but after a shift (or roll) down by 10 pixels. In the next two rows the
transformation prior to evaluation is instead a 180 degree rotation. The key
take home message from these first six rows is that all of these
transformations do not affect the network score, and so the network is
axisymmetric as desired. The lower six rows repeat the upper six row, but
after an initial left-right parity flip. In this case we see that there is (as
desired) an immediate flip in the output from positive to negative (or in rare
cases vice versa) because the network is, by construction, parity-odd.

Figure 4 demonstrates that our network function is indeed parity-odd and
axisymmetric as desired.

##  5 Discussion.

The strategy has achieved its goals. What is perhaps more interesting to
reflect on is what this buys in comparison to other methods for analysing
multi-jet events for evidence of parity violation in colliders (such as the
Large Hadron Collider) without polarised beams.

Only one work (Ref. [5]) has previously tried to tackle this task in a
systematic way. To give a chastening example: when that work considered
collisions of the form $S_{ab\rightarrow jjj+X}$ (three important jets, and
perhaps some other stuff $X$ in the final state, resulting from the collision
of an ‘$a$’ with a ‘$b$’) it was found that 19-dimensional continuous ideal
vector parity (to use the nomenclature of [3]) would, in principle, detect any
parity-violating signature. However: (i) over 60 pages of algebra were
necessary to generate that result, and (ii) there seems little prospect of
generalising that result to four final-state particles (due to factorial
increases in computation complexity) let alone to $n$ final-state particles.

Why do we mention these limitations of [5]?

The reason is that the present work can be considered to be doing a very
similar job to [5] but for events with ‘up to’ 784 = 28x28 jets in the final
state, yet in a way that scales much better.333The intensity of each pixel
can, in principle, represent a very narrow jet, and so an image with (say)
only three black pixels could represent a three-jet event. Moreover, the
permutation symmetry modelled between jets in [5] carries over to our images,
since the pixels to not somehow list jets in a particular order. The one major
difference between [5] and the present work is, however, that [5] provides
guarantees of ideality (see definition of this term in [3]) which the present
work never can.

##  6 Conclusion

It was demonstrated that functions with appropriate symmetries (in our case
(i) parity oddness and (ii) invariance with respect to transformations
encoding isotropy) could demonstrate an underlying symmetry (in our case
parity symmetry) was violated by being able to learn to label data with
positive values.

The worked example presented here was intended to be pedagogical in nature and
so considered only a simple symmetry (parity) in an ideal detector. The issues
which were ignored and would concern real detectors (notably acceptance and
inefficiency) are considered in detail in a sister paper [7].

##  7 Acknowledgements

The authors greatly acknowledge discussions with Ben Nachman related to the
present work and to [1], and with Shikma Bressler related to the present work
and to [8]. Useful discussions with Radha Mastandrea, Daniel Noel and Vidhi
Lalchand are all very much appreciated, as has been the patience of Tom Gillam
in assisting with numerous questions about pytorch.

## References

  * [1] Jack Collins, Kiel Howe and Benjamin Nachman  “Anomaly Detection for Resonant New Physics with Machine Learning”  In _Phys. Rev. Lett._ 121 American Physical Society, 2018, pp. 241803  DOI: 10.1103/PhysRevLett.121.241803
  * [2] Yann LeCun and Corinna Cortes  “MNIST handwritten digit database”, http://yann.lecun.com/exdb/mnist/, 2010  URL: http://yann.lecun.com/exdb/mnist/
  * [3] Christopher G. Lester  “Chiral measurements”, 2021  arXiv:2111.00623 [hep-ph]
  * [4] Christopher G. Lester  “Code to support this paper”  Zenodo, 2021  URL: https://github.com/kesterlester/ColabForParityPaper3
  * [5] Christopher G. Lester, Ward Haddadin and Ben Gripaios  “Lorentz and permutation invariants of particles III: constraining non-standard sources of parity violation”, 2020  arXiv:2008.05206 [hep-ph]
  * [6] Adam Paszke et al.  “PyTorch: An Imperative Style, High-Performance Deep Learning Library”  In _Advances in Neural Information Processing Systems 32_ Curran Associates, Inc., 2019, pp. 8024–8035  URL: http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * [7] Rupert Tombs and Christopher G. Lester  “A method to challenge symmetries in data with self-supervised learning”, 2021  arXiv:2111.05442 [hep-ph]
  * [8] Sergey Volkovich, Federico De Vito Halevy and Shikma Bressler  “The Data-Directed Paradigm for BSM searches”, 2021  arXiv:2107.11573 [hep-ex]


#  Through the lens of Sgr A∗: identifying and resolving strongly lensed
Continuous Gravitational Waves beyond the Einstein radius

Stefano Savastano<EMAIL_ADDRESS>Max Planck Institute for
Gravitational Physics (Albert Einstein Institute)
Am Mühlenberg 1, D-14476 Potsdam-Golm, Germany  Filippo Vernizzi
<EMAIL_ADDRESS>Institut de Physique Théorique, Université Paris
Saclay CEA, CNRS, 91191 Gif-sur-Yvette, France  Miguel Zumalacárregui
<EMAIL_ADDRESS>Max Planck Institute for Gravitational
Physics (Albert Einstein Institute)
Am Mühlenberg 1, D-14476 Potsdam-Golm, Germany

###### Abstract

Lensed gravitational waves will offer new means to probe the distribution of
matter in the universe, complementary to electromagnetic signals. Lensed
continuous gravitational waves provide new challenges and opportunities beyond
those of transient compact binary coalescence. Here we consider continuous
gravitational waves emitted by isolated neutron stars and lensed by Sgr A∗,
the supermassive black hole at the center of our galaxy, a system observable
by the next generation of gravitational wave detectors. We analyze the
signatures of this system in detail, addressing parameter estimation and model
selection. Future detectors can distinguish lensed continuous waves and
measure their parameters with precision $\sim 1-10\%$ for sources within $2-4$
Einstein radii of Sgr A∗, depending on the source distance, thanks to the
relative motion of the observer-lens-source system.

The chances of observing strongly-lensed neutron stars increase by one order
of magnitude relative to previous estimates, thanks to the possibility of
detecting lensed systems at several Einstein radii. Multiple images can be
resolved with an angular accuracy $\sim 10$mas, comparable to the best optical
telescopes. Image localization probes deviations from axial symmetry and the
existence of companions to Sgr A∗ in regions complementary stellar orbits and
black hole imaging. Our methods and many of our results extend to other lenses
(e.g. galactic substructure) and sources (e.g. long-lived inspiralling
binaries), rendering lensed continuous gravitational waves into versatile
probes of astrophysics and fundamental physics.

Introduction. Gravitational lensing has become a powerful tool for
astrophysics and cosmology, for example in the search for dark objects,
exploration of the universe’s large-scale structure and measurement of
cosmological parameters Bartelmann (2010); Oguri (2019); Liao et al. (2022).

With the recent rise of gravitational wave (GW) astronomy, the lensing of GWs
emitted by coalescing binary black holes (BBHs) and neutron stars (NSs) has
become the subject of intense research. The coherence, low frequency and
frequency evolution of these sources enables the observation of diffraction
Takahashi and Nakamura (2003); Dai et al. (2018); Çalışkan et al. (2022);
Tambalo et al. (2022a, b) and phase Dai and Venumadhav (2017); Ezquiaga et al.
(2021); Vijaykumar et al. (2022) effects that are challenging to observe with
electromagnetic (EM) waves. While no detection of lensed GWs has yet been made
Hannuksela et al. (2019); Dai et al. (2020); Abbott et al. (2021); Basak et
al. (2022), the increasing rate of GW observations offers a promising future
Ng et al. (2018); Xu et al. (2022).

Besides transient signals produced by explosive binary coalescences, detectors
can observe long-lived signals with a slow frequency evolution. Sources of
these quasi-monochromatic signals fall into two main categories: 1) stellar-
mass binaries well before coalescence, which will be detectable by LISA Amaro-
Seoane et al. (2017); Sesana (2016); Wagg et al. (2022), and exotic binaries,
such as sub-solar primordial BBHs, which can be searched with ground-based
detectors Kalogera et al. (2021); Miller et al. (2021); 2) rapidly rotating
non-axisymmetric neutron stars with a quadrupolar deformation, producing
continuous GWs (CWs), observable by ground-based detectors. These signals have
been searched for in the LIGO-Virgo-Kagra (LVK) data Abbott et al. (2022a);
Steltner et al. (2023); Dergachev and Papa (2022); see Refs. Riles (2017);
Sieniawska and Bejger (2019); Piccinni (2022); Riles (2022) for reviews.

The lensing of CWs presents interesting and distinct characteristics.
Specifically, due to the coherent nature of CWs, a gravitational lens can act
as a diffractive barrier, resulting in interference fringes on the detected
signal when the observer-lens-source system undergoes relative transverse
motion Liao et al. (2019); Suvorov (2022); Moylan et al. (2007); Jung and Kim
(2023); Takahashi et al. (2023). Additionally, CWs allow for very precise sky
localization, impossible with binary coalescence observations. Despite the
promising prospects for distinguishing between lensed and unlensed CWs, the
ability of observations to accurately infer the lens parameters remains an
open question.

While no clear detection of CWs has been reported yet, the sheer number of NSs
in the Milky Way (MW) galaxy, estimated in $\sim 10^{9}$ from population
synthesis studies Narayan (1987), suggests that detection of lensed signals is
plausible for future observations. Motivated by such a prospect, we study the
phenomenology of strongly-lensed CWs, the inference of lens parameters and the
potential for detection by next-generation GW observatories. For concreteness,
we focus on monochromatic, isolated, rotating NS moving at constant velocity
and lensed by Sgr A∗, the supermassive black hole (BH) at the center of the MW
Ghez et al. (2000); Genzel et al. (2010); Akiyama et al. (2022). Third-
generation (3G) interferometers will be able to observe up to $\sim 6$
strongly lensed signals of rapidly spinning NSs within the Einstein cone of
Sgr A∗ Basak et al. (2023), the precise number depending on the NS properties
and their spatial distribution in the MW Reed et al. (2021). While the Sgr A∗
scenario is compelling, our methods and many of our conclusions extend to
other systems involving lensed CWs.

Strong lensing imprints on CWs. Let us summarize the effect of gravitational
lensing on a CW. We work in the geometric optics (GO) limit and denote the
observed and true angular position of the source by $\bm{\theta}_{l}$ and
$\bm{\theta}_{s}$ and the observer-lens, observer-source and lens-source
angular-diameter distances by $r_{l}$, $r_{s}$ and $r_{ls}$; see Fig. 1.

Figure 1: Strong lensing of a quasi-monochromatic source. Top: Gravitational
lensing by Sgr A∗ forms two images of a rotating NS. Here,
$\bm{\eta}=\bm{\theta}_{s}r_{s}$ is the physical impact parameter in the
source plane, related to $y$ by $|\bm{\eta}|=y\theta_{E}r_{s}$. Middle: static
case, $h_{+},h_{-}$ are magnified and delayed. Their interference at the
detector is again a monochromatic waveform. Bottom: moving case, the images
are additionally slightly red/blue-shifted. Their interference produces a
modulated signal. (here exaggerated for illustrative purposes).

Lensed CWs can be determined via 1) source frequency evolution 2) modulation
due to transverse motion, or 3) spatially resolved images. Typically, rapidly
rotating NSs have small period variations $\dot{P}\lesssim 10^{-18}$s/s
Harding (2013), precluding frequency evolution from revealing lensed systems
(see App. A). Below we discuss the transverse motion modulation, leaving the
discussion of spatially resolved images for later.

In the presence of transverse motion, the time delay of an image $j$, $t_{j}$,
can be expanded at linear order in $t$ around a reference time $t_{0}$, hence

$t_{j}(t)=t_{j}(t_{0})+\bm{\alpha}_{j}(t_{0})\cdot\bm{v}(t-t_{0}),$ (1)

where
$\bm{\alpha}_{j}\equiv\frac{r_{l}}{r_{ls}}(\bm{\theta}_{j}-\bm{\theta}_{s})$
is the deflection angle for the $j$-th image, and $\bm{v}$ is the projected
transverse velocity, given by Kayser et al. (1986)

$\bm{v}=\bm{v}_{l}-\frac{r_{l}}{r_{s}}\bm{v}_{s}-\frac{r_{ls}}{r_{s}}\bm{v}_{o}\;,$
(2)

where $\bm{v}_{o}$, $\bm{v}_{l}$ and $\bm{v}_{s}$ are the transverse velocity
of the observer, lens and source in their respective planes. In general, the
magnification and deflection angle of each image also acquire a time
dependency from the transverse motion, but the effect is negligible for heavy
lenses not closely aligned with the sources, as considered here.

The lensed waveform is proportional to the unlensed one,
$\tilde{h}_{l}=F\tilde{h}_{0}$, where the amplification factor reads

$F\equiv\frac{\tilde{h}_{l}(f)}{\tilde{h}_{0}(f)}=\sum_{I}\sqrt{|\mu_{j}|}\exp\big{[}2\pi
i\left(ft_{j}+n_{j}\right)\big{]}\,.$ (3)

Here $j$ labels the GO images, $\mu_{j}$, $n_{j}$ are the magnification and
Morse phase of the $j$-th image (see App.A for details). Moreover, it is
convenient to define the dimensionless impact parameter,
$y=|\bm{\theta}_{s}|/\theta_{E}$, where
${\theta_{E}\equiv\sqrt{\frac{4GM_{l}r_{ls}}{r_{l}r_{s}}}}$ is the Einstein
angle.

The time-varying time delay in Eq.(1) induces an effective blue/red-shift on
each image’s frequency by a factor
$z_{j}=\frac{d}{dt}{\ln(t_{j})}=\bm{v}\cdot\bm{\alpha}_{j}$ . When two images
$i$ and $j$ interfere in the detector, the strain exhibits an amplitude and
phase modulation arising from the difference $z_{i}-z_{j}$. Given a lensing
signal with lens mass $M_{l}$, effective lens velocity $v\equiv|\bm{v}|$, and
impact parameter $y_{0}\equiv y(t_{0})$, potential information on these
parameters is contained in three features: 1) the modulation period,

$T_{ij}=\frac{1}{f_{0}v\alpha_{ij}}=2.4\,{\rm d}\left(\frac{1{\rm
kHz}}{f_{0}}\right)\left(\frac{10^{-3}}{v}\right)\left(\frac{1^{\prime\prime}}{\Delta\alpha_{ij}}\right)\,,$
(4)

where $\Delta\alpha_{ij}=|\bm{\alpha}_{i}-\bm{\alpha}_{j}|$; 2) the modulation
amplitude $\sqrt{\mu_{j}/\mu_{i}}$, which depends purely on $y(t)$; and 3) the
modulation phase, given by $\Delta t_{ij}\propto M_{l}$. Additional effects on
CWs (detector motion, source’s frequency evolution or orbital motion, two
overlapping sources), can be corrected for, or distinguished from, the
modulation signature (see App. C).

Figure 2: MC results. For the fiducial setup, two cases are considered:
uniform prior on the lens mass (solid green line); fixed lens mass (dashed
green line). For variable mass, the accuracy remains for very small relative
velocity, $v=10^{-5}$ (dotted). Top three panels: 68% C.L. marginalized error
on $M_{l}$ (top), $y_{0}$ (middle) and $v$ (bottom), as a function of the
initial impact parameter $y_{0}$. Bottom panel: Bayes’ ratio of lensed vs
unlensed hypothesis against the initial impact parameter for multiple distance
ratios, $r_{ls}/r_{l}$. $\log{\cal B}_{LU}>10$ corresponds to _very strong
evidence_ for the lensing hypothesis. The purple line shows the case of
advanced LIGO for $r_{ls}/r_{l}=1$.  Figure 3: 1-d and 2-d marginalized
1-$\sigma$ C.L. posteriors of the lens parameters, for different initial
impact parameters (fiducial setup).

Lens parameter estimation. To study the prospect of measuring $M_{l}$, $y_{0}$
and $v$, we follow a Bayesian approach and compute their marginalized
posteriors through numerical Monte Carlo (MC) sampling. Parameter inference
and sampling methods are discussed in App. D and E. We assume a single
Einstein Telescope detector and in-band signal duration $T_{\rm obs}=10\,{\rm
yr}$. Parameter estimation will improve using a global network of 3G detectors
Kalogera et al. (2021), with increased SNR and sky coverage. As a fiducial
setup, we consider a point-like lens with the measured mass of Sgr A∗,
$M_{l}=4.154\cdot 10^{6}\,M_{\odot}$ and R. Abuter et al. (2019) and lens
velocity $v=10^{-3}$, at a distance $r_{l}=8.178\,{\rm kpc}$ and R. Abuter et
al. (2019). A point lens produces a brighter and a fainter image, indicated
respectively as $+$ and $-$. The source is a spinning NS emitting CWs at a
frequency $f_{0}=800\,{\rm Hz}$ and relative distance $r_{ls}=r_{l}$. We will
discuss how the results change away from these fiducial values.

The lens and source parameters can be accurately measured even at large impact
parameter, despite degenerate posteriors. The three top panels in Fig. 2 show
68% C.L. limits on the lensing parameters from the 1-dimensional marginalized
posteriors, obtained from the MC samples, as a function of the initial impact
parameter $y_{0}$ (fiducial setup). Since the mass of Sgr A∗ is known with
sub-percent accuracy and R. Abuter et al. (2019), we run two analyses: one
with the lens mass $M_{l}$ treated as a free parameter, and another with
$M_{l}$ fixed to its known value. In both cases, the lens parameters can be
extracted from the signal beyond $y_{0}=1$. For example, for $y_{0}=3$ all
parameters are constrained with a relative error smaller than $50\%$ at
$1\sigma$ C.L. At large enough $y_{0}$, the sampled contours become consistent
with $\mu_{2}=0$ (i.e., no second image detected) and lens parameters cannot
be constrained. Note that, for strictly monochromatic sources, $M_{l}$ and $v$
only enter the GW signal phase and can be constrained only up to a periodic
factor, see Eq. (3). Figure 3 shows the 1- and 2-d 95% C.L. marginalized lens
parameter posteriors from the MC. For large impact parameter, particularly
noticeable in the $y_{0}=3$ contour, the posteriors deviate from Gaussian
behaviour. For additional information and details, see App. E.

To assess our capacity to differentiate between the lensed and null (i.e.,
unlensed) hypotheses, we compute the Bayes ratio between the two models,
denoted as $\mathcal{B}_{LU}$. The bottom panel in Fig. 2 shows its variation
as a function of the initial impact parameter for different distance ratios
$r_{ls}/r_{l}$ (fiducial setup). Very strong evidence for the lensed signal
(i.e. $\log\mathcal{B}_{LU}>10$) can be robustly established for $y_{0}\leq
y_{\rm max}\simeq 5$ (2.2) for $r_{ls}/r_{l}=0.1$ ($2$). A similar analysis
for advanced LIGO shows a more modest gain, with $y_{\rm max}\simeq 1.5$, for
$r_{ls}/r_{l}=1$. Note that the evidence for the unlensed and static lens
hypotheses are identical, as they are both degenerate for a monochromatic
signal.

We explored different setups varying the fiducial parameters. In the Gaussian
limit, the covariance matrix scales as the inverse of the ${\rm
SNR}\sim\sqrt{T_{\rm obs}}\mathcal{A}/r_{s}$ (cf. Eq. (27) in App. D). For
fixed SNR, changing the frequency of the source influences our results
marginally, as long as the modulation period in Eq. (4) is smaller than or
comparable to the observation time. Similarly, while parameter estimation
requires a non-zero $v$, the error and the evidence for the lensed signal is
rather insensitive to its exact value as long as $v\gtrsim 10^{-5}\approx
3\,{\rm km/s}$, which is high enough to observe the modulation period. Thus,
small velocities are unlikely to prevent the observation of lensed NS.
Moreover, the sensitivity to very small velocities suggests that CWs could
also measure lens accelerations. The parameter with the largest influence on
the precision is the distance ratio $r_{ls}/r_{l}$, with smaller values
leading to better parameter estimation and larger evidence.

Resolving individual images. CWs also enable accurate sky localization. In
order to identify a CW in a blind search, it is necessary to “undo” the
detector’s motion by redefining the time variable:
$t=t^{\prime}+\hat{n}\cdot\bm{x}(t^{\prime})+\Delta T$ Brady et al. (1998).
Here $t^{\prime}$ is the observation time, $\hat{n}$ is the source’s
direction, $\bm{x}$ is the detector’s position relative to the Solar System
barycenter and $\Delta T$ is the relativistic time delay. Both Earth’s
rotation and orbital motion contribute to $\bm{x}$, but the latter becomes
dominant after an integration period of a few days Maggiore (2007). Once a
signal is identified, analysis exploiting the complete coherence of the signal
can accurately determine an image’s position, depending on the source
location.

We determine the accuracy in the image position using Eq. (7.151) of Ref.
Maggiore (2007). This generally depends on the relative inclination between
the Earth’s orbital plane and the source direction, so the resulting skymap is
not isotropic. We have computed the size and shape of the skymap for each
image, which is shown in the insets of Fig. 4 (details on the derivation in
App. C). The heart-shaped maps are tilted with respect to the optical axis by
an angle $\phi_{\rm sky}$ that encodes the source-galactic plane relative
inclination, and their size scales as $\propto 1/(\sqrt{\mu_{i}}{\rm SNR})$. A
simple analytical estimate can be derived considering the Earth to be on a
circular orbit around the Sun, with radius $1\,{\rm AU}$, and the source
direction on the orbital plane. Then, an observation time of $T_{\rm
obs}>1\,{\rm yr}$ corresponds to an angular accuracy of
${\sigma_{\theta}\simeq 21\,{\rm mas}\left(60/{\rm SNR}\right)\left(1\,{\rm
kHz}/f_{0}\right)}$. This is comparable to the mean angular accuracy of the
anisotropic setup. In general, we expect to be able to individually resolve
the images for lenses with mass $M_{l}>1.2\times
10^{3}M_{\odot}\left(r_{l}/8.18\,{\rm
kpc}\right)\left(r_{ls}/r_{s}\right)\left(60/{\rm SNR}\right)\left(1\,{\rm
kHz}/f_{0}\right)$, at $y=1$. [MZ: check]

Figure 4: Resolving individual images. The insets show the 68,95% C.L.
localization regions (filled shapes), the unpertubed images (crosses) and the
projected lens-image axes (filled bands). A $\sim 4\cdot 10^{4}M_{\odot}$
perturber at $\bm{x}_{m}=(1.5,-0.5)$ displaces the positive-parity image by
$\gtrsim\sim 3\sigma$ away from the optical axis. Einstein radii are shown in
gray. The projected orbit of the S2 star Abuter et al. (2020) and the image of
Sgr A∗ Akiyama et al. (2022) are shown for comparison.

Probing Sgr A∗ companions. Let us consider the prospect of searching for
additional objects via their effect on the image positions. A smoking gun for
such objects is a misalignment between Sgr A∗ and the two images, as this
would require breaking the axial symmetry. Let us consider a point-lens
perturber with mass $M_{\rm pert}=m\cdot M_{\rm SgrA^{*}},\,(m\ll 1)$ at
$\bm{\theta}_{m}=(\theta_{m}^{\parallel},\theta_{m}^{\perp})$ and focus on the
off-axis displacement of the images
${\Delta\theta^{\perp}_{i}\approx\frac{-m\,\theta_{E}}{1-(\theta_{i}/\theta_{E})^{-2}}\frac{\theta_{m}^{\perp}}{|\bm{\theta}_{i}-\bm{\theta}_{m}|^{2}}}$
(assumed small), where $\bm{\theta}_{i}$ is the position of the unperturbed
image. We will ignore other effects, such as time delays and the introduction
of a faint image with magnification $\propto m^{2}$ (see App. B for details).

Following the sky localization accuracy discussed above, each image defines an
axis through the known lens position. Figure 4 illustrates this setup with an
example perturber leading to a $\sim 3\sigma$ measurable displacement. The
blue/red shaded regions represent the 1-2$\sigma$ confidence bands on the
lens-image axis. In practice, the axis of the positive parity image $x_{+}$
determined more precisely due to better intrinsic localization and larger
level arm (i.e. both $|\mu_{+}/\mu_{-}|,|\theta_{+}/\theta_{-}|>0$).
Therefore, we consider $\bm{\theta}_{+}$ to fix the optical axis and its
uncertainty, while $\bm{\theta}_{-}$ determines the sensitivity to off-axis
image displacement.

Figure 5: Minimum perturber mass producing a 2-$\sigma$ lens-image
misalignment (cf. Fig. 4), depending on its position in the lens plane. The
source is our fiducial setup at $y=0.5,0.1$ (upper/lower), assuming the median
error in sky localization. Lines indicate $10\%,1\%$ and $0.1\%$ of the mass
of Sgr A∗. Orange lines centered around Sgr A∗ show the corresponding 95% C.L.
limits from the Schwarzschild precession of the star S2 Abuter et al. (2020),
assuming $r_{ls}=r_{l}$.

Figure 5 shows the 95% C.L. minimum detectable mass via the offset between Sgr
A∗ and the two images of a rotating NS at $y=0.5,0.1$ (with fiducial
parameters), as a function of the perturber’s position. Precision degrades at
large $y$, both because the negative parity image becomes faint (poorly
localized) and close to Sgr A∗. For $y\lesssim 1$ the area over which a
perturber can be detected scales as $\sim\frac{m}{y^{3}}$, increasing greatly
for closely aligned systems. This analysis is conservative, as combining sky
localization with time-domain information (cf. Fig. 2) will improve
sensitivity. This method is complementary to other probes of the galactic
center: analogue constraints form pericenter passage of the S2 star Abuter et
al. (2020) probe the region near the star’s orbit, rather than the GO images.
Thus, lensed NSs are sensitive to different regions in the galactic center. In
addition, gravitational lensing is sensitive to perturbers at intermediate
distances (i.e. between the observer and lens or between the lens and source).

Discussion and prospects. Lensed CWs offer strong complementarity to lensed EM
signals and GWs from compact binary coalescences, bringing new challenges and
opportunities. We have established the prospect of reconstructing the lens
parameters and further probing the lens by individually resolving the images.
We have focused on rotating NSs lensed by Sgr A∗, but our conclusions extend
to other systems.

The prospect of identifying lens systems at several Einstein radii leads to
enhanced detection probabilities. The number of potentially lensed sources
scales with the square of $y_{\rm max}$, i.e. the maximum impact parameter
within which a detection can be established with strong evidence. In Ref.
Basak et al. (2023), it was shown that 3G detectors can detect up to $\sim 6$
NSs within the Einstein radius of Sgr A∗, corresponding to $y_{\rm max}=1$.
Even when applying a stringent threshold for strong evidence (cf. Fig. 2), we
find that $y_{\rm max}$ can be as large as $\sim 3$. This increases in the
number of detectable sources by an order of magnitude compared to the estimate
based on $y_{\rm max}=1$. (See the App. F for more details). While the
ultimate prospects depend on unknown astrophysics, a targeted search of NSs
lensed by Sgr A∗ is warranted.

Individual images can be resolved with $\sim 10$mas accuracy (Fig. 4),
comparable to the best optical telescopes. As an application, we showed how
lens-image misalignment provides a smoking gun for additional structure (Fig.
5), probing companion objects in regions complementary to stellar orbits Naoz
et al. (2020); Abuter et al. (2022); Will et al. (2023) and BH imaging Akiyama
et al. (2022). Future work will address additional signatures of these objects
and the benefits of combining timing and sky-localization information.

CWs lensed by Sgr A∗ provide a novel probe of the galactic center. Detection
of lensed EM radiation from objects closely aligned with Sgr A∗ is challenging
due to abundant stars, gas and dust in the central region of the galaxy Bozza
and Mancini (2012). In contrast, lensed CWs are negligibly absorbed, providing
a pristine view of the region near our supermassive BH. CWs will complement
other approaches to probe the matter distribution around Sgr A*, test dark
matter scenarios Gondolo and Silk (1999); Hui et al. (2017) and GW propagation
Ezquiaga and Zumalacárregui (2020); Goyal et al. (2023) in an extreme
environment.

This study is a first step towards understanding parameter reconstruction in
lensed CWs. We have focused on NSs lensed by Sgr A∗, but our results apply to
other lenses and sources. Lensing of CWs from NSs in all-sky searches will
probe compact objects and galactic substructure, enabling novel tests of
intermediate-mass black holes and the dark matter distribution. Our treatment
can be extended to other long-lived sources, such as inspiraling compact
binaries observable by planned or proposed space detectors Amaro-Seoane et al.
(2017); Gong et al. (2021); Sedda et al. (2020); Baibhav et al. (2021).
Expanding beyond our simplifying assumptions—i.e. constant magnifications,
geometric optics and point-like lenses— will unveil the full potential of
lensed CWs and enable novel tests of astrophysics and fundamental physics.

###### Acknowledgements.

We are very grateful to Soummyadip Basak, Jing Ming, Ornella Juliana Piccinni,
Giovanni Tambalo and Hector Villarrubia-Rojo for discussions, Lorenzo Speri,
Jonathan Gair and Ollie Burke for help with aspects of the parameter
estimation, and the anonymous referees for valuable and constructive
criticism.

##  Appendix A Gravitational lensing

In this section, we will introduce some concepts of gravitational lensing that
will be beneficial for the main text. Using spherical coordinates centered at
the observer, $\vec{r}=(r,\theta,\varphi)$, we define the two-dimensional
vector $\bm{\theta}=\theta(\cos\varphi,\sin\varphi)$. We denote by
$\bm{\theta}_{l}$ and $\bm{\theta}_{s}$ the observed and true angular
positions of the source, respectively. Additionally, we use $r_{l}$, $r_{s}$,
and $r_{ls}$ to represent the observer-lens, observer-source, and lens-source
angular-diameter distances, respectively.

In wave optics, the frequency-domain amplification factor due to the presence
of a lens at rest, $F\equiv\tilde{h}_{l}/\tilde{h}_{0}$, is given in the form
of a diffraction integral Nakamura and Deguchi (1999); Takahashi and Nakamura
(2003),

$F(\omega,\bm{\theta}_{s})=\frac{\omega}{2\pi i}\frac{r_{l}r_{s}}{r_{ls}}\int
d^{2}\bm{\theta}\exp[i\omega t(\bm{\theta},\bm{\theta}_{s})]\;.$ (5)

Here we have defined

$t(\bm{\theta},\bm{\theta}_{s})=\frac{r_{l}r_{s}}{2r_{ls}}|\bm{\theta}-\bm{\theta}_{s}|^{2}-\hat{\psi}(\bm{\theta})\,,$
(6)

with the lensing potential given as $\hat{\psi}(\bm{\theta})\equiv\int
dr\,U(r,\bm{\theta})$. In the limit of geometric optics (GO), which applies
when $\omega t(\bm{\theta},\bm{\theta}_{s})\gg 1$, the diffraction integral is
dominated by the stationary points of $t(\bm{\theta},\bm{\theta}_{s})$ and
each image position $\bm{\theta}_{j}$ can be obtained by solving the lens
equation Schneider et al. (1992),

$\nabla_{\bm{\theta}}t=0\;.$ (7)

Moreover, for each image $j$,
$\bm{\alpha}_{j}\equiv\frac{r_{l}}{r_{ls}}(\bm{\theta}_{j}-\bm{\theta}_{s})$
denotes its deflection angle,
$\mu_{j}(\bm{\theta}_{s})\equiv\det(\partial\bm{\theta}_{s}/\partial\bm{\theta}_{j})^{-1}$
its magnification, and $t_{j}(\bm{\theta}_{s})\equiv
t(\bm{\theta}_{j},\bm{\theta}_{s})$ its time of arrival, which is a solution
of the above equation. In this limit, the amplification factor simplifies to a
discrete sum over images, as shown in Eq. (3) in the main text.

Figure 6: Amplification factor of a point lens as a function of the
dimensionless frequency $w=8\pi GM_{l}\,f$ and impact parameter
$y=|\bm{\theta}_{s}|/\theta_{E}$. Dashed grey contours correspond to the
transition between wave and geometric optics regimes (less than $10\%$
relative difference). Diffraction and wave effects can be observed for
chirping binaries (varying $w$ and keeping $y$ constant) or moving lens
systems (varying $y$ and keeping $w$ constant) as explored in this work. Our
fiducial system lies in geometric optics with $w=8\pi GM_{\rm
SgrA^{*}}f_{0}\simeq 10^{3}-10^{4}$.

For convenience, we convert angles to dimensionless coordinates by normalizing
by the Einstein angle
${\theta_{E}\equiv\sqrt{\frac{4GM_{l}r_{ls}}{r_{l}r_{s}}}}$, i.e.
$\bm{x}\equiv{\bm{\theta}}/{\theta_{E}}$. In particular, we indicate the
$j$-image’s position as $\bm{x}_{j}\equiv{\bm{\theta}_{j}}/{\theta_{E}}$ and
the dimensionless impact parameter as
$\bm{y}\equiv\bm{\theta}_{s}/\theta_{E}$. We choose the horizontal axis of our
coordinate system to coincide with the optical axis, i.e. $\bm{y}=(y,0)$.
Moreover, we define the dimensionless frequency, $w\equiv 4GM_{l}\,\omega$,
where $M_{l}$ is the lens mass.

Figure 6 shows $|F|$ as a function of $w$ and $y$ for a point lens and
highlights the region where geometric optics applies. A point-like lens always
splits the source into two images, identified by $+$ and $-$, with Schneider
et al. (1992)

$\displaystyle\frac{\mu_{-}}{\mu_{+}}$
$\displaystyle=\frac{2-y(\sqrt{y^{2}+4}-y)}{2+y(\sqrt{y^{2}+4}+y)}\;,$
$\displaystyle\Delta t_{+-}$
$\displaystyle={2GM_{l}}\left[y\sqrt{y^{2}+4}+2\log\left(\frac{\sqrt{y^{2}+4}+y}{\sqrt{y^{2}+4}-y}\right)\right]\;,$
$\displaystyle x_{\pm}$
$\displaystyle=\frac{1}{2}\left(y\pm\sqrt{y^{2}+4}\right)$
$\displaystyle\alpha_{\pm}$
$\displaystyle=\sqrt{\frac{4GM_{l}r_{s}}{r_{l}r_{ls}}}x_{\pm}\,,$

where we denote by $\Delta t_{+-}$ the time delay elapsed between the two
images. For example, at $y=1$ and $r_{ls}/r_{l}=1$, lensing by Sgr A∗ induces
a $\Delta t_{+-}=163.9\,{\rm s}$ and $\Delta\alpha_{+-}=6.3^{\prime\prime}$.
We will not consider the effects of external convergence and shear An and
Evans (2006), as they will produce negligible corrections for galactic lenses
at $y\sim\mathcal{O}(1)$.

Projected transverse motion results in a time variation of the lensing
functions. In fact, this is equivalent to considering a source position that
varies with time (as per Eq. (2) in the main manuscript). Specifically, we
have
${y(t)=\sqrt{y_{0}^{2}+\left(\frac{v(t-t_{0})}{r_{l}\theta_{E}}\right)^{2}}}$,
where $y_{0}$ is the impact parameter at $t=t_{0}$.

##  Appendix B Perturbations by a point lens

Let us introduce the dimensionless lensing potential, the so-called Fermat
potential,

$T(\bm{x})=\frac{t(\bm{x})}{4GM_{l}}\;,$ (8)

and consider the effect of a point lens perturber with mass $M_{\rm
pert}=m\cdot M_{\rm SgrA^{*}},$ at
$\bm{x}_{m}=(x_{m}^{\parallel},x_{m}^{\perp})$ (here $\bm{x}$ is defined by
normalizing to the Einstein angle of the unperturbed lensed, as defined
above). The total lensing potential is

$T(\bm{x})=T^{(0)}(\bm{x})+m\log(|\bm{x}-\bm{x}_{m}|)\,,$ (9)

where $T^{(0)}(\bm{x})$ is the unperturbed Fermat potential.

We will assume $m\ll 1$ sufficiently small for a perturbative treatment and
$y$ large enough for the source to be outside the caustic network of the
binary lens. Under our assumptions, the system forms an additional image at

$\bm{x}_{3}\approx\bm{x}_{m}-\bm{M}^{-1}\cdot\bm{\Delta}_{ym}\;,$ (10)

where the vector $\bm{\Delta}_{ym}$ is defined as

$\bm{\Delta}_{ym}=\bm{y}-\bm{x}_{m}\left(1-x_{m}^{-2}\right)\;,$ (11)

and the elements of the matrix $\bm{M}$ are given by

$M_{ij}=\delta_{ij}\bigg{(}\frac{\Delta^{2}_{ym}}{m}+\frac{1}{x_{m}^{2}}\bigg{)}-\frac{2{x}_{m}^{i}{x}_{m}^{j}}{x_{m}^{4}}\;.$
(12)

The third image will undergo a frequency shift and induce an additional
amplitude modulation, similar to the interference effect described in the main
text. However, the amplitude of the third image is typically very small, as
$\mu_{3}\propto m^{2}/|\bm{y}-\bm{x}_{m}|^{4}$. Even for GWs,
$\sqrt{|\mu_{3}|}h_{0}$ will most likely be below the detection threshold.
Moreover, the existence of an additional image at arbitrary $y$ is a feature
of point lenses and other very compact matter distributions, while it is
absent in generic extended lenses (cf. Ref. Tambalo et al. (2022b), Sec IIIA).
We will thus focus on the perturber’s effect on the position and time delay
between the main images. (Given that $\delta\mu_{\pm}/\mu_{\pm}\propto
m/|\bm{y}-\bm{x}_{m}|^{2}$, the effect on the magnification is negligible in
most cases.)

The leading-order effect on the Fermat potential of each image is

$T_{\pm}-T_{\pm}^{(0)}=m\log(|\bm{x}_{\pm}-\bm{x}_{m}|)\,,$ (13)

where we have used Fermat’s principle on the unperturbed lens, i.e.
$\bm{\nabla}T_{\pm}^{(0)}=0$. (A more accurate result follows from evaluating
the full Fermat potential (9) on the perturbed image positions, Eqs. (16,
17).) Hence, the time delay between the images is proportional to

$\Delta T_{\pm}\approx
m\log\left(\frac{|\bm{x}_{+}-\bm{x}_{m}|}{|\bm{x}_{-}-\bm{x}_{m}|}\right)\,.$
(14)

Restoring the lens’ mass gives us the effect of a perturber on the relative
phase between the GO images:

$\Delta k=8\pi GM_{l}f_{0}\Delta T_{\pm}\sim 5.14\cdot
10^{5}\frac{m}{x_{m}}\left(\frac{f_{0}}{{\rm kHz}}\right)\,.$ (15)

Therefore, a perturber with mass $m$ can affect the phase at observable levels
up to $x_{m}\lesssim 8\cdot 10^{4}m/\sigma_{k}$, where $\sigma_{k}\sim 0.3$
(see analysis) is the sensitivity to the phase. The expected number of such
objects scales as the projected density times $x_{m}^{2}$ and it’s likely to
be significant: a perturber with $m\sim 2\cdot 10^{-6}$ (corresponding to a
mass $\sim 10M_{\odot}$, for which $8\pi GM_{l}f\sim 1$, i.e. at the WO
diffraction limit for a $1$kHz source) can influence the signal beyond the
Einstein radius. Nevertheless, the actual sensitivity will be much degraded by
the degeneracies between the lens mass and motion parameters (velocity,
acceleration), cf. Fig. 3 in the main text. Therefore, the contribution of
light objects to the time delay will result in a slight shift of these
parameters posteriors, with relative magnitude $\sim\sum_{i}m_{i}\ll 1$.

Let us now turn to the effects of a perturber on the image positions. In the
absence of perturber, the positive/negative parity images lie in the optical
axis (i.e. the source-lens axis) at $\bm{x}_{\pm}=(x_{\pm},0)$. The perturber
displaces the main images by

$\displaystyle\Delta x^{\parallel}_{\pm}$
$\displaystyle\approx\frac{m}{1+x_{\pm}^{-2}}\frac{x_{\pm}-{x}_{m}^{\parallel}}{{(\bm{x}_{\pm}-\bm{x}_{m})}^{2}}\,,$
(16) $\displaystyle\Delta x^{\perp}_{\pm}$
$\displaystyle\approx\frac{-m}{1-x_{i}^{-2}}\frac{{x}_{m}^{\parallel}}{{(\bm{x}_{\pm}-\bm{x}_{m})}^{2}}\,,$
(17)

The image displacement along the optical axis (16) can be probed, but it is
degenerate with the main lens parameters ($y$ and $M_{l},r_{s}$ via $x_{E}$).
In contrast, the off-axis displacement (17) reveals additional structure.

A perturber will generically produce an off-axis displacement of the images,
$\Delta x_{\pm}^{\perp}\neq 0$. Because the optical axis is not known, one can
only measure a misalignment between the two images and the black hole. If the
sky localization uncertainty were isotropic, the uncertainty in the $\pm$
image displacement with respect to the optical axis is given by

$\sigma^{\perp}_{\pm}=\sqrt{\sigma_{\pm}^{2}+({x_{\mp}}/{x_{\pm}})\sigma_{\mp}^{2}}=\sigma_{\pm}\sqrt{1+\left({x_{\pm}}/{x_{\mp}}\right)^{3}}\,.$
(18)

The first equality follows from adding in quadrature the uncertainty on the
image and the projected uncertainty of the other image along the direction of
the lens (see Fig. 4 in main text). The second equality employs the scaling of
an image’s localization accuracy with the magnification,
$\sigma_{\pm}^{2}=\sigma_{\mp}^{2}\mu_{\mp}/\mu_{\pm}$, and the fact that for
a point lens $\mu_{-}/\mu_{+}=x_{-}/x_{+}$. Precision degrades at large $y$
because the negative parity image becomes faint (poorly localized) and close
to Sgr A∗, thus reducing the level arm.

To account for the anisotropic sky-localization uncertainty
$\bm{\sigma}_{\pm}$ (see discussion in App. C), we consider the projected
uncertainty in the image’s axis orientation (through the angle with respect to
Sgr A∗) as $\varphi_{\pm}(\varphi,\phi_{\rm
sky})=\arctan\left(\frac{(\bm{x}_{\pm}+\bm{\sigma}_{\pm}(\phi_{\rm
sky}))_{2}}{(\bm{x}_{\pm}+\bm{\sigma}_{\pm}(\phi_{\rm sky}))_{1}}\right)$.
Here, $\varphi$ is the angle that defines the orientation of a vector
$\bm{\theta}=(\theta_{1},\theta_{2})$ in the lens plane, and $\phi_{\rm sky}$
is the angle that set the orientation of the optical axis with respect to the
Earth’s orbital plane and of the shape of the anistoropic sky-map. Then the
uncertainty in the projection of the axes is given by
$\sigma_{\pm}^{\perp}\in\left(\min_{\varphi}(\phi_{\pm}),\max_{\varphi}(\phi_{\pm})\right)x_{\pm}$.
In practice the $+$ image is always better localized (both because
$\mu_{+}>|\mu_{-}|$ and $x_{+}>|x_{-}|$). Therefore, we can interpret the
uncertainty $\sigma_{+}^{\perp}$ as associated with the optical axis (i.e.
measurement of $\phi_{\rm sky}$). The (larger) uncertainty in the negative-
parity image $\sigma_{-}^{\perp}$ gives the sensitivity to deviations from
axial symmetry and the presence of perturbers. Because uncertainties depend
strongly on the fiducial value of $\phi_{\rm sky}$, we will quote the median
sensitivities. Note, however, that the sensitivity improves substantially for
certain configurations, e.g. $\phi_{\rm sky}\sim\pm\pi/2$ (this information is
available from the localization of the positive-parity image). Combining sky-
localization and time-domain analysis will further improve the sensitivity.

##  Appendix C CWs lensed by Sgr A∗

Isolated spinning NSs are expected to produce CWs in the
$10^{2}-10^{3}\,\textrm{Hz}$ band through a variety of different mechanisms
Glampedakis and Gualtieri (2018); Sieniawska and Bejger (2019); Piccinni
(2022), and could be observed by future ground-based detector campaigns. For
an estimate of the total number of neutron stars that can be probed with
current and future detectors, we refer to Reed et al. (2021); Pagliaro et al.
(2023). Neglecting the effect of the frequency evolution and the detector’s
motion and response, the NS signal is monochromatic:

$h_{0}(t,r_{s})=\frac{\cal A}{r_{s}}\,e^{-i\phi(t)}\;,\qquad\phi(t)=2\pi
f_{0}t+\phi_{0}\;,$ (19)

where ${\cal A}$ depends on the source’s orientation and $\phi_{0}$ is the
phase value at $t=0$. The dominant mode of a CW emitted by a NS has amplitude
${\rm\cal A}=4\pi^{2}G{I_{zz}f_{\mathrm{0}}^{2}}\epsilon\,$, where $I_{zz}$ is
the moment of inertia of a perpendicular biaxial rotor spinning with axis
$\hat{z}$ and the ellipticity parameter, $\epsilon=(I_{xx}-I_{yy})/I_{zz}$,
describes the degree of anisotropy around the NS rotation axis. To date, no
CWs have been detected but frequency-dependent upper limits on their amplitude
and on $\epsilon$ at a fixed distance were inferred Abbott et al. (2022a, b);
Dergachev and Papa (2022). We adopt $\epsilon=10^{-7}$ and $I_{zz}=$
$10^{38}{\rm\,kg\,m^{2}}$ Abbott et al. (2022a), compatible with data,
theoretical expectations Suvorov et al. (2016); Gittins et al. (2020) and the
lower bound from population-based studies, $\epsilon\gtrsim 10^{-9}$ Woan et
al. (2018).

Depending on their distribution, a fraction of NSs in our galaxy is expected
to lie close to the line of sight of Sgr A∗. Moreover, they are expected to
have sizeable projected transverse velocities (cf. Eq.(2) in the manuscript)
Sartore et al. (2010) due to either natal kicks ($\sim 400$ km/s) Hobbs et al.
(2005), or the Solar System’s motion in the galaxy ($\sim 200$ km/s) Sofue et
al. (2009); Mró z et al. (2019); Eilers et al. (2019). These objects are the
subject of targeted CWs searches Abbott et al. (2022b). Sgr A∗ can act as a
foreground lens for sources sufficiently aligned with our line of sight.

CWs searches I think the aim of this subsection should be to show that lensing
does not interfere with or spoil CW searches. Now, we will briefly summarize
the methodology of CW searches and their interplay with the occurrence of
lensed signals in the datastream. Our aim is to show that standard searches
are able to detect lensed events.

In blind searches the data is pre-processed assuming certain source properties
(location, frequency evolution) and then the monochromatic CW signals are
searched via a fast-Fourier transform Brady et al. (1998) (we will follow Ref.
Maggiore (2007), Ch. 7 Maggiore (2007)). In particular, the pre-processing is
necessary to clean the signal from other sources of distortion such as the
orbital motion and frequency evolution of the signal. This scheme is referred
to as _coherent search_. Under this type of search, strongly lensed CWs would
first appear as two signals very close in frequency and sky localization.

The large number of unknown parameters and the large observation time required
for pre-processing signals make coherent searches computationally infeasible.
Instead, searches are typically performed using a semi-coherent method
Piccinni et al. (2020), in which data streams are divided into segments of
duration $T_{\rm coh}$. These segments are processed as described below and
then combined incoherently, i.e. neglecting the relative phase between stacks.
Therefore, GW detectors are capable of observing NSs up to distances of
Maggiore (2007)

$r_{\rm hor}=35\,{\rm
kpc}\times\left(\frac{\epsilon}{10^{-7}}\right)\left(\frac{f_{0}}{800\,{\rm
Hz}}\right)^{2}\left(\frac{T_{\rm obs}}{10\,{\rm yr}}\right)^{1/2}\gamma\;,$
(20)

where we have defined an efficiency factor,
${\gamma\equiv\left(\frac{3750}{\mathcal{N}}\right)^{1/4}}\left(\frac{4}{{\rm
SNR}_{\rm thr}}\right)\left(\frac{4.8\,10^{-25}}{\sqrt{S_{n}(f_{0})\,{\rm
Hz}}}\right)$ using the reference numbers expected in next-generation
detectors Maggiore (2007). Here, ${\rm SNR}_{\rm thr}$ is the signal-to-noise
ratio (SNR) threshold of the search and ${\mathcal{N}={T_{\rm obs}}/{T_{\rm
coh}}}$ is the number of stacks in which the data is divided. The signal’s
amplitude has been sky-averaged over the solid angle and polarization angle.
We stress that once a CW detection from the Galactic center has been
confirmed, follow-up analyses can exploit the full coherence of the signal,
with no SNR loss. The outcomes of our study hinge on this premise.

The first step in obtaining a monochromatic waveform is a resampling of the
signal to account for the detector’s motion. This is done by a redefinition of
the time variable: $t=t^{\prime}+\hat{n}\cdot\vec{x}(t^{\prime})+\Delta T$
Brady et al. (1998). Here, $t^{\prime}$ is the observation time, $\hat{n}$ is
the source’s direction, $\vec{x}$ is the detector’s position relative to the
Solar System barycenter and $\Delta T$ is the relativistic time delay. The
antenna pattern (given by $\hat{n}$) has been factored out of ${\cal A}$. As
discussed in Ref. Maggiore (2007) (see Eq. 7.151), resampling requires
accuracy on the source’s position at the level of $\delta\theta\lesssim
22^{\prime}\left(1{\rm d}/{T_{\rm coh}}\right)^{2}\left({\rm
kHz}/{f_{0}}\right)$, with $1\,{\rm d}\lesssim T_{\rm coh}\lesssim
10^{2}\,{\rm d}$. This limit corresponds to all-sky NS searches, where a
relatively low $T_{\rm coh}$ allows computational efficiency Brady et al.
(1998). Since the separation between images is $|\alpha_{\pm}|\lesssim
10^{\prime\prime}$ data resampling is common to all images for $T_{\rm
coh}\lesssim 11$ d. For larger coherence time, they will show up in two
different sky bins. As discussed in the main manuscript, once a source is
found, follow-up analyses with larger $T_{\rm coh}$ allow for fine sky
localization, with an accuracy ${\sigma_{\theta}\simeq
0.0021^{\prime\prime}\left(60/{\rm SNR}\right)\left(1\,{\rm
kHz}/f_{0}\right)}$, for $T_{\rm obs}>1\,{\rm yr}$. This sensitivity is enough
to resolve the individual lensing images, as discussed later on in this
section.

An additional step is to correct the source’s frequency evolution Brady et al.
(1998). This is done by demodulating the signal via a time coordinate
redefinition $t\to t-t_{0}+\dot{f}_{0}/(2f_{0})(t-t_{0})^{2}+\cdots$.
Frequency evolution affects the CW’s phase on a timescale
$(\dot{f}_{0}/2)^{-1/2}=16.4\,{\rm d}\left(\frac{1{\rm pHz}/{\rm
s}}{\dot{f}_{0}}\right)^{1/2}$, making demodulation a necessary step of the
analysis. However, the interplay between lensing and $\dot{f}_{0}$ for rapidly
rotating NSs is typically negligible: it can produce additional modulation,
due to a phase difference between images on a timescale of $(\dot{f}_{0}\Delta
t_{\pm})^{-1}=317{\rm yr}\left(\frac{100{\rm s}}{\Delta
t_{\pm}}\right)\left(\frac{1{\rm pHz}/{\rm s}}{\dot{f}_{0}}\right)$. For ms
pulsars ($\dot{f}\lesssim 1$pHz/s) Woan et al. (2018), this is much longer
than $T_{\rm obs}$ and treating frequency evolution as common to all images is
an excellent approximation. Hence, demodulation does not affect the search of
lensed CW or their analysis.

Resolving individual images As anticipated in the section above, after
detection, accurate sky localization of the source can be achieved by
exploiting the full coherence of the observed signal to extract the source
direction $\hat{n}$ from the modulation produced by Earth’s rotation and
orbital motion. We will neglect the contribution from rotation, as it becomes
subdominant after a few days. To estimate the sky localisation accuracy, we
compute the mismatch between the two waveforms $(\delta h|\delta h)$. In our
case, $\delta h$ is the difference between two monochromatic waveforms
resampled with different source directions, $({\hat{n}},{\hat{n}}^{\prime})$,
so

$(\delta h|\delta h)\simeq
2\rho^{2}\left[1-\frac{1}{T}\int_{0}^{T}dt\cos{\left(2\pi
f_{0}\bm{x}(t)\cdot({\hat{n}}-{\hat{n}}^{\prime})\right)}\right].$ (21)

Imposing the condition $(\delta h|\delta h)=1$ returns the sky localization
variance Lindblom et al. (2008), $\sigma_{\theta}$, as a function of the sky
orientation of the vector ${\hat{n}}-{\hat{n}}^{\prime}$.

In the context of lensed sources, the accuracy of determining their sky
localization varies by the SNR of the observed images,
$\rho_{i}=\sqrt{\mu_{i}}\rho$. Fig. 7 shows the 1-$\sigma$ angular accuracy on
the image sky position that can be achieved for a CW signal lensed by SgrA∗
for our fiducial setup. The shape of the skymap is a squeezed cardioid that is
rescaled and flipped (because of spatial parity) between the two images. We
employed the module Astropy Astropy Collaboration et al. (2022) to model
Earth’s orbital motion relative to Sgr A∗ position in Eq.(21). Sgr A∗ is put
at a right ascension of 17h 45m 40.0409s and a declination of -29∘ 0’ 28.118”
Reid and Brunthaler (2004).

Figure 7:  Accuracy in resolving image sky positions. The figure displays the
1-$\sigma$ angular precision for determining the sky position of images
resulting from CWs sources lensed by Sgr A∗. The blue curve corresponds to the
brighter image, while the red curve to the fainter image. Both are rescaled by
the respective magnification factor. The resulting skymaps are shown in the
inset for illustrative purposes.

Lensing mimickers Situations that mimic strongly lensed sources can be ruled
out. A source in a binary may exhibit an amplitude modulation via precession
Apostolatos et al. (1994); Breton et al. (2008). However, this can be
distinguished from lensing through a periodic phase difference, which is
absent for lensed sources moving at constant velocity.

| spherical | disk
---|---|---

spatial | 0.059 | 771.6
frequency | 7.02
both | $4.1\cdot 10^{-9}$ | $5.4\cdot 10^{-5}$
Table 1: Expected number of overlapping signals that may mimic strong lensing
signatures. All these quantities scale with the number of detectable NSs as
$N_{0}/10^{8}$.

Finally, two unrelated sources may be close enough in frequency and sky
localization to appear as two images of a strongly lensed signal. However, the
chances of such a coincidence are negligible: to mimic strong lensing, their
angular separation must be $\lesssim 10^{\prime\prime}$ (cf. Fig. 4 in the
main manuscript) and their frequencies might differ by no more than $\Delta
f/f\lesssim 10^{-7}(\alpha_{ij}/10^{\prime\prime})(v/3\times 10^{-3})$, see
Eq. (3) in the main manuscript. Table 1 gives the expected number of
overlapping signals out of a total of $N_{0}=10^{8}$ detectable sources. For
the spatial overlap, we consider sources randomly distributed in the sky,
whose rate of overlap is $\sim N_{0}\delta\theta^{2}/4$, as well as sources
confined to the galactic disk $\sim N_{0}\theta/2\pi$ (both cases corresponds
to the limit of large and small $\Delta z$ in Eq. (32) below). Regarding the
frequency, we will consider sources distributed homogeneously in the range
$\log(f_{0})=[0.1,1]$ kHz, so the overlap rate is $\sim N_{0}\Delta
f/(f\log(10))$. While the chance of spatial or frequency overlap is sizeable,
the probability of both occurring simultaneously is negligible.

Additional information can be further used to constrain this possibility,
including the alignment of the images relative to Sgr A*, the signals’
relative amplitude, and the frequency evolution.

##  Appendix D Parameter inference

The signal observed at the detector, $d(t)$, is a superposition between
background noise and the true gravitational waves signal,

$d(t)=h(t;\bm{\theta})+n(t).$ (22)

The vector $\bm{\theta}$ includes the parameters that determine the propagated
waveform. For a stationary and Gaussian distributed noise, a realisation
$n_{0}$ has a probability

$p(n_{0})\propto\exp\left[-\frac{1}{2}\int
df\frac{|\tilde{n}_{o}(f)|^{2}}{S_{n}(f)}\right]\,,$ (23)

where the one-sided power spectral density of noise, $S_{n}(f)$, is the
variance associated with the noise distribution.

The probability of observing $d$ given $\bm{\theta}$, i.e. the _likelihood_ of
the parameters $\bm{\theta}$, is

$\log
p(d\mid\bm{\theta})\propto-\frac{1}{2}(d-h(\bm{\theta}),d-h(\bm{\theta}))\,,$
(24)

where the inner product for two signals $f$ and $g$, in the Fourier space, is
defined as

$(f,g)=4\,\text{Re}\left[\int_{0}^{\infty}df\,\frac{\tilde{h}^{\star}(f)\,\tilde{g}(f)}{S_{n}(f)}\right]\,.$
(25)

For quasi-monochromatic sources, the inner product in Eq. (25) can be
equivalently defined in the time-domain, following Ref. Takahashi and Seto
(2002); Seto (2002), as

$(h,g)=\frac{2}{S_{n}(f_{0})}\int_{0}^{T_{\mathrm{o}bs}}h(t)g(t)\,,$ (26)

with $T_{\mathrm{o}bs}$ the observational time. Hence, the SNR, which measures
the loudness of the signal, reads

$\rho\equiv(h|h)\simeq\sqrt{\frac{\mathcal{A}^{2}}{r_{s}^{2}}\frac{T_{\rm
obs}}{S_{n}(f_{0})}}\;.$ (27)

The posterior of $\bm{\theta}$ is obtained from the likelihood using Bayes’
theorem, i.e.,

$p(\bm{\theta}\mid
d)=\frac{p(d\mid\bm{\theta})p(\bm{\theta})}{p(d|\mathcal{M}_{i})}\;,$ (28)

where $p(\bm{\theta})$ is the parameter prior and $p(d|\mathcal{M}_{i})$,
namely the _evidence_ , is the marginal likelihood for a given model
$\mathcal{M}_{i}$:

$p(d|\mathcal{M}_{i})=\int d\bm{\theta}p(d|\bm{\theta})p(\bm{\theta})\,.$ (29)

In the Bayesian context, the marginalized posteriors are the probability
distribution functions of the parameters.

Given two models $\mathcal{M}_{i}$ and $\mathcal{M}_{j}$, the ratio of their
evidence called the Bayes factor, is used as an index to test different
interpretations of the same data:

$\mathcal{B}_{ij}=\frac{p(d|\mathcal{M}_{i})}{p(d|\mathcal{M}_{j})}\,.$ (30)

Following the empirical Jeffrey’s prescription Kass and Raftery (1995), in
Tab. 2 we report the range of value of $\mathcal{B}_{ij}$ and the
corresponding degree of evidence of the model $\mathcal{M}_{i}$ over
$\mathcal{M}_{j}$.

$\log\mathcal{B}_{10}$ |  $\mathcal{B}_{10}$ |  Evidence against $H_{0}$
---|---|---
2 to 6 |  3 to 20 |  Positive
6 to 10 |  20 to 150 |  Strong
$>$10 |  $>$150 |  Very Strong
Table 2: Evidence against the null hypothesis compared to ranges of Bayes
factor.

##  Appendix E Sampling methods

Combining the lensed waveform in Eq.(19) with Eq. (26) we get the analytical
likelihood function for CWs events lensed by SgrA∗. We use the dynamic nested
sampling package Dynesty Speagle (2020) to infer parameter posteriors and the
evidence of fitting models through numerical MC simulations.

Our analysis is restricted to the source and lens parameters that are
potentially correlated. In particular, the source is assumed to be purely
monochromatic and modelled by $({{\cal A}/r_{s},f_{0},\phi_{0}})$. For the MC
sampling, the lens is modelled by :

$u=\sqrt{\frac{\mu_{2}}{\mu_{1}}},\quad k=\Delta\phi_{21},\quad z=\frac{\Delta
f_{21}}{f}\,,$ (31)

where $\Delta\phi_{21}$ is the constant phase difference between the two
images and $\Delta f_{21}$ is the relative frequency shift induced by lens
motion. Initially, we use these model-independent parameters to avoid
complications that can arise from the periodicity of the likelihood with
respect to $M_{l}$ and $v$ and from the unboundedness of $y_{0}$. The
numerical samples are then converted to the point-lens parameters
$\\{M_{l},y_{0},v\\}$.

In Fig. 9, we present the complete set of 1-d and 2-d marginalized 2-$\sigma$
posteriors for the fiducial setup, already partially presented in Fig. 3 of
the main text. Figure 10 shows the posteriors of the corresponding original
sampling parameters $\\{u,k,z\\}$. For impact parameter $y\gtrsim 4$, the lack
of resolution of the small amplitude modulation, i.e. the $r$ contour being
compatible with $r=0$, prevents the reconstruction of the lensing parameters.

##  Appendix F Lensing probability

We will now estimate the probability of strong lensing and compare it with
previous results. We will follow the source distribution proposed in Ref. Reed
et al. (2021) (also considered in Ref. Basak et al. (2023)),

$\frac{dP_{s}(r,z)}{dV}=\frac{1}{2\pi}\frac{1}{\sigma_{r}^{2}}e^{-\frac{r^{2}}{2\sigma_{r}^{2}}}\frac{1}{2\Delta
z}e^{-\frac{|z|}{\Delta z}}\,.$ (32)

Here $dV=drrd\phi dz$ is the volume element, $\sigma_{r}=5$kpc is the radial
scatter of sources and $\Delta z$ is the scatter perpendicular to the galactic
plane. We will assume that the Earth and any detectable lensed source deviate
negligibly from the galactic plane, relative to $\Delta z\sim 0.1-1$kpc, and
set $z\approx 0$.

The fraction of strongly lensed sources is

$f_{l}\equiv\frac{\bar{N}_{l}}{N_{0}}\approx\int_{0}^{r_{\rm
hor}}dr_{ls}r_{ls}(r_{l}+r_{ls})\left(\theta_{E}y_{\rm
max}\right)^{2}\frac{dP_{s}}{dV}\,,$ (33)

were $\bar{N}_{l}$ is the average number of lensed sources and $N_{0}$ is the
number of observable sources Basak et al. (2023). The integral is performed up
to a detection horizon $r_{\rm hor}$ (cf. Eq. 20) and $y_{\rm max}$ is the
highest value of the impact parameter for which very strong evidence for
lensing can be established, see Fig. 2 in the main manuscript. For
$r_{ls}/r_{l}\geq 0.5$, $y_{\rm max}$ obeys a linear relation, which we have
extrapolated to the whole domain to compute the integral above. This
underestimates $y_{\rm max}$ close to the lens (stars vs line in Fig. 8) and
is thus a conservative assumption.

Figure 8: Expected fraction of NSs lensed by Sgr A∗ using the strong-evidence
cutoff (thick solid red) and the Einstein radius cutoff (thick solid grey) as
a function of $r_{\rm hor}/r_{l}-1$ for $\Delta z=0.1$kpc. Thin lines
represent other quantities appearing in Eq. (33), plotted as a function of
$r_{ls}/r_{l}$: Einstein angle (dashed), source density (dash-dotted), values
of $y_{\rm max}$ interpolated (dashed) from the results of the MC sampling
(stars).

Figure 8 shows the fraction of lensed NSs with strong evidence, $f_{l}[y_{\rm
max}]$, as a function of $r_{\rm hor}$, where we also show the fraction of
sources within the Einstein radius, ${f_{l}[y=1]}$, obtained from setting
$y_{\rm max}\to 1$ in Eq. (33). The Einstein radius, source distribution and
strong-evidence parameter entering the calculation are also shown. For $r_{\rm
hor}\gtrsim 1.5\,r_{l}$, the expected number of observed sources is

${\bar{N}_{L}}=3.36\left(\frac{0.1{\rm kpc}}{\Delta
z}\right)\left(\frac{N_{0}}{10^{8}}\right)\,.$ (34)

This is a factor $\sim 9.04$ larger than the estimate obtained setting $y_{\rm
max}=1$ in Eq. (33), as in the analysis of Ref. Basak et al. (2023). The
number of detected sources is given by Poisson statistics, with

$P_{l}(k)=\frac{\bar{N}_{l}^{k}}{k!}e^{-\bar{N}_{l}}\,,$ (35)

where $k$ is the number of detections. Hence, the probability of detecting at
least a single lensed event is $P_{l}(k\geq 1)=1-e^{-N_{l}}$. The analysis in
Ref. Basak et al. (2023) uses $N_{0}=10^{9}$ as a fiducial number of
observable sources. In that case, and for a narrow spread $\Delta z=0.1$ kpc,
our study suggests that an average of 34 strong-evidence lensed events would
be observed, with $\sim 2\cdot 10^{-15}$ chance of observing none.

## References

  * Bartelmann (2010) M. Bartelmann, Class. Quant. Grav. 27, 233001 (2010), eprint 1010.3829.
  * Oguri (2019) M. Oguri, Rept. Prog. Phys. 82, 126901 (2019), eprint 1907.06830.
  * Liao et al. (2022) K. Liao, M. Biesiada, and Z.-H. Zhu, Chin. Phys. Lett. 39, 119801 (2022), eprint 2207.13489.
  * Takahashi and Nakamura (2003) R. Takahashi and T. Nakamura, Astrophys. J. 595, 1039 (2003), eprint astro-ph/0305055.
  * Dai et al. (2018) L. Dai, S.-S. Li, B. Zackay, S. Mao, and Y. Lu, Phys. Rev. D 98, 104029 (2018), eprint 1810.00003.
  * Çalışkan et al. (2022) M. Çalışkan, L. Ji, R. Cotesta, E. Berti, M. Kamionkowski, and S. Marsat (2022), eprint 2206.02803.
  * Tambalo et al. (2022a) G. Tambalo, M. Zumalacárregui, L. Dai, and M. H.-Y. Cheung (2022a), eprint 2210.05658.
  * Tambalo et al. (2022b) G. Tambalo, M. Zumalacárregui, L. Dai, and M. H.-Y. Cheung (2022b), eprint 2212.11960.
  * Dai and Venumadhav (2017) L. Dai and T. Venumadhav (2017), eprint 1702.04724.
  * Ezquiaga et al. (2021) J. M. Ezquiaga, D. E. Holz, W. Hu, M. Lagos, and R. M. Wald, Phys. Rev. D 103, 064047 (2021), eprint 2008.12814.
  * Vijaykumar et al. (2022) A. Vijaykumar, A. K. Mehta, and A. Ganguly (2022), eprint 2202.06334.
  * Hannuksela et al. (2019) O. A. Hannuksela, K. Haris, K. K. Y. Ng, S. Kumar, A. K. Mehta, D. Keitel, T. G. F. Li, and P. Ajith, Astrophys. J. Lett. 874, L2 (2019), eprint 1901.02674.
  * Dai et al. (2020) L. Dai, B. Zackay, T. Venumadhav, J. Roulet, and M. Zaldarriaga (2020), eprint 2007.12709.
  * Abbott et al. (2021) R. Abbott et al. (LIGO Scientific, VIRGO), Astrophys. J. 923, 14 (2021), eprint 2105.06384.
  * Basak et al. (2022) S. Basak, A. Ganguly, K. Haris, S. Kapadia, A. K. Mehta, and P. Ajith, Astrophys. J. Lett. 926, L28 (2022), eprint 2109.06456.
  * Ng et al. (2018) K. K. Y. Ng, K. W. K. Wong, T. Broadhurst, and T. G. F. Li, Phys. Rev. D 97, 023012 (2018), eprint 1703.06319.
  * Xu et al. (2022) F. Xu, J. M. Ezquiaga, and D. E. Holz, Astrophys. J. 929, 9 (2022), eprint 2105.14390.
  * Amaro-Seoane et al. (2017) P. Amaro-Seoane et al. (LISA) (2017), eprint 1702.00786.
  * Sesana (2016) A. Sesana, Phys. Rev. Lett. 116, 231102 (2016), eprint 1602.06951.
  * Wagg et al. (2022) T. Wagg, F. S. Broekgaarden, S. E. de Mink, L. A. C. van Son, N. Frankel, and S. Justham, Astrophys. J. 937, 118 (2022), eprint 2111.13704.
  * Kalogera et al. (2021) V. Kalogera et al. (2021), eprint 2111.06990.
  * Miller et al. (2021) A. L. Miller, S. Clesse, F. De Lillo, G. Bruno, A. Depasse, and A. Tanasijczuk, Phys. Dark Univ. 32, 100836 (2021), eprint 2012.12983.
  * Abbott et al. (2022a) R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA) (2022a), eprint 2201.00697.
  * Steltner et al. (2023) B. Steltner, M. A. Papa, H. B. Eggenstein, R. Prix, M. Bensch, B. Allen, and B. Machenschalk, Astrophys. J. 952, 55 (2023), eprint 2303.04109.
  * Dergachev and Papa (2022) V. Dergachev and M. A. Papa (2022), eprint 2202.10598.
  * Riles (2017) K. Riles, Mod. Phys. Lett. A 32, 1730035 (2017), eprint 1712.05897.
  * Sieniawska and Bejger (2019) M. Sieniawska and M. Bejger, Universe 5, 217 (2019), URL https://doi.org/10.3390%2Funiverse5110217.
  * Piccinni (2022) O. J. Piccinni, Galaxies 10, 72 (2022), eprint 2202.01088.
  * Riles (2022) K. Riles (2022), eprint 2206.06447.
  * Liao et al. (2019) K. Liao, M. Biesiada, and X.-L. Fan, Astrophys. J. 875, 139 (2019), eprint 1903.06612.
  * Suvorov (2022) A. G. Suvorov, Astrophys. J. 930, 13 (2022), eprint 2112.01670.
  * Moylan et al. (2007) A. J. Moylan, D. E. McClelland, S. M. Scott, A. C. Searle, and G. V. Bicknell, in _11th Marcel Grossmann Meeting on General Relativity_ (2007), pp. 807–823, eprint 0710.3140.
  * Jung and Kim (2023) S. Jung and S. Kim, JCAP 07, 042 (2023), eprint 2210.02649.
  * Takahashi et al. (2023) R. Takahashi, S. Morisaki, and T. Suyama (2023), eprint 2304.08220.
  * Narayan (1987) R. Narayan, Astrophys. J. 319, 162 (1987).
  * Ghez et al. (2000) A. Ghez, M. Morris, E. E. Becklin, T. Kremenek, and A. Tanner, Nature 407, 349 (2000), eprint astro-ph/0009339.
  * Genzel et al. (2010) R. Genzel, F. Eisenhauer, and S. Gillessen, Rev. Mod. Phys. 82, 3121 (2010), eprint 1006.0064.
  * Akiyama et al. (2022) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 930, L12 (2022).
  * Basak et al. (2023) S. Basak, A. K. Sharma, S. J. Kapadia, and P. Ajith, Astrophys. J. Lett. 942, L31 (2023), eprint 2205.00022.
  * Reed et al. (2021) B. T. Reed, A. Deibel, and C. J. Horowitz, Astrophys. J. 921, 89 (2021), eprint 2104.00771.
  * Harding (2013) A. K. Harding, Front. Phys. (Beijing) 8, 679 (2013), eprint 1302.0869.
  * Kayser et al. (1986) R. Kayser, S. Refsdal, and R. Stabell, Symposium - International Astronomical Union 166, 36 (1986).
  * and R. Abuter et al. (2019) and R. Abuter, A. Amorim, M. Bauböck, J. P. Berger, H. Bonnet, W. Brandner, Y. Clé net, V. C. du Foresto, P. T. de Zeeuw, J. Dexter, et al., Astronomy & Astrophysics 625, L10 (2019), URL https://doi.org/10.1051%2F0004-6361%2F201935656.
  * Brady et al. (1998) P. R. Brady, T. Creighton, C. Cutler, and B. F. Schutz, Phys. Rev. D 57, 2101 (1998), eprint gr-qc/9702050.
  * Maggiore (2007) M. Maggiore, _Gravitational Waves. Vol. 1: Theory and Experiments_ , Oxford Master Series in Physics (Oxford University Press, 2007), ISBN 978-0-19-857074-5, 978-0-19-852074-0.
  * Abuter et al. (2020) R. Abuter et al. (GRAVITY), Astron. Astrophys. 636, L5 (2020), eprint 2004.07187.
  * Naoz et al. (2020) S. Naoz, C. M. Will, E. Ramirez-Ruiz, A. Hees, A. M. Ghez, and T. Do, Astrophys. J. Lett. 888, L8 (2020), eprint 1912.04910.
  * Abuter et al. (2022) R. Abuter et al. (GRAVITY), Astron. Astrophys. 657, L12 (2022), eprint 2112.07478.
  * Will et al. (2023) C. M. Will, S. Naoz, A. Hees, A. Tucker, E. Zhang, T. Do, and A. Ghez (2023), eprint 2307.16646.
  * Bozza and Mancini (2012) V. Bozza and L. Mancini, Astrophys. J. 753, 56 (2012), eprint 1204.2103.
  * Gondolo and Silk (1999) P. Gondolo and J. Silk, Phys. Rev. Lett. 83, 1719 (1999), eprint astro-ph/9906391.
  * Hui et al. (2017) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D 95, 043541 (2017), eprint 1610.08297.
  * Ezquiaga and Zumalacárregui (2020) J. M. Ezquiaga and M. Zumalacárregui, Phys. Rev. D 102, 124048 (2020), eprint 2009.12187.
  * Goyal et al. (2023) S. Goyal, A. Vijaykumar, J. M. Ezquiaga, and M. Zumalacarregui (2023), eprint 2301.04826.
  * Gong et al. (2021) Y. Gong, J. Luo, and B. Wang, Nature Astronomy 5, 881 (2021).
  * Sedda et al. (2020) M. A. Sedda et al., Classical and Quantum Gravity 37, 215011 (2020).
  * Baibhav et al. (2021) V. Baibhav et al., Experimental Astronomy 51, 1385 (2021).
  * Nakamura and Deguchi (1999) T. T. Nakamura and S. Deguchi, Progress of Theoretical Physics Supplement 133, 137 (1999), ISSN 0375-9687, eprint https://academic.oup.com/ptps/article-pdf/doi/10.1143/PTPS.133.137/5283012/133-137.pdf, URL https://doi.org/10.1143/PTPS.133.137.
  * Schneider et al. (1992) P. Schneider, J. Ehlers, and E. E. Falco, _Gravitational Lenses_ (1992).
  * An and Evans (2006) J. H. An and N. W. Evans, Mon. Not. Roy. Astron. Soc. 369, 317 (2006), eprint astro-ph/0601457.
  * Glampedakis and Gualtieri (2018) K. Glampedakis and L. Gualtieri, in _The Physics and Astrophysics of Neutron Stars_ (Springer International Publishing, 2018), pp. 673–736, URL https://doi.org/10.1007%2F978-3-319-97616-7_12.
  * Pagliaro et al. (2023) G. Pagliaro, M. A. Papa, J. Ming, J. Lian, D. Tsuna, C. Maraston, and D. Thomas, Astrophys. J. 952, 123 (2023), eprint 2303.04714.
  * Abbott et al. (2022b) R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA) (2022b), eprint 2204.04523.
  * Suvorov et al. (2016) A. G. Suvorov, A. Mastrano, and U. Geppert, Monthly Notices of the Royal Astronomical Society 459, 3407 (2016), ISSN 0035-8711, eprint https://academic.oup.com/mnras/article-pdf/459/3/3407/8110424/stw909.pdf, URL https://doi.org/10.1093/mnras/stw909.
  * Gittins et al. (2020) F. Gittins, N. Andersson, and D. I. Jones, Monthly Notices of the Royal Astronomical Society 500, 5570 (2020), URL https://doi.org/10.1093%2Fmnras%2Fstaa3635.
  * Woan et al. (2018) G. Woan, M. D. Pitkin, B. Haskell, D. I. Jones, and P. D. Lasky, Astrophys. J. Lett. 863, L40 (2018), eprint 1806.02822.
  * Sartore et al. (2010) N. Sartore, E. Ripamonti, A. Treves, and R. Turolla, Astron. Astrophys. 510, A23 (2010), eprint 0908.3182.
  * Hobbs et al. (2005) G. Hobbs, D. R. Lorimer, A. G. Lyne, and M. Kramer, Mon. Not. Roy. Astron. Soc. 360, 974 (2005), eprint astro-ph/0504584.
  * Sofue et al. (2009) Y. Sofue, M. Honma, and T. Omodaka, Publ. Astron. Soc. Jap. 61, 227 (2009), eprint 0811.0859.
  * Mró z et al. (2019) P. Mró z, A. Udalski, D. M. Skowron, J. Skowron, I. Soszyński, P. Pietrukowicz, M. K. Szymański, R. Poleski, S. Kozłowski, and K. Ulaczyk, The Astrophysical Journal 870, L10 (2019), URL https://doi.org/10.3847%2F2041-8213%2Faaf73f.
  * Eilers et al. (2019) A.-C. Eilers, D. W. Hogg, H.-W. Rix, and M. K. Ness, The Astrophysical Journal 871, 120 (2019), URL https://doi.org/10.3847%2F1538-4357%2Faaf648.
  * Piccinni et al. (2020) O. J. Piccinni, P. Astone, S. D’Antonio, S. Frasca, G. Intini, I. La Rosa, P. Leaci, S. Mastrogiovanni, A. Miller, and C. Palomba, Phys. Rev. D 101, 082004 (2020), eprint 1910.05097.
  * Lindblom et al. (2008) L. Lindblom, B. J. Owen, and D. A. Brown, Phys. Rev. D 78, 124020 (2008), eprint 0809.3844.
  * Astropy Collaboration et al. (2022) Astropy Collaboration, A. M. Price-Whelan, P. L. Lim, N. Earl, N. Starkman, L. Bradley, D. L. Shupe, A. A. Patil, L. Corrales, C. E. Brasseur, et al., Astrophys. J. 935, 167 (2022), eprint 2206.14220.
  * Reid and Brunthaler (2004) M. J. Reid and A. Brunthaler, Astrophys. J. 616, 872 (2004), eprint astro-ph/0408107.
  * Apostolatos et al. (1994) T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Phys. Rev. D 49, 6274 (1994).
  * Breton et al. (2008) R. P. Breton, V. M. Kaspi, M. Kramer, M. A. McLaughlin, M. Lyutikov, S. M. Ransom, I. H. Stairs, R. D. Ferdman, F. Camilo, and A. Possenti, Science 321, 104 (2008), eprint 0807.2644.
  * Takahashi and Seto (2002) R. Takahashi and N. Seto, Astrophys. J. 575, 1030 (2002), arXiv:astro-ph/0204487.
  * Seto (2002) N. Seto, Mon. Not. Roy. Astron. Soc. 333, 469 (2002), arXiv:astro-ph/0202364.
  * Kass and Raftery (1995) R. E. Kass and A. E. Raftery, J. Am. Statist. Assoc. 90, 773 (1995).
  * Speagle (2020) J. S. Speagle, Monthly Notices of the Royal Astronomical Society 493, 3132 (2020), URL https://doi.org/10.1093%2Fmnras%2Fstaa278.

Figure 9: 1-d and 2-d marginalized 2-$\sigma$ posteriors of the full set of
source and lens parameters, for the fiducial setup discussed in the text.
Here, $A\equiv\mathcal{A}/r_{s}$. Figure 10: 1-d and 2-d marginalized
2-$\sigma$ posteriors of the full set of source and lens parameters, with the
letter expressed in terms of the original sampling parameters, for the
fiducial setup discussed in the text. Here, $A\equiv\mathcal{A}/r_{s}$.