# Exploring the hidden Universe: A novel phenomenological approach for
recovering arbitrary gravitational-wave millilensing configurations

Anna Liu,1 Isaac C. F. Wong,1 Samson H. W. Leong, 1 Anupreeta More, 2,3 Otto
A. Hannuksela,1 and Tjonnie G. F. Li1,4,5
1 Department of Physics, The Chinese University of Hong Kong, Shatin, NT, Hong
Kong
2 The Inter-University Centre for Astronomy and Astrophysics (IUCAA), Post Bag
4, Ganeshkhind, Pune 411007, India
3 Kavli Institute for the Physics and Mathematics of the Universe (IPMU),
5-1-5 Kashiwanoha, Kashiwa-shi, Chiba 277-8583, Japan
4 Institute for Theoretical Physics, KU Leuven, Celestijnenlaan 200D, B-3001
Leuven, Belgium
5 Department of Electrical Engineering (ESAT), KU Leuven, Kasteelpark Arenberg
10, B-3001 Leuven, Belgium E-mail<EMAIL_ADDRESS>(AL)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Since the first detection of gravitational waves in 2015, gravitational-wave
astronomy has emerged as a rapidly advancing field that holds great potential
for studying the cosmos, from probing the properties of black holes to testing
the limits of our current understanding of gravity. One important aspect of
gravitational-wave astronomy is the phenomenon of gravitational lensing, where
massive intervening objects can bend and magnify gravitational waves,
providing a unique way to probe the distribution of matter in the universe, as
well as finding applications to fundamental physics, astrophysics, and
cosmology. However, current models for gravitational-wave millilensing - a
specific form of lensing where small-scale astrophysical objects can split a
gravitational wave signal into multiple copies - are often limited to simple
isolated lenses, which is not realistic for complex lensing scenarios. In this
paper, we present a novel phenomenological approach to incorporate
millilensing in data analysis in a model-independent fashion. Our approach
enables the recovery of arbitrary lens configurations without the need for
extensive computational lens modeling, making it a more accurate and
computationally efficient tool for studying the distribution of matter in the
universe using gravitational-wave signals. When gravitational-wave lensing
observations become possible, our method could provide a powerful tool for
studying complex lens configurations in the future.

###### keywords:

gravitational waves – gravitational lensing: micro – gravitational lensing:
strong

††pubyear: 2023††pagerange: Exploring the hidden Universe: A novel
phenomenological approach for recovering arbitrary gravitational-wave
millilensing configurations–Two-signal example Figure 1: An example
configuration of millilenses within a galaxy producing multiple millisignals
and influencing each other gravitationally (figure not to scale); (a) as a
gravitational wave passes near a massive object such as a galaxy or a galaxy
cluster, its path gets bent which can result in multiple strongly-lensed GW
signals (represented as curved black lines) reaching the observer; (b)
additionally, if the signals encounter smaller compact objects, e.g., stars,
dark matter subhalos, or massive compact halo objects, acting as millilenses,
further splitting of the signals occurs, which results in multiple lensed BBH
images arising in the lens plane. The circular curves around each lens are the
(tangential) critical curves. Throughout this work, we use the thin lens
approximation in which the source and lens objects are confined to two-
dimensional planes.

##  1 Introduction

Gravitational waves (GWs) were predicted by Albert Einstein’s theory of
general relativity in 1916. However, it took nearly a century before their
direct detection by the LIGO and Virgo Collaborations in 2015, opening a new
era of gravitational-wave astronomy (Abbott et al., 2016a, b). Since then,
numerous GW detections have been made, including binary black hole (BBH)
mergers, binary neutron star (BNS) mergers, and the merger of a BBH with a
neutron star (BH-NS) (Abbott et al., 2019, 2021d, 2021a) with the Advanced
LIGO (Aasi et al. (2015)), Advanced Virgo (Acernese et al. (2014)), and KAGRA
detectors (Akutsu et al. (2019)).

Gravitational lensing, the bending of light or GWs by a massive object, was
first observed in 1919 by Sir Arthur Eddington during a total solar eclipse
(Dyson et al., 1920), which confirmed Einstein’s theory of general relativity.
In recent years, the search for GW lensing signatures has become a fast-
developing field. Multiple lensed signals can arrive at the observer,
described as time-separated de/magnified GW signals with a phase shift
relative to the unlensed GW signal (Takahashi & Nakamura, 2003; Cao et al.,
2014; Dai & Venumadhav, 2017). The current gravitational-wave data has been
used to search for GW lensing, but so far, there has not yet been any widely
accepted detection.

The challenges in making GW lensing detections are significant, including the
need for new tools to detect strong lensing (Haris et al., 2018; Dai et al.,
2020; Liu et al., 2020; Lo & Hernandez, 2021; Janquart et al., 2021a, 2022),
micro- and millilensing (Lai et al., 2018; Dai et al., 2018; Liao et al.,
2018; Christian et al., 2018; Pagano et al., 2020; Kim et al., 2020; Seo et
al., 2021; Wright & Hendry, 2021), and wave optics lensing (Lai et al., 2018;
Christian et al., 2018; Bulashenko & Ubach, 2021; Oguri & Takahashi, 2022;
Basak et al., 2022b; Tambalo et al., 2022a), beyond the need for statistical
modelling of both gravitational lenses and binary black holes (Dai et al.,
2016; Ng et al., 2017; Smith et al., 2017; Smith et al., 2018a; Oguri, 2018;
Wierda et al., 2021a; Xu et al., 2021; More & More, 2021; Smith et al., 2022).
Furthermore, the low expected rates of lensed GW events make detection
difficult (Çalışkan et al., 2022). Nevertheless, recent studies have shown
that the detection of lensed GW events is possible with current and future GW
detectors, assuming the current best models for the strong lensing galaxy
population and presuming that binary black holes trace the star-formation rate
density (Ng et al., 2017; Li et al., 2018; Oguri, 2018; Xu et al., 2021;
Wierda et al., 2021a; Wempe et al., 2022; Smith et al., 2022; Ma et al.,
2022). Indeed, the era of gravitational wave lensing is imminent, and the
detection of lensed GW events would provide valuable information about the
properties of the lensing objects, test the theory of general relativity in a
new regime, probe high-redshift cosmology, study dark matter, search for
primordial black holes and MACHOs, precisely localize merging black holes, and
more (Takahashi, 2005; Itoh et al., 2009; Baker & Trodden, 2016; Collett &
Bacon, 2017; Liao et al., 2017; Fan et al., 2017; Lai et al., 2018; Dai et
al., 2018; Mukherjee et al., 2019, 2020; Diego, 2019; Oguri & Takahashi, 2020;
Goyal et al., 2020; Hannuksela et al., 2020; Finke et al., 2021; Iacovelli et
al., 2022; Chung & Li, 2021; Sereno et al., 2010; Bolejko et al., 2012; Magaña
Hernandez, 2022; Tambalo et al., 2022b; Basak et al., 2022b; Basak et al.,
2022a) and the searches for gravitational-wave lensing have started recently
Hannuksela et al. (2019); Abbott et al. (2021b); Li et al. (2019); Dai et al.
(2020); Liu et al. (2020); Pang et al. (2020); The LIGO Scientific
Collaboration et al. (2021); Kim et al. (2022).

Depending on the lens mass, gravitational lenses can cause different types of
lensing. Massive objects such as galaxies and galaxy clusters can produce
strong gravitational lensing, resulting in multiple signals with varying time
separations and effects on the GW signal (Takahashi & Nakamura, 2003; Dai &
Venumadhav, 2017; Smith et al., 2017; Smith et al., 2018b; Haris et al., 2018;
Liu et al., 2020; Robertson et al., 2020; Ryczanowski et al., 2020; Dai et
al., 2020; Wang et al., 2021; Ezquiaga et al., 2021; Lo & Hernandez, 2021;
Janquart et al., 2021a, b; Vijaykumar et al., 2022; Çalışkan et al., 2022; Cao
et al., 2022). On the other hand, smaller-mass lenses such as stars or
stellar-mass compact objects can cause microlensing, resulting in potentially
observable beating patterns in the frequency evolution of the GW waveform
(Deguchi & Watson, 1986; Nakamura, 1998; Takahashi & Nakamura, 2003; Christian
et al., 2018; Jung & Shin, 2019; Diego et al., 2019; Mishra et al., 2021;
Meena & Bagla, 2020; Cheung et al., 2020; Bulashenko & Ubach, 2021; Cremonese
et al., 2021; Seo et al., 2021; Yeung et al., 2021; Qiu et al., 2022; Wright &
Hendry, 2021; Kim et al., 2022). In this paper, microlensing refers to lens
masses of tens to hundreds of solar masses, while millilensing refers to lens
masses in the range of $10^{2}-10^{6}M_{\odot}$, covering dark matter
subhalos, primordial black holes, massive compact halo objects (MACHOs), and
lenses with an Einstein radius approximately the size of a milliarcsecond. The
distinction between microlensing and millilensing is based on the lens mass
range and the corresponding time delays between individual lensed GW signals.
The analysis of microlensing requires the wave optics approximation in the
low-mass regime, where the GW wavelength matches the Schwarzschild radius of
the lens, while we assume millilensing follows the geometrical optics
approximation in this work, which holds for the corresponding mass range
(Takahashi & Nakamura, 2003).

To date, several lens models have been suggested for GW micro- and
millilensing studies. These models typically assume a particular lens mass
distribution. Examples include the point-mass lens (PML) model and singular
isothermal sphere (SIS) model (Nakamura (1998); Takahashi & Nakamura (2003);
Congdon & Keeton (2018)), where the lenses are assumed to be isolated from any
astrophysical objects. However, past electromagnetic lensing observations have
shown that astrophysical objects, typically galaxies, can have a significant
impact on the lens morphology, making it difficult to treat the lens as an
isolated object (e.g. Diego et al., 2019; Seo et al., 2021; Oguri & Takahashi,
2022). For example, when a gravitational wave is strongly lensed by a galaxy
and micro- or millilensed by the small-scale structure within the galaxy, the
millilens would experience additional effects, including gravitational shear
due to the galaxy’s potential (see Fig. 1) (Mishra et al., 2021; Seo et al.,
2021; Oguri & Takahashi, 2022). This scenario could lead to a non-symmetric
distribution of the millilens mass, as shown in Fig. 2. The widely used
isolated spherically symmetric lens assumption is then physically unrealistic,
at least in scenarios with significant shearing effects. In such a scenario,
the current parameter estimation tools could obtain a significantly biased
result or even miss the signal altogether (Yeung et al., 2021).

One potential way to improve our ability to recover these more complex lenses
is to implement a collection of more intricate lens models in our parameter
estimation framework, including extensions that account for galaxy shearing
effects and fields of lenses, as well as several millilensing mass profiles.
However, this approach would require a significant amount of effort to extend
the parameter estimation framework to every plausible lensing scenario.
Additionally, to implement this approach, we would need to perform parameter
estimation on each signal for each lens model, making the process
computationally expensive. Instead of relying on complex and computationally
intensive lens models to account for physically realistic effects of
astrophysical objects on millilensed gravitational-wave signals, we propose a
different approach in this work.

In particular, we suggest using a general, model-independent description of
millilensed signals that parameterizes "image parameters" instead of the
lensing system itself. Every millilensing system produces an integer number of
gravitational-wave milli-images, each with their own set of magnifications,
time delays, and arrival times, regardless of the lens configuration
(Takahashi & Nakamura, 2003; Dai & Venumadhav, 2017). By parameterizing the
properties of these milli-images, we can target any millilensing
configuration, covering signals from different types of lenses and
configurations without imposing limits on the system setup or the number of
lensed GW signals formed. This phenomenological approach also allows us to
recycle the results from the analysis of multiple data sets, without
performing separate analyses with different lens models independently.
Furthermore, we demonstrate that the results obtained from the
phenomenological analysis parameterizing the lensed GW signal can be mapped to
specific lens models to select the most favorable lens mass profile.

Figure 2: Schematic representation of point-mass lens model and
phenomenological approach in the lens plane: (a) in the point-mass lens model,
the lens is represented as a point-like mass (purple) producing two GW
millisignals from the background BBH source. This model is often used in
microlensing analysis, where the size of the lens is small compared to the
relative size of the system and the critical curve of the lens is circular;
(b) in the phenomenological approach presented here, we assume multiple
millilenses can influence each other gravitationally, producing an overall
gravitational shearing effect, hence breaking the spherical symmetry of the
millilenses. Additionally, unlike the point-mass lens model, we assume each
millilens can produce an arbitrary integer number of millisignals.

##  2 Methods

###  2.1 Lens models

In order to demonstrate the effectiveness of our phenomenological model for
searching generic millilensing configurations, we need to test it by injecting
lensed GW signals into LVK detector data. To simulate these signals, we must
choose a specific lens model, and for this exercise, we will focus on well-
understood models that can be easily tested.

Therefore, we consider the point mass lens (PML) and singular isothermal
sphere (SIS) models as our millilenses (Nakamura, 1998; Takahashi & Nakamura,
2003; Congdon & Keeton, 2018). We also incorporate an external shearing effect
when these models are embedded into a galaxy. This effect stretches the
caustic curves and can lead to additional milli-images (Diego et al., 2019;
Diego, 2019). While the PML model with a shearing effect is commonly used to
approximate compact objects acting as gravitational lenses, the SIS model
assumes a spherically symmetric mass distribution and is often used as a toy
model to represent galaxies, star clusters, or dark matter subhalos. Note that
more complex models, such as tidally stripped Navarro Fenk White lenses, may
be required to accurately model dark matter subhalos, but our primary goal is
to demonstrate a phenomenological search that can be applied to any lens
model. We use the Lenstronomy package, a multi-purpose Python package for
gravitational lensing, to perform all lens modeling (Birrer & Amara, 2018).

In addition to choosing a lens model, simulating the lensed GW waveform is
also necessary. A lensed GW can be mathematically expressed as
$\tilde{h}_{\textrm{L}}(f)\\!=\\!F(f)\tilde{h}(f)$, where $F(f)$ is the
amplification factor, which is a function of frequency $f$ and specific lens
parameters, and $\tilde{h}(f)$ represents the unlensed GW. In the case of the
PML and SIS models, the amplification factor is dependent on two lensed
parameters: the redshifted lens mass $M{Lz}$ and the relative position of the
source with respect to the lens, denoted by $y$. The addition of an external
shearing effect, as described in the previous paragraph, would introduce three
additional parameters: two components of the shearing effect, denoted as
$\gamma_{1,2}$, and a convergence value denoted as $\kappa$.

In the case of an isolated millilens, using the thin lens approximation and
solving the lens equation for the PML model gives rise to two solutions of
lensed GW signals (see the left panel in Fig.2). For relatively short time
delays, the millilensed signals will overlap, leading to a single lensed GW
signal arriving at the detector (see Fig.3). However, assuming the lens is
embedded within a galaxy and including the internal shear, the number of the
component millilensed GW signals can be different (Diego et al., 2019; Diego,
2019). Indeed, a generic lens potential may even include a field of small-
scale lenses.

###  2.2 Model-independent inference of millilensed GW signals

When considering millilensing, there are a multitude of potential lensing
configurations, each with varying levels of complexity. This can include
different types of millilenses, such as PMLs, SISs, or more complex lens
structures, as well as a varying number of such lenses. Additionally, the
effect of the strong lensing galaxy can further complicate the system. Due to
this complexity, it is impractical to implement every gravitational lens model
directly in parameter estimation. However, regardless of the specific
millilens model, the observable signal will always consist of an integer
number of GW millisignals, or "milli-images," each with their own properties
such as time delay, magnification, and Morse phase. These signals overlap at
the GW detector. In this work, we present a phenomenological search approach
that directly targets the milli-images, allowing for an arbitrary
configuration of the lensing system.

Figure 3: Illustration of the millilensing effect on the GW waveform. A GW
signal passing near a millilens is split into two millisignals that overlap,
resulting in a single GW signal at the detector with frequency-dependent
beating patterns. One can discern the millilensing effect by locating the
beating patterns in GWs. The thin lens approximation used confines the source,
the lens and the observer to a 2-dimensional plane each. The deflection of the
GW occurs instantaneously at the point where the GW crosses the lens plane.

When considering a millilens embedded in a macrosystem such as a galaxy, it is
important to take into account the effects of the macrosystem. In contrast to
previous gravitational wave analyses that have focused primarily on isolated
point mass lenses and singular isothermal sphere profiles described by
physical parameters of the lensing system (such as lens mass $M_{Lz}$ and
source position $y$) (e.g. The LIGO Scientific Collaboration et al., 2021), we
parameterize each individual millilensed GW image using a set of millilensing
parameters: relative magnification $\mu$, time delay $t$, and Morse factor
$n$. This parameterization can account for non-symmetric effects due to the
macrosystem without making prior assumptions about the system configuration
and mass distribution. Each individual millilensed GW signal is treated
independently of the other millisignals.

Before the individual millisignals reach the detector, they overlap due to
relatively short time separations (order of milliseconds). The resultant GW
signal can therefore be described as a sum of the millisignals with the
waveform shape characterised by frequency-dependent beating patterns (see Fig.
3).

To allow for the most general case where any integer number of millisignals
may be formed due to lensing, we introduce an additional lensing parameter,
$K$, which corresponds to the number of millisignals. Unlike traditional
models that assume a fixed number of millisignals based on a specific lensing
scenario, our phenomenological approach does not limit the number of signals.
By relaxing this assumption, we provide a more flexible framework for
characterizing the complex physical configuration of the lensing system. The
mathematical formulation of our approach is described below.

###  2.3 A phenomenological formulation

A gravitationally lensed waveform can be separated into an unlensed part
(unlensed complex strain amplitude) multiplied by an amplification factor
$F(f,\;\boldsymbol{\theta_{L}})$ which contains lensing information described
by lensing parameters $\boldsymbol{\theta_{L}}$. A millilensed waveform
$\tilde{h}_{L}(f;\boldsymbol{\theta},\;\boldsymbol{\theta_{L}})$ can thus be
expressed as

$\displaystyle\begin{aligned}
\tilde{h}_{L}(f;\boldsymbol{\theta},\;\boldsymbol{\theta_{L}})={F(f,\boldsymbol{\theta_{L}})}\cdot\tilde{h}_{U}(f;\boldsymbol{\theta}),\end{aligned}$
(1)

where $\tilde{h}_{U}(f;\boldsymbol{\theta})$ represents a frequency-domain GW
waveform in the absence of millilensing with $\boldsymbol{\theta}$
corresponding to unlensed BBH source parameters,
$F(f,\boldsymbol{\theta}_{L})$ is an amplification function dependent on
parameters associated with millilensing
$\boldsymbol{\theta_{L}}=(\boldsymbol{\mu},\boldsymbol{t},\boldsymbol{n})$ and
defined as:

$\displaystyle
F(f,\boldsymbol{\mu},\boldsymbol{t},\boldsymbol{n})=\sum_{j=1}^{K_{\mathrm{max}}}\left|\mu_{j}\right|^{1/2}\exp\left[2\pi
ift_{j}-i\pi n_{j}\right],$ (2)

where the expression is summed over the total number of millisignals up to a
chosen maximum number $K_{\mathrm{max}}$ and ($\mu_{j}$, $t_{j}$, $n_{j}$)
correspond to the lensing parameters of the $j^{\mathrm{th}}$ millisignal.
Such a summation corresponds to the geometrical optics approximation, which
applies when the GW wavelength $\lambda_{\mathrm{GW}}$ is longer than the
Schwarzschild radius $R_{S}$ corresponding to the lens mass:
$\lambda_{\mathrm{GW}}>R_{S}$ (Deguchi & Watson (1986); Takahashi & Nakamura
(2003)). For ground-based GW detectors, the geometrical optics approximation
can be applied down to masses $M_{Lz}\sim\mathcal{O}(10^{2})\,M_{\odot}$ (see
Appendix A for an explanation of the validity of geometrical optics
approximation). Throughout this work, we use geometrical units ($c=G=1$).

The notation introduced in Eq. (2) uses the conventional notation of lensing
magnification $\mu$. However, to overcome degeneracy between consecutive
millilensed signals, we use the effective luminosity distance notation which
relates to the true luminosity distance from the source $d_{L}$ and the
magnification of $j^{\mathrm{th}}$ lensed GW signal by:

$d^{\mathrm{eff}}_{j}=\frac{d_{L}}{\sqrt{\mu_{j}}}.$ (3)

This effective luminosity distance notation allows for a better separation of
individual millilensed signals and reduces degeneracy between the binary
orbital and millilensing parameters.

The individual millilensed GW signals must be time-ordered to avoid degeneracy
between their arrival times. To achieve this, we define the time delay between
each consecutive signal such that the arrival time of each signal increases
with $j$, that is, $t_{j+1}>t_{j}$. We can describe the millilensed GW signals
with respect to the first signal, which arrives earliest, by choosing
$t_{1}=0$. Using the effective luminosity distance notation introduced in Eq.
(3), we can express the amplification factor

$\begin{split}&F(f,d^{\mathrm{eff}}_{j},t_{j},n_{j})\\\
&=d_{L}\left(\frac{1}{d_{1}^{\mathrm{eff}}}\exp[-i\pi
n_{1}]+\frac{1}{d^{\mathrm{eff}}_{2}}\exp\left[2\pi ft_{2}-i\pi
n_{2}\right]+\cdots\right)\end{split}$ (4)

where the first signal (arriving the earliest) is described by its Morse
factor $n_{1}$ and the consecutive millilensed signals ($j\geq 2$) are
described by their effective luminosity distance $d^{\mathrm{eff}}_{j}$
(magnification) and time delay $t_{j}$ w.r.t the arrival time and the
luminosity distance of the first signal.

Consequently, the millilensed GW signal can be expressed as

$\begin{split}\tilde{h}_{L}(f)&=\sum_{j=1}^{K_{\mathrm{max}}}\frac{d_{L}}{d^{\mathrm{eff}}_{j}}\exp\left[2\pi
ift_{j}-i\pi
n_{j}\right]\tilde{h}_{U}(f)(f;\boldsymbol{\theta},d_{L},t_{\mathrm{col}})\\\
&=\sum_{j=1}^{K_{\mathrm{max}}}\exp\left[-i\pi
n_{j}\right]\tilde{h}_{U}(f)^{\mathrm{eff}}(f;\boldsymbol{\theta},d_{L},t_{\mathrm{col}},d^{\mathrm{eff}}_{\mathrm{j}},t_{\mathrm{j}})\end{split}$
(5)

where in the second line we have defined an effective GW signal dependent on
the source parameters, as well as the millilensing parameters
$({d^{\mathrm{eff}}_{j}},{t}_{j})$ of $j^{\mathrm{th}}$ component signal.

The resultant millilensed GW signal, which is a sum of $K$ individual
millilensed GW signals, can hence be expressed in terms of the amplification
function as

$\tilde{h}_{L}(f;\boldsymbol{\theta},\boldsymbol{\theta}^{L}_{K_{\mathrm{max}}},K_{\mathrm{max}})=F(f;\boldsymbol{\theta}^{L}_{K_{\mathrm{max}}},K_{\mathrm{max}})\cdot\tilde{h}(f;\boldsymbol{\theta})$
(6)

where $\boldsymbol{\theta}^{L}_{K_{\mathrm{max}}}$ represents the lensing
parameters of the sum of all $K_{\mathrm{max}}$ millilensed signals.

###  2.4 Problem statement in Bayesian framework

To perform parameter estimation and model selection on real GW data, we adopt
a Bayesian framework. Our goal is to obtain a result with a specified number
of millilensed GW signals, $K$, while considering varying dimensionality of
the models due to the changing number of signals. To address this problem, we
use joint Bayesian model selection, where the posterior distribution is
defined on a union of subspaces with different dimensions, each corresponding
to a model with a fixed number of signals. Assuming a countable set of
$K_{\mathrm{max}}$ models with $k$ models, we define the total union of
parameters as the BBH parameters $\boldsymbol{\theta}$ and lensing parameters
$\boldsymbol{\theta}^{L}_{K_{\mathrm{max}}}$ with $K_{\mathrm{max}}$ being the
maximum number of allowed signals specified by the user in each parameter
estimation run. Our aim is to estimate the unknowns $k$ and
$\boldsymbol{\theta}^{L}_{k}$, where $k\in{1,...,K_{\mathrm{max}}}$, given the
data set $\boldsymbol{d}$. We use the nested sampling algorithm (Skilling
(2006)) to draw these unknowns from corresponding prior distributions.

We express the joint distribution of all variables as

$p(k,\boldsymbol{\theta},\boldsymbol{\theta}^{L}_{k},\boldsymbol{d})=p(\boldsymbol{d}|\boldsymbol{\theta},\boldsymbol{\theta}^{L}_{k},k)p(\boldsymbol{\theta})p(\boldsymbol{\theta}^{L}_{k}|k)p(k)$
(7)

where $p(\boldsymbol{d}|\boldsymbol{\theta},\boldsymbol{\theta}^{L}_{k},k)$ is
the likelihood, $p(\boldsymbol{\theta}^{L}_{k}|k)$ is the prior distribution
of lensing parameters, $p(\boldsymbol{\theta})$ is the prior distribution of
the BBH source parameters and $p(k)$ is the prior distribution of the number
of millilensed signals.

We can express the marginal likelihood as follows (for detailed derivation,
see Appendix (B)):

$p(\boldsymbol{d})=\sum_{k=1}^{K_{\max}}\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\boldsymbol{\theta}^{L}_{k},k\right)p(\boldsymbol{\theta})p\left(\boldsymbol{\theta}^{L}_{k}\mid
k\right)p(k)d\boldsymbol{\theta}d\boldsymbol{\theta}^{L}_{k}.$ (8)

We can obtain the joint posterior distribution
$p(\boldsymbol{\theta}^{L}_{k},k|\boldsymbol{d})$ and the posterior
distribution for the number of images can be expressed as

$p(k\mid\boldsymbol{d})=\frac{\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\boldsymbol{\theta}^{L}_{k},k\right)p(\boldsymbol{\theta})p(\boldsymbol{\theta}^{L}_{k}\mid
k)p(k)d\boldsymbol{\theta}d\boldsymbol{\theta}^{L}_{k}}{\sum_{k=1}^{K_{\max}}\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\boldsymbol{\theta}^{L}_{k},k\right)p(\boldsymbol{\theta})p(\boldsymbol{\theta}^{L}_{k}\mid
k)p(k)d\boldsymbol{\theta}d\boldsymbol{\theta}^{L}_{k}}.$ (9)

##  3 Model-independent simulations

###  3.1 Simulations

We begin by simulating a lensed GW waveform with varying numbers of
millisignals using the lenstronomy package to model the physical setup of the
system. The simulated signal is injected with BBH parameters from GW190408 and
added to the detector Gaussian noise of the network of three ground-based
detectors (LIGO Livingston, LIGO Hanford, and Virgo) using the bilby Bayesian
inference library. We perform parameter estimation of both the source and
millilensing parameters using the phenomenological waveform approximant
IMRPhenomPv2 and the nested sampling technique with the dynesty sampler. The
priors used in the simulation are listed in Table 2, with the time delay and
magnification specified with continuous uniform priors and the Morse factor
and number of images having discrete uniform priors. We simulate three
injections with different signal-to-noise ratios (SNRs): 20, 30, and 50.

Table 1: Injected parameters for the SNR$=20$ parameter estimation run. BBH source parameters’ values correspond to the median values of GW190408 (Abbott et al., 2021c), millilensing parameters are obtained from simulating a lensing system with lenstronomy package (Birrer & Amara, 2018; Birrer et al., 2021) BBH source parameter | Value
---|---
Mass 1, $m_{1}$ | $31.6\,M_{\odot}$
Mass 2, $m_{2}$ | $23.7\,M_{\odot}$
Luminosity distance $d_{L}$ | 1598 Mpc
Dimensionless Spin 1, $a_{1}$ | 0.35
Dimensionless Spin 2, $a_{2}$ | 0.36
Tilt Angle 1, $\theta_{1}$ | 1.7 rad
Tilt Angle 2, $\theta_{2}$ | 1.6 rad
Right Ascension, $\mathrm{RA}$ | 6.088 rad
Declination, $\delta$ | 0.92 rad
Polarization, $\psi$ | 3\. 16 rad
Inclination, $\theta_{\mathrm{jn}}$ | 0.73 rad
Azimuthal Angle of $\vec{L}$, $\theta_{\mathrm{jl}}$ | 2.9 rad
Azimuthal Angle Difference, $\phi_{12}$ | 3.1 rad
Phase $\phi$ | 3.134 rad
Millilensing parameter | Value
Effective luminosity distance $d^{\mathrm{eff}}_{1}$ | 1598 Mpc
Effective luminosity distance $d^{\mathrm{eff}}_{2}$ | 1577 Mpc
Effective luminosity distance $d^{\mathrm{eff}}_{3}$ | 2570 Mpc
Effective luminosity distance $d^{\mathrm{eff}}_{4}$ | 2758 Mpc
Time delay $t_{2}$ | 0.0066 s
Time delay $t_{3}$ | 0.0467 s
Time delay $t_{4}$ | 0.0512 s
Morse phase $n_{1}$ | 0
Morse phase $n_{2}$ | 0
Morse phase $n_{3}$ | 0.5
Morse phase $n_{4}$ | 0.5
Number of millisignals $K$ | 4
Table 2: Prior distributions used for injection runs Parameter | Prior distribution
---|---
Luminosity distance $d_{L}$ | Uniform $\mathcal{U}$(50 Mpc, 20000 Mpc)
Effective luminosity distance $d^{\mathrm{eff}}$ | Uniform $\mathcal{U}$(50 Mpc, 20000 Mpc)
Time delay $dt$ | Uniform $\mathcal{U}$($10^{-3}$s, $10^{-1}$s)
Morse factor $n$ | Discrete Uniform $\mathcal{U}\\{0,0.5,1\\}$
Number of millisignals $K$ | Discrete Uniform $\mathcal{U}\\{1,2,3,4,5,6\\}$

###  3.2 Lens configuration

We start with a strongly lensed GW signal by a galaxy resulting in two GW
signals lensed at different deflection angles, and each of the two strongly
lensed signals passes near a millilens located within the same plane as the
strong lens, which experiences a gravitational shearing effect due to the
presence of the strong lens galaxy. The modelling is divided into two steps:
first, solving the lens equation for the strong lens to obtain the shear
components ($\gamma_{1},\gamma_{2}$), convergence components
($\kappa_{1},\kappa_{2}$) and magnifications ($\mu_{1}^{S},\mu_{2}^{S}$) of
the two strongly lensed GW signals; secondly, using the shear values from a
strong lens as a macro-environment input into millilens models to obtain the
magnifications, time delays and Morse factors for the resultant millisignals
which are then estimated in injection run parameter estimation. In this
example, we assume that one of the strongly lensed signals is further split
into four millilensed signals, but the choice of the number of millisignals is
arbitrary and can take different integer values depending on the system setup,
convergence, and shear values.

The strong lens is modelled as a galaxy with a singular isothermal ellipsoid
(SIE) mass distribution located at redshift $z_{\mathrm{lens}}=0.3$, with
eccentricity components $e_{1}=0.05$ and $e_{2}=0$ and an Einstein radius of
$\theta_{E}=1$. The millilens is modelled as an SIS mass distribution, and the
source is a BBH emitting GWs at typical redshift $z_{\mathrm{source}}=1.5$
(Wierda et al. (2021b)). A flat $\Lambda$-CDM cosmological model is assumed
throughout the simulation with $H_{0}=70;\mathrm{kms}^{-1}\mathrm{Mpc}^{-1}$
and $\Omega_{\mathrm{M}}=0.3$. Results of the injection run are presented in
the following section.

###  3.3 Example: 4 millisignals

The case reported here illustrates a GW signal lensed into 4 millisignals
injected into detector noise. Fig. 4 shows the posterior distribution for
parameter $K$, corresponding to the number of millilensed GW signals
recovered. All three runs with three different SNR values lead to consistent
results, but for the clarity of the figure, only results from the injection
run with SNR 20 are shown. As can be seen from the plot, the number of
millisignals has been recovered in agreement with the injected value of $K=4$.

Fig. 5(a) presents a corner plot of the effective luminosity distance
parameters obtained for injection run with SNR 20. It can be seen from the
figure that the injected values (orange lines) are well recovered in the
posterior distributions for each parameter. Consistent results for higher SNR
values were obtained. Similarly, the posterior distribution of time delays of
the consecutive millilensed component signals are presented in Fig. 5(b). The
injected values are recovered within $1\sigma$ region of the posterior
distributions for $t_{2}$ and $t_{4}$, and within $2\sigma$ for $t_{3}$. As
tested with higher SNR runs, the accuracy of the recovered parameters
increases with SNR. The recovered BBH source parameters are shown in Fig 6.

Figure 4: Number of millilensed signals, $K$, recovered from the injection run
with uniform prior. The posterior distribution (shown in purple) is in
agreement with the injected value $K=4$ (orange line).

(a)

(b)

Figure 5: Results recovered from parameter estimation of the injected
millilensed GW with an SNR 20 and GW190408 BBH parameters: (a) effective
luminosity distances, (b) time delays of component millilensed signals w.r.t.
the earliest arriving signal ($t_{1}=0)$. The millilens model is an SIS
embedded into a strong lensing galaxy, which we model with the SIE lens model.
The colours of the two-dimensional corner plots indicate 1$\sigma$, 2$\sigma$
and 3$\sigma$ credible regions. The orange lines show injected values for each
parameter and the dashed lines in individual histogram plots correspond to
1$\sigma$ credible intervals. The effective luminosity distances and time
delays are recovered with good accuracy and peak at the appropriate injected
values. Figure 6: Black hole source parameters posterior distributions
recovered from parameter estimation of the injected millilensed GW: chirp mass
$\mathcal{M}$, dimensionless spin parameters $a_{1}$, $a_{2}$, tilt angles
$\theta_{1}$, $\theta_{2}$. The orange lines represent the injected values for
each parameter. The posterior distributions recovered are consistent with
injected values.

##  4 lens mapping

The methodology developed so far aims to provide a physically realistic
picture of the lensing system, accounting for the fact that the lens located
within a galaxy can be gravitationally affected by it. Having performed
parameter estimation with the phenomenological approach, the results obtained
can be mapped to specific lens models. We present an example mapping to an SIS
model below. We choose the SIS model for its simplicity, which allows for
analytical mapping from the phenomenological to SIS model parameters and vice
versa. We use it as a testing example to demonstrate the accuracy of the
method, noting that the method is not sensitive to physical assumptions and
can be applied to more complex models. However, the drawback associated with
the use of the SIS model is that it is a simplified model, which can not be
applied to physically generic mass distributions if we want to take into
account other physical effects, such as galactic shear.

###  4.1 Singular Isothermal Sphere

In the geometrical optics limit, the time delay and magnification of two
millilensed signals can be related to SIS lens parameters ($y$, $M_{Lz}$) by

$\displaystyle t_{d}=8M_{Lz}y$ (10) $\displaystyle\mu_{\pm}=\pm 1+1/y$

where $t_{d}$ is the time delay between two signals and $\mu_{\pm}$ are the
magnifications of the two lensed signals which can be related to effective
luminosity distances following Eq. (3), all expressed in $c=G=1$ units
(Takahashi & Nakamura (2003)).

The first step of the mapping is to extract the necessary data from the multi-
signal millilensing analysis. If we consider $y\leq 1$, the SIS model predicts
that double-lensed GW signals are formed. Restricting ourselves to this case,
we perform a millilensing injection run injecting two millilensed GW signals
into detector noise, but allowing their number $K$ to take on values up to and
including 6. From the analysis result, we select the posteriors of
$(d^{\mathrm{eff}}_{2},t_{2})$ corresponding to $K=2$ case in order to map it
to the predicted two-signal SIS model. Then, we construct the likelihood
marginalized over all BBH parameters $\theta$ except for
$(d^{\mathrm{eff}}_{2},t_{2})$:

$p(\boldsymbol{d}|d^{\mathrm{eff}}_{2},t_{2})=\int
p(\boldsymbol{d}|d^{\mathrm{eff}}_{2},t_{2},\theta)p(\theta)d\theta.$ (11)

Knowing the relation between two sets of parameters (Eq. (10)), we assume that
the likelihood can be expressed as

$p(\boldsymbol{d}|d^{\mathrm{eff}}_{2},t_{2})=p(\boldsymbol{d}|d^{\mathrm{eff}}_{2}(y),t_{2}(M_{Lz},y))=p(\boldsymbol{d}|y,M_{Lz}).$
(12)

Then, using the likelihood, as well as assuming uniform prior distributions
for ($y$, $M_{Lz}$), we perform a nested sampling over the two SIS parameters
($y$, $M_{Lz}$) in order to obtain corresponding posterior distributions and
evidence for the SIS model. For the purpose of implementing the mapped
likelihood into the Bayesian inference library bilby, we used the likelihood
ratio to perform the nested sampling with the results presented in the
following section.

Figure 7: Posterior distribution of the redshifted lens mass $M_{Lz}$ (top)
and source position parameter $y$ (bottom), obtained from mapping
phenomenological results to an SIS lens model. The mapping was performed: i)
analytically (histogram in yellow), ii) by nested sampling algorithm
(histogram in purple). The results confirm that we are able to map the results
from observables to lensing parameters with nested sampling.

Following Eq. (10), it is also possible to write down the inverse relations

$\displaystyle y=\frac{\mu_{\mathrm{rel}}-1}{\mu_{\mathrm{rel}}+1}$ (13)
$\displaystyle
M_{Lz}=\frac{t_{d}}{8}\frac{\mu_{\mathrm{rel}}+1}{\mu_{\mathrm{rel}}-1}$

where $\mu_{\mathrm{rel}}$ is the relative magnification $\mu_{+}/\mu_{-}$
between two lensed GW signals. For validation of the results obtained from
nested sampling, we construct posterior distributions for the SIS parameters
($y,M_{Lz}$) by taking samples from the posterior distributions
$p(t_{d}|\boldsymbol{d})$, $p(d^{\mathrm{eff}}_{2}|\boldsymbol{d})$ and using
the analytical relations Eq. (13).

###  4.2 Results

Figure 7 represents the posterior distributions of SIS parameters $M_{Lz}$ and
$y$, obtained from mapping phenomenological results. As described above, we
used two methods to obtain the posterior distributions of the SIS parameters.
Firstly, we performed a nested sampling algorithm (plots shown in purple in
Fig. 7). Secondly, to validate the results, we also used analytical relations
from Eq. (13) to construct the posteriors for $M_{Lz}$ and $y$ (plots in
yellow in Fig. 7). The prior ranges of ($M_{Lz},y$) used in the two methods
were different, hence the shape and position of the peak of the two
distributions are not exactly the same (Fig. 7). Nevertheless, the
distributions obtained with the two methods lead to consistent results.

##  5 Discussion and conclusion

###  5.1 Discussion

In this work, we have presented a new phenomenological approach to GW
millilensing analysis, applicable to astrophysical lenses with masses in the
range $10^{3}-10^{6}\,M_{\odot}$. The novel approach, unlike currently most
widely used models, parameterizes the lensed GW signals, which not only
enables to include the physical effects due to the presence of other massive
objects in the vicinity of the gravitational lens but also provides an
efficient lens model selection tool. Having tested the feasibility of the
approach we applied the method to simulated injections of millilensed GW
signals into detector noise. The parameter estimation of the lensed BBH
parameters recovered results that are in agreement with injected values.

The question raised by this study aims to provide a physically realistic
description of an astrophysical lens within a galaxy. However, the
generalizability of the results is subject to the limitation of assuming a
single-lens system. More broadly, the study should be repeated addressing
systems with multiple lenses located in the vicinity of each other (within the
lens plane) or located along the line of sight. This would be a fruitful area
for further work, providing a more generic description of gravitationally
millilensed systems in different physical configurations.

Moreover, further research should be undertaken to investigate the distinction
of millilensing effects on the gravitational waveform from other physical
processes, such as spin precession or non-GR effects, which may also lead to
frequency-dependent beating patterns in the waveform, mimicking gravitational
lensing. It would be of crucial importance to distinguish between those
effects once a GW signal with a potentially lensed waveform is detected.

The methodology can be applied to a range of problems within theoretical
physics, astrophysics and cosmology. A potential interesting example
application is the problem of dark matter subhalos. Despite many observations
supporting the presence of dark matter, the nature of dark matter remains one
of the key open questions within the field (Zackrisson & Riehm (2010); Ellis
(2010); Bertone & Hooper (2016); Bertone & Tait (2018)). To date, a number of
dark matter models have been proposed and some of them predict dark matter
halos to be formed hierarchically from smaller subhalos which can further be
formed from even smaller subhalos (Moore et al. (1999); Metcalf & Madau
(2001); Liao et al. (2018); Dai et al. (2018)). However, current dark matter
models show discrepancies at the smallest scales, often referred to as the
substructure crisis (Metcalf & Madau (2001); Somerville (2002); Moore et al.
(2006)). In particular, some models predict dark matter subhalos cannot be
present below a certain scale (Kravtsov (2010); Oguri & Takahashi (2020)).
Therefore, in order to test the most feasible dark matter models, it is
necessary to study dark matter subhalos down to the smallest scales. Moreover,
if dark matter subhalos with masses of order $10^{3}-10^{6}\;M_{\odot}$ exist
within galaxies as predicted by some models, they could potentially act as
millilenses, directly influencing distant GW signals and producing lensing
distortions. The probe of dark matter subhalos with gravitational millilensing
has been proposed within EM studies of lensed quasars (Wambsganss & Paczynski,
1992; Mao & Schneider, 1997; Metcalf & Madau, 2001), however, due to telescope
resolution and propagation effects, EM observations can be subject to
uncertainties. GW millilensing is subject to different systematics and could
potentially be used as a direct measurement of the unknown dark matter subhalo
mass function. Therefore, the dark matter substructures problem is a promising
direction to study with the millilensing approach developed, which could
become an alternative probe of the nature of dark matter at small scales.
However, more investigation is needed to study the detailed capabilities of GW
millilensing as a probe of dark matter subhalos.

###  5.2 Conclusions

The present study was designed to develop a phenomenological approach to
gravitational millilensing studies which accounts for effects not included in
currently used lens models. The simulations confirmed the feasibility of the
method and the results support the idea that it is possible to study millilens
configurations embedded in macro systems with mutual gravitational
interactions. These findings suggest that in general, we can analyse non-
isolated and non-symmetric gravitational lenses. Furthermore, results obtained
from millilensing analysis can be mapped to existing lens models, providing a
useful tool to distinguish between the most feasible models. These results add
to the rapidly expanding field of gravitational-wave lensing and will prove
useful in developing analysis tools for observational GW data with lensed GW
detections predicted to take place in the coming years. The major limitation
of this study is the geometrical optics approximation which limits the
millilens mass range considered and the target lens population.
Notwithstanding the relatively limited lens sample, this work offers valuable
insights into the studies of gravitational millilensing analysis and can be
further developed into population studies of potential millilens candidates.
Non-observation of millilensing of predicted lenses at typical redshifts could
also shed light on lens populations. The approach can also be expanded into
generic millilensing studies of multiply lensed GW signals. Moreover, the
millilensing framework could be applied to existing problems in astrophysics
and cosmology, such as studies of the dark matter subhalos and primordial
black holes with masses expected to lie within the corresponding millilensing
mass range. More work will be needed, however, to study the detailed science
case.

## Acknowledgements

The work is partially supported by grants from the Research Grants Council of
the Hong Kong (CUHK 14306218), The Research Foundation – Flanders (G086722N)
and KU Leuven (STG/21/061). The analysed data and the corresponding power
spectral densities are publicly available at the online Gravitational-Wave
Open Science Center (Abbott et al., 2021e). The authors are grateful for
computational resources provided by the CIT cluster of the LIGO Laboratory and
supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.
This material is based upon work supported by NSF’s LIGO Laboratory which is a
major facility fully funded by the National Science Foundation. This
manuscript has LIGO DCC number P2200365.

## Data Availability

This research has made use of data or software obtained from the Gravitational
Wave Open Science Center (gwosc.org), a service of LIGO Laboratory, the LIGO
Scientific Collaboration, the Virgo Collaboration, and KAGRA. LIGO Laboratory
and Advanced LIGO are funded by the United States National Science Foundation
(NSF) as well as the Science and Technology Facilities Council (STFC) of the
United Kingdom, the Max-Planck-Society (MPS), and the State of
Niedersachsen/Germany for support of the construction of Advanced LIGO and
construction and operation of the GEO600 detector. Additional support for
Advanced LIGO was provided by the Australian Research Council. Virgo is
funded, through the European Gravitational Observatory (EGO), by the French
Centre National de Recherche Scientifique (CNRS), the Italian Istituto
Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions
by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan,
Monaco, Poland, Portugal, Spain. KAGRA is supported by Ministry of Education,
Culture, Sports, Science and Technology (MEXT), Japan Society for the
Promotion of Science (JSPS) in Japan; National Research Foundation (NRF) and
Ministry of Science and ICT (MSIT) in Korea; Academia Sinica (AS) and National
Science and Technology Council (NSTC) in Taiwan.

## References

  * Aasi et al. (2015) Aasi J., et al., 2015, Classical and quantum gravity, 32, 074001
  * Abbott et al. (2016a) Abbott B. P., et al., 2016a, Physical Review X, 6, 041014
  * Abbott et al. (2016b) Abbott B. P., et al., 2016b, Physical review letters, 116, 061102
  * Abbott et al. (2019) Abbott B., et al., 2019, Physical Review X, 9, 031040
  * Abbott et al. (2021a) Abbott R., et al., 2021a, arXiv preprint arXiv:2111.03606
  * Abbott et al. (2021c) Abbott R., et al., 2021c, Open data from the first and second observing runs of Advanced LIGO and Advanced Virgo
  * Abbott et al. (2021b) Abbott R., et al., 2021b, Search for lensing signatures in the gravitational-wave observations from the first half of LIGO-Virgo’s third observing run (arXiv:2105.06384)
  * Abbott et al. (2021d) Abbott R., et al., 2021d, Physical Review X, 11, 021053
  * Abbott et al. (2021e) Abbott R., et al., 2021e, SoftwareX, 13, 100658
  * Acernese et al. (2014) Acernese F. a., et al., 2014, Classical and Quantum Gravity, 32, 024001
  * Akutsu et al. (2019) Akutsu T., Ando M., Arai K., Arai Y., et al., 2019, Nature Astronomy, 3, 35
  * Baker & Trodden (2016) Baker T., Trodden M., 2016, Physical Review D, 95
  * Basak et al. (2022a) Basak S., Sharma A. K., Kapadia S. J., Ajith P., 2022a, http://arxiv.org/abs/2205.00022
  * Basak et al. (2022b) Basak S., Ganguly A., Haris K., Kapadia S., Mehta A. K., Ajith P., 2022b, The Astrophysical Journal Letters, 926, L28
  * Bertone & Hooper (2016) Bertone G., Hooper D., 2016, Reviews of Modern Physics, 90
  * Bertone & Tait (2018) Bertone G., Tait T. M., 2018, arXiv preprint arXiv:1810.01668
  * Birrer & Amara (2018) Birrer S., Amara A., 2018, Physics of the Dark Universe, 22, 189
  * Birrer et al. (2021) Birrer S., et al., 2021, Journal of Open Source Software, 6, 3283
  * Bolejko et al. (2012) Bolejko K., Clarkson C., Maartens R., Bacon D., Meures N., Beynon E., 2012, Physical Review Letters, 110
  * Bulashenko & Ubach (2021) Bulashenko O., Ubach H., 2021, Journal of Cosmology and Astroparticle Physics, 2022
  * Çalışkan et al. (2022) Çalışkan M., Ezquiaga J. M., Hannuksela O. A., Holz D. E., 2022, arXiv preprint arXiv:2201.04619
  * Cao et al. (2014) Cao Z., Li L.-F., Wang Y., 2014, Physical Review D, 90, 62003
  * Cao et al. (2022) Cao S., Qi J., Cao Z., Biesiada M., Cheng W., Zhu Z.-H., 2022, Astronomy and Astrophysics, 659
  * Cheung et al. (2020) Cheung M. H. Y., Gais J., Hannuksela O. A., Li T. G. F., 2020, Monthly Notices of the Royal Astronomical Society, 503, 3326
  * Christian et al. (2018) Christian P., Vitale S., Loeb A., 2018, Physical Review D, 98
  * Chung & Li (2021) Chung A. K. W., Li T. G., 2021, Physical Review D, 104, 124060
  * Collett & Bacon (2017) Collett T. E., Bacon D., 2017, Physical Review Letters, 118, 091101
  * Congdon & Keeton (2018) Congdon A. B., Keeton C. R., 2018, Principles of Gravitational Lensing: Light Deflection as a Probe of Astrophysics and Cosmology. Springer Cham, doi:10.1007/978-3-030-02122-1
  * Cremonese et al. (2021) Cremonese P., Ezquiaga J. M., Salzano V., 2021, Physical Review D, 104, 023503
  * Dai & Venumadhav (2017) Dai L., Venumadhav T., 2017, arxiv preprint http://arxiv.org/abs/1702.04724
  * Dai et al. (2016) Dai L., Venumadhav T., Sigurdson K., 2016, Physical Review D, 95
  * Dai et al. (2018) Dai L., Li S.-S., Zackay B., Mao S., Lu Y., 2018, Physical Review D, 98
  * Dai et al. (2020) Dai L., Zackay B., Venumadhav T., Roulet J., Zaldarriaga M., 2020, arxiv preprint http://arxiv.org/abs/2007.12709
  * Deguchi & Watson (1986) Deguchi S., Watson W. D., 1986, Phys. Rev. D, 34, 1708
  * Diego (2019) Diego J. M., 2019, Physical Review D, 101
  * Diego et al. (2019) Diego J. M., Hannuksela O. A., Kelly P. L., Pagano G., Broadhurst T., Kim K., Li T. G. F., Smoot G. F., 2019, Astronomy & Astrophysics, 627, A130
  * Dyson et al. (1920) Dyson F. W., Eddington A. S., Davidson C., 1920, Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 220, 291
  * Ellis (2010) Ellis R. S., 2010, Philosophical Transactions of the Royal Society of London Series A, 368, 967
  * Ezquiaga et al. (2021) Ezquiaga J. M., Holz D. E., Hu W., Lagos M., Wald R. M., 2021, Physical Review D, 103, 064047
  * Fan et al. (2017) Fan X. L., Liao K., Biesiada M., Piórkowska-Kurpas A., Zhu Z. H., 2017, Physical Review Letters, 118, 091102
  * Finke et al. (2021) Finke A., Foffa S., Iacovelli F., Maggiore M., Mancarella M., 2021, Physical Review D, 104
  * Goyal et al. (2020) Goyal S., Haris K., Mehta A. K., Ajith P., 2020, Physical Review D, 103
  * Hannuksela et al. (2019) Hannuksela O. A., Haris K., Ng K. K. Y., Kumar S., Mehta A. K., Keitel D., Li T. G. F., Ajith P., 2019, The Astrophysical Journal, 874, L2
  * Hannuksela et al. (2020) Hannuksela O. A., Collett T. E., Çalışkan M., Li T. G. F., 2020, Monthly Notices of the Royal Astronomical Society, 498, 3395
  * Haris et al. (2018) Haris K., Mehta A. K., Kumar S., Venumadhav T., Ajith P., 2018, http://arxiv.org/abs/1807.07062
  * Iacovelli et al. (2022) Iacovelli F., Finke A., Foffa S., Maggiore M., Mancarella M., 2022, arxiv preprint http://arxiv.org/abs/2203.09237
  * Itoh et al. (2009) Itoh Y., Futamase T., Hattori M., 2009, Physical Review D - Particles, Fields, Gravitation and Cosmology, 80
  * Janquart et al. (2021a) Janquart J., Hannuksela O. A., Haris K., Van Den Broeck C., 2021a, Monthly Notices of the Royal Astronomical Society, 506, 5430
  * Janquart et al. (2021b) Janquart J., Seo E., Hannuksela O. A., Li T. G. F., Broeck C. V. D., 2021b, The Astrophysical Journal Letters, 923, L1
  * Janquart et al. (2022) Janquart J., Hannuksela O. A., Haris K., Broeck C. V. D., 2022, arXiv preprint arXiv:2203.06444
  * Jung & Shin (2019) Jung S., Shin C. S., 2019, Phys. Rev. Lett., 122, 041103
  * Kim et al. (2020) Kim K., Lee J., Yuen R. S. H., Hannuksela O. A., Li T. G. F., 2020, The Astrophysical Journal, 915, 119
  * Kim et al. (2022) Kim K., Lee J., Hannuksela O. A., Li T. G. F., 2022, arxiv preprint https://arxiv.org/abs/2206.08234v2
  * Kravtsov (2010) Kravtsov A., 2010, Advances in Astronomy, 2010
  * Lai et al. (2018) Lai K.-H., Hannuksela O. A., Herrera-Martín A., Diego J. M., Broadhurst T., Li T. G. F., 2018, Physical Review D, 98
  * Li et al. (2018) Li S.-S., Mao S., Zhao Y., Lu Y., 2018, Monthly Notices of the Royal Astronomical Society, 476, 2220
  * Li et al. (2019) Li A. K. Y., Lo R. K. L., Sachdev S., Chan C. L., Lin E. T., Li T. G. F., Weinstein A. J., 2019, arxiv preprint http://arxiv.org/abs/1904.06020
  * Liao et al. (2017) Liao K., Fan X. L., Ding X., Biesiada M., Zhu Z. H., 2017, Nature Communications, 8, 1
  * Liao et al. (2018) Liao K., Ding X., Biesiada M., Fan X.-L., Zhu Z.-H., 2018, The Astrophysical Journal, 867, 69
  * Liu et al. (2020) Liu X., Hernandez I. M., Creighton J., 2020, The Astrophysical Journal, 908, 97
  * Lo & Hernandez (2021) Lo R. K. L., Hernandez I. M., 2021, arxiv preprint http://arxiv.org/abs/2104.09339
  * Ma et al. (2022) Ma H., Lu Y., Guo X., Zhang S., Chu Q., 2022, Monthly Notices of the Royal Astronomical Society, 518, 6183
  * Magaña Hernandez (2022) Magaña Hernandez I., 2022, arXiv preprint arXiv:2211.01272
  * Mao & Schneider (1997) Mao S., Schneider P., 1997, Monthly Notices of the Royal Astronomical Society, 295, 587
  * Meena & Bagla (2020) Meena A. K., Bagla J. S., 2020, Monthly Notices of the Royal Astronomical Society, 492, 1127
  * Metcalf & Madau (2001) Metcalf R. B., Madau P., 2001, The Astrophysical Journal, 563, 9
  * Mishra et al. (2021) Mishra A., Meena A. K., More A., Bose S., Bagla J. S., 2021, Monthly Notices of the Royal Astronomical Society, 508, 4869
  * Moore et al. (1999) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, The Astrophysical Journal, 524, L19
  * Moore et al. (2006) Moore B., Diemand J., Madau P., Zemp M., Stadel J., 2006, Monthly Notices of the Royal Astronomical Society, 368, 563
  * More & More (2021) More A., More S., 2021, arxiv preprint http://arxiv.org/abs/2111.03091
  * Mukherjee et al. (2019) Mukherjee S., Wandelt B. D., Silk J., 2019, Physical Review D, 101
  * Mukherjee et al. (2020) Mukherjee S., Wandelt B. D., Silk J., 2020, Monthly Notices of the Royal Astronomical Society, 494, 1956
  * Nakamura (1998) Nakamura T. T., 1998, Phys. Rev. Lett., 80, 1138
  * Ng et al. (2017) Ng K. K. Y., Wong K. W. K., Li T. G. F., Broadhurst T., 2017, Physical Review D, 97
  * Oguri (2018) Oguri M., 2018, Monthly Notices of the Royal Astronomical Society, 480, 3842
  * Oguri & Takahashi (2020) Oguri M., Takahashi R., 2020, The Astrophysical Journal, 901, 58
  * Oguri & Takahashi (2022) Oguri M., Takahashi R., 2022, Physical Review D, 106
  * Pagano et al. (2020) Pagano G., Hannuksela O. A., Li T. G. F., 2020, Astronomy and Astrophysics, 643
  * Pang et al. (2020) Pang P. T. H., Hannuksela O. A., Dietrich T., Pagano G., Harry I. W., 2020, Monthly Notices of the Royal Astronomical Society, 495, 3740
  * Qiu et al. (2022) Qiu Y., Wang K., He J.-h., 2022, arxiv preprint http://arxiv.org/abs/2205.01682
  * Robertson et al. (2020) Robertson A., Smith G. P., Massey R., Eke V., Jauzac M., Bianconi M., Ryczanowski D., 2020, Monthly Notices of the Royal Astronomical Society, 495, 3727
  * Ryczanowski et al. (2020) Ryczanowski D., Smith G. P., Bianconi M., Massey R., Robertson A., Jauzac M., 2020, Monthly Notices of the Royal Astronomical Society, 495, 1666
  * Seo et al. (2021) Seo E., Hannuksela O. A., Li T. G. F., 2021, arxiv preprint http://arxiv.org/abs/2110.03308
  * Sereno et al. (2010) Sereno M., Sesana A., Bleuler A., Jetzer P., Volonteri M., Begelman M. C., 2010, Physical Review Letters, 105, 251101
  * Skilling (2006) Skilling J., 2006, Bayesian Analysis, 1, 833
  * Smith et al. (2017) Smith G. P., Jauzac M., Veitch J., Farr W. M., Massey R., Richard J., 2017, Monthly Notices of the Royal Astronomical Society, 475, 3823
  * Smith et al. (2018a) Smith G. P., et al., 2018a, Proceedings of the International Astronomical Union, 13, 98
  * Smith et al. (2018b) Smith G. P., et al., 2018b, Monthly Notices of the Royal Astronomical Society, 485, 5180
  * Smith et al. (2022) Smith G. P., et al., 2022, Monthly Notices of the Royal Astronomical Society
  * Somerville (2002) Somerville R. S., 2002, The Astrophysical Journal, 572, L23
  * Takahashi (2005) Takahashi R., 2005, The Astrophysical Journal, 644, 80
  * Takahashi & Nakamura (2003) Takahashi R., Nakamura T., 2003, The Astrophysical Journal, 595, 1039
  * Tambalo et al. (2022b) Tambalo G., Zumalacárregui M., Dai L., Cheung M. H.-Y., 2022b, https://arxiv.org/abs/2212.11960v1
  * Tambalo et al. (2022a) Tambalo G., Zumalacárregui M., Dai L., Cheung M. H.-Y., 2022a, https://arxiv.org/abs/2210.05658v1
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration et al., 2021, The Astrophysical Journal, 923, 14
  * Vijaykumar et al. (2022) Vijaykumar A., Mehta A. K., Ganguly A., 2022, http://arxiv.org/abs/2202.06334
  * Wambsganss & Paczynski (1992) Wambsganss J., Paczynski B., 1992, ApJ, 397, L1
  * Wang et al. (2021) Wang Y., Lo R. K. L., Li A. K. Y., Chen Y., 2021, Physical Review D, 103
  * Wempe et al. (2022) Wempe E., Koopmans L. V. E., Wierda A. R. A. C., Hannuksela O. A., Broeck C. v. d., 2022, arxiv preprint http://arxiv.org/abs/2204.08732
  * Wierda et al. (2021a) Wierda A. R. A. C., Wempe E., Hannuksela O. A., Koopmans L. V. E., Broeck C. V. D., 2021a, The Astrophysical Journal, 921, 154
  * Wierda et al. (2021b) Wierda A. R. A. C., Wempe E., Hannuksela O. A., Koopmans L. e. V. E., Van Den Broeck C., 2021b, Astrophys. J., 921, 154
  * Wright & Hendry (2021) Wright M., Hendry M., 2021, The Astrophysical Journal, 935, 68
  * Xu et al. (2021) Xu F., Ezquiaga J. M., Holz D. E., 2021, The Astrophysical Journal, 929, 9
  * Yeung et al. (2021) Yeung S. M. C., Cheung M. H. Y., Gais J. A. J., Hannuksela O. A., Li T. G. F., 2021, arxiv preprint http://arxiv.org/abs/2112.07635
  * Zackrisson & Riehm (2010) Zackrisson E., Riehm T., 2010, Advances in Astronomy, 2010

##  Appendix A Validity of geometrical optics approximation

Throughout this work, we have assumed the geometrical optics approximation. In
order to verify the validity of the approximation in the corresponding mass
range considered, we compute the match between the lensed GW in geometric and
wave optics for a PML with redshifted lens mass
$M_{Lz}\in[10,10^{3}]\,M_{\odot}$ and source position $y\in[0.01,1]$. We
define the match as

$\mathcal{M}(\tilde{h}^{L}_{\mathrm{wave}},\tilde{h}^{L}_{\mathrm{geo}}):=\frac{\left(\tilde{h}^{L}_{\mathrm{wave}},\tilde{h}^{L}_{\mathrm{geo}}\right)}{\sqrt{(\tilde{h}^{L}_{\mathrm{wave}},\tilde{h}^{L}_{\mathrm{wave}})(\tilde{h}^{L}_{\mathrm{geo}},\tilde{h}^{L}_{\mathrm{geo}})}}$
(14)

where the numerator is the noise-weighted inner product between the waveforms
in two approximations

$\left(\tilde{h}^{L}_{\mathrm{wave}},\tilde{h}^{L}_{\mathrm{geo}}\right)=4\;\Re\int_{f_{\mathrm{low}}}^{f_{\mathrm{high}}}\frac{\tilde{h}^{L}_{\mathrm{wave}}(f)\tilde{h}^{\ast
L}_{\mathrm{geo}}(f)}{S_{n}(f)}df$ (15)

with $S_{n}(f)$ corresponding to the power spectral density (PSD) of the noise
(in this case we use aLIGO design sensitivity PSD with frequency range [10,
2048] Hz), and the lensed GWs described as
$\tilde{h}^{L}_{\mathrm{wave}}(f)=F_{\mathrm{wave}}(f)\tilde{h}(f)$ for wave
optics effects and
$\tilde{h}^{L}_{\mathrm{geo}}(f)=F_{\mathrm{geo}}(f)\tilde{h}(f)$ in
geometrical optics approximation where the amplification functions are
described respectively by

$F_{\mathrm{wave}}(f)=\exp{\left[\frac{\pi
w}{4}+i\frac{w}{2}\left(\ln\left(\frac{w}{2}\right)-2\phi_{m}(y)\right)\right]}\times\Gamma\left(1-\frac{i}{2}w\right){}_{1}F_{1}\left(\frac{i}{2}w,1;\frac{i}{2}wy^{2}\right)$
(16)

where $w=8\pi M_{Lz}f;\phi_{m}(y)=\left(x_{m}-y\right)^{2}/2-\ln x_{m}$ with
$x_{m}=\left(y+\sqrt{y^{2}+4}\right)/2;M_{Lz}=$ $M_{L}\left(1+z_{L}\right)$ is
the redshifted lens mass and ${}_{1}F_{1}$ is the confluent hypergeometric
function (Takahashi & Nakamura (2003)); in the geometrical optics limit, the
amplification function is

$F_{\mathrm{geo}}(f)=\left|\mu_{+}\right|^{1/2}-i\left|\mu_{-}\right|^{1/2}e^{2\pi
if\Delta t_{d}}$ (17)

where the magnifications of the two lensed images are
$\mu_{\pm}=1/2\pm\left(y^{2}+2\right)/\left(2y\sqrt{y^{2}+4}\right)$ and the
time delay between two lensed signals is

$\Delta
t_{d}=4M_{Lz}\left[y\sqrt{y^{2}+4}/2+\ln\left(\left(\sqrt{y^{2}+4}+y\right)/\left(\sqrt{y^{2}+4}-y\right)\right)\right].$

From the match between the two waveforms, we can conclude that the geometric
and wave optics effects match over 95% for lens masses of order 200
$M_{\odot}$ (see Fig. 8). Therefore, we can conclude that geometrical optics
approximation is valid for the mass range considered in the millilensing
analysis presented in this paper: $M\in[10^{3},10^{6}]M_{\odot}$.

Figure 8: Contour plot representing the match between the plus polarisation of
a lensed GW in geometrical optics approximation and the one of a lensed GW
accounting for wave optics effects, plotted as a function of the relative
source position $y$ and redshifted lens mass $M_{\mathrm{Lz}}$. The lens is
assumed to be a PML.

##  Appendix B Derivations

The following section presents derivation for the marginal likelihood of the
millilensing methodology.

The marginal likelihood of the millilensing model with an unknown number of
images takes the general form

$\begin{split}p(\boldsymbol{d})&=\sum_{k=1}^{K_{\max}}\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\boldsymbol{\theta}^{L}_{k},k\right)p(\boldsymbol{\theta})p\left(\boldsymbol{\theta}^{L}_{k}\mid
k\right)p\left(\boldsymbol{\theta}^{L}_{-k}\mid
k\right)p(k)\,d\boldsymbol{\theta}d\boldsymbol{\theta}^{L}_{k}d\boldsymbol{\theta}^{L}_{-k}\\\
&=\sum_{k=1}^{K_{\max}}\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\boldsymbol{\theta}^{L}_{k},k\right)p(\boldsymbol{\theta})p\left(\boldsymbol{\theta}^{L}_{k},\boldsymbol{\theta}^{L}_{-k}\mid
k\right)p(k)\,d\boldsymbol{\theta}d\boldsymbol{\theta}^{L}_{k}d\boldsymbol{\theta}^{L}_{-k}\\\
&=\sum_{k=1}^{K_{\max}}\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\boldsymbol{\theta}^{L}_{k},k\right)p(\boldsymbol{\theta})p(\boldsymbol{\theta}^{L}_{k}\mid
k)p(k)\,d\boldsymbol{\theta}d\boldsymbol{\theta}^{L}_{k}\end{split}$ (18)

where $\boldsymbol{\theta}^{L}_{k}$ corresponds to the lensing parameters of
$k$ millilensed GW signals, $\boldsymbol{\theta}^{L}_{-k}$ represent
millilensing parameters of those millilensed signals for which the actual
number of millilensed signals formed is greater than $k$. The two terms are
assumed to be independent of each other in the prior distribution. We provide
an illustrative example in the following section.

## Two-signal example

Let us consider an example with $K_{\mathrm{max}}=2$. Following equation (18),
we can write

$\begin{split}p(\boldsymbol{d})&=\sum_{k=1,2}\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\theta_{1}^{L},\theta_{2}^{L},k\right)p(\boldsymbol{\theta})p(\theta_{1}^{L})p(\theta_{2}^{L})p(k)d\boldsymbol{\theta}d\theta_{1}^{L}d\theta_{2}^{L}\\\
&=\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\theta_{1}^{L},\theta_{2}^{L},k=1\right)p(\boldsymbol{\theta})p(\theta_{1}^{L})p(\theta_{2}^{L})p(k=1)d\boldsymbol{\theta}d\theta_{1}^{L}d\theta_{2}^{L}+\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\theta_{1}^{L},\theta_{2}^{L},k=2\right)p(\boldsymbol{\theta})p(\theta_{1}^{L})p(\theta_{2}^{L})p(k=2)d\boldsymbol{\theta}d\theta_{1}^{L}d\theta_{2}^{L}\\\
&=\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\theta_{1}^{L},k=1\right)p(\boldsymbol{\theta})p(\theta_{1}^{L})p(k=1)d\boldsymbol{\theta}d\theta_{1}^{L}\textcolor{blue}{\int
p(\theta_{2}^{L})d\theta_{2}^{L}}+\int
p\left(\boldsymbol{d}\mid\boldsymbol{\theta},\theta_{1}^{L},\theta_{2}^{L},k=2\right)p(\boldsymbol{\theta})p(\theta_{1}^{L})p(\theta_{2}^{L})p(k=2)d\boldsymbol{\theta}d\theta_{1}^{L}d\theta_{2}^{L}\end{split}$
(19)

where in the second line we have written out explicitly that each of the terms
corresponding to a different value of $k$, two terms in this case for $k=1$
and $k=2$. Furthermore, the integral corresponding to $k=1$ term is split into
two parameter spaces: $\theta_{1}$ and $\theta_{2}$, where the integral over
$\theta_{2}$ indicated in blue gives 1 by definition. Hence, the conditional
marginal likelihood of $k=1$ is not affected by the lensing parameters for
$k>1$.


g}{\partial\dot{y}_{a}}(\dot{y}_{1},...,\dot{y}_{m},0,...,0).\delta_{bc}$
(5.3) $\displaystyle=$ $\displaystyle 0$

We define the $\mathbb{C}^{\ast}$-equivariant space $\check{V}$ and the
$\mathbb{C}^{\ast}$-equivariant étale morphisms
$\check{\jmath},\check{\alpha},\check{\beta}$ by the
$\mathbb{C}^{\ast}$-equivariant Cartesian square:

${\check{V}}$${\tilde{V}}$${\tilde{U}\times\mathbb{C}^{n}}$${\mathbb{C}^{m+n}}$$\scriptstyle{\check{\jmath}}$$\scriptstyle{\check{\alpha}\times\check{\beta}}$
$\scriptstyle(\tilde{y}_{1},...,\tilde{y}_{m},\tilde{z}_{1},...,\tilde{z}_{n})$
$\scriptstyle(\tilde{x}_{1},...,\tilde{x}_{m})\times{\rm Id}_{\mathbb{C}^{n}}$

where $\mathbb{C}^{m+n}$ is provided with the $\mathbb{C}^{\ast}$-action
induced from the action of $\mathbb{C}^{\ast}$ on the étale coordinates. There
is a unique $\check{v}\in\check{V}$ with $\check{\alpha}(\check{v})=i(x)$,
$\check{\beta}(\check{v})=(0,...,0)$ and $\check{\jmath}(\check{v})=j(x)$.).
We can regard $\check{\jmath}:\check{V}\to\tilde{V}\subset V$ as an étale
$\mathbb{C}^{\ast}$-invariant open set in V. Define
$\mathbb{C}^{\ast}$-equivariant étale coordinates
$(\check{y}_{1},...,\check{y}_{m},\check{z}_{1},...,\check{z}_{n}):\check{V}\to\mathbb{C}^{m+n}$
by
$\check{y}_{a}=\tilde{y}_{a}\circ\check{\jmath},\check{z}_{b}=\tilde{z}_{b}\circ\check{\jmath}$,
and define the $\mathbb{C}^{\ast}$-invariant functions
$\check{f},\check{g}:\check{V}\to\mathbb{C}$ by
$\check{f}=f\circ\check{\alpha},\check{g}=g\circ\check{\jmath}$, and
$\check{h}=\check{g}-\check{f}:\check{V}\to\mathbb{C}$. The previous argument
now shows that on the smooth subscheme $\check{U}\subset\check{V}$ defined by
$\check{z}_{1}=...=\check{z}_{n}=0$ we have $\check{h}|_{\check{U}}$ and
$\frac{\partial\check{h}}{\partial\check{z}_{b}}|_{\check{U}}=0$ for
$b=1,...,n$. Therefore the $\mathbb{C}^{\ast}$ invariant function $\check{h}$
lies in the ideal $(\check{z}_{1},...,\check{z}_{n})^{2}$ generated by
$\mathbb{C}^{\ast}$-equivariant functions on $\check{V}$. So making
$\check{U},\check{V},\tilde{U},\tilde{V}$ smaller but still
$\mathbb{C}^{\ast}$-invariant, we may write
$\check{h}=\sum_{b,c=1}^{n}\check{z}_{b}\check{z}_{c}Q_{bc}$ for some
$\mathbb{C}^{\ast}$-invariant $Q_{bc}:\check{V}\to\mathbb{C}$ with
$Q_{bc}=Q_{cb}$. Up to a $\mathbb{C}^{\ast}$-equivariant linear change of
coordinates $\check{z}$, we can write:

(5.4)
$\displaystyle\sum_{b,c=1}^{n}\check{z}_{b}\check{z}_{c}Q_{bc}(0,...,0)=2\sum_{b=1}^{r}\check{z}_{2b-1}\check{z}_{2b}+\sum_{b=2r+1}^{n}\check{z}_{b}^{2}$

with $n_{2i-1}=-n_{2i}>0$ for $i\leq r$ and $n_{i}=0$ for $i\geq 2r+1$.
Suppose that $n\geq 2r+1$: the $\mathbb{C}^{\ast}$-invariant function
$Q_{nn}(\check{y}_{1},...\check{y}_{m},\check{z}_{1},...,\check{z}_{n})$ is
invertible near $\check{v}$, we can take a double étale cover of
$\check{V}-Q_{nn}^{-1}(0)$ such that $Q_{nn}$ has a square root
$Q_{nn}^{1/2}$. Because $Q_{nn}$ is $\mathbb{C}^{\ast}$-invariant,
$\check{V}-Q_{nn}^{-1}(0)$ is stable under $\mathbb{C}^{\ast}$, hence we can
extend the $\mathbb{C}^{\ast}$ action to the double cover, such that the cover
is $\mathbb{C}^{\ast}$-equivariant. We replace then $\check{V}$ by this étale
neighboorhood of $\check{v}$, and all the coordinates and $Q_{bc}$, and
$Q_{nn}^{1/2}$, are still $\mathbb{C}^{\ast}$-equivariant. We can then write:

$\displaystyle h$
$\displaystyle=\sum_{b,c=1}^{n}\check{z}_{b}\check{z}_{c}Q_{bc}$
$\displaystyle=\sum_{b,c=1}^{n-1}\check{z}_{b}\check{z}_{c}(Q_{bc}-Q_{nn}^{-1}Q_{bn}Q_{cn})+(Q_{nn}^{1/2}\check{z}_{n}+\sum_{b=1}^{n-1}Q_{nn}^{-1/2}Q_{bn}\check{z}_{b})^{2}$
(5.5)
$\displaystyle=\sum_{b,c=1}^{n-1}\check{z}_{b}\check{z}_{c}\hat{Q}_{bc}+z_{n}^{2}$

where $\hat{Q}_{bc}=Q_{bc}-Q_{nn}^{-1}Q_{bn}Q_{cn}$ is $\mathbb{C}^{\ast}$
equivariant of character $-n_{b}-n_{c}$ and
$z_{n}=Q_{nn}^{1/2}\check{z}_{n}+\sum_{b=1}^{n-1}Q_{nn}^{-1/2}Q_{bn}\check{z}_{b}$
is $\mathbb{C}^{\ast}$-invariant, with $\frac{\partial
z_{n}}{\partial\check{z}_{b}}(\check{v})=\delta_{nb}$ and
$\hat{Q}_{bc}(0,...,0)=\delta_{bc}$. By recursion we can find an étale open
neighborhood of $\check{v}$ in $V$, $\mathbb{C}^{\ast}$-invariant functions on
it $z_{b}$ for $b\geq 2r+1$ with $\frac{\partial
z_{b}}{\partial\check{z}_{c}}(\check{v})=\delta_{bc}$, and
$\mathbb{C}^{\ast}$-equivariant functions on it $\check{Q}_{bc}$ for $b,c\geq
2r$ with $\check{Q}_{bc}(0,...,0)=Q_{bc}(0,...,0)$, such that:

(5.6) $\displaystyle
h=\sum_{b,c=1}^{2r}\check{z}_{b}\check{z}_{c}\check{Q}_{bc}+\sum_{b=2r+1}^{n}z_{b}^{2}$

On a $\mathbb{C}^{\ast}$-invariant open neighborhood of $\check{v}$ such that
$\check{Q}_{2r-1,2r}$ is invertible:

$\displaystyle h=$
$\displaystyle\sum_{b,c=1}^{2r-2}\check{z}_{b}\check{z}_{c}(\check{Q}_{bc}-\sum_{b,c=1}^{2r-2}\check{Q}_{2r-1,2r}^{-1}\check{Q}_{b,2r}\check{Q}_{2r-1,c})+\sum_{b=2r+1}^{n}z_{b}^{2}$
$\displaystyle+2(\check{z}_{2r-1}+\sum_{b=1}^{2r-2}\check{Q}_{2r-1,2r}^{-1}\check{Q}_{b,2r}\check{z}_{b})(\check{Q}_{2r-1,2r}\check{z}_{2r}+\sum_{c=1}^{2r-2}\check{Q}_{2r-1,c}\check{z}_{c})$
(5.7) $\displaystyle=$
$\displaystyle\sum_{b,c=1}^{2r-2}\check{z}_{b}\check{z}_{c}\tilde{Q}_{bc}+2z_{2r-1}z_{2r}+\sum_{b=2r+1}^{n}z_{b}^{2}$

where
$\tilde{Q}_{bc}=\check{Q}_{bc}-\sum_{b,c=1}^{2r-2}\check{Q}_{2r-1,2r}^{-1}\check{Q}_{b,2r}\check{Q}_{2r-1,c}$
is $\mathbb{C}^{\ast}$-equivariant of character $-n_{b}-n_{c}$ with
$\tilde{Q}_{bc}(0,...,0)=Q_{bc}(0,...,0)$,
$z_{2r-1}=\check{z}_{2r-1}+\sum_{b=1}^{2r-2}\check{Q}_{2r-1,2r}^{-1}\check{Q}_{b,2r}\check{z}_{b}$
is $\mathbb{C}^{\ast}$ equivariant of character $n_{2r-1}$ with
$\frac{\partial z_{2r-1}}{\partial\check{z}_{b}}(\check{v})=\delta_{2r-1,b}$,
and
$z_{2r}=\check{Q}_{2r-1,2r}\check{z}_{2r}+\sum_{c=1}^{2r-2}\check{Q}_{2r-1,c}\check{z}_{c}$
is $\mathbb{C}^{\ast}$-equivariant of weight $n_{2r}$ with $\frac{\partial
z_{2r}}{\partial\check{z}_{b}}(\check{v})=\delta_{2r,b}$. By recursion, there
is then on an a $\mathbb{C}^{\ast}$-invariant étale open neighborhood
$\jmath^{\prime}:V^{\prime}\to\check{V}$ of $\check{v}$ and
$\mathbb{C}^{\ast}$-equivariant functions $z_{b}$ with $\frac{\partial
z_{b}}{\partial\check{z}_{c}}(\check{v})=\delta_{bc}$ (i.e. by restricting
$V^{\prime}$ to a $\mathbb{C}^{\ast}$-invariant open neighborhood of
$\check{v}$,
$(y_{1}:=y_{1}\circ\jmath^{\prime},...,y_{m}:=y_{m}\circ\jmath^{\prime},z_{1},...,z_{n}$
forms a system of étale coordinates of $V^{\prime}$), such that:

(5.8) $\displaystyle
h=2\sum_{b=1}^{r}z_{2b-1}z_{2b}+\sum_{b=2r+1}^{n}z_{b}^{2}$

Define the $\mathbb{C}^{\ast}$-invariant subset $U^{\prime}=\\{v^{\prime}\in
V^{\prime}:z_{1}(v^{\prime})=...=z_{n}(v^{\prime})=0$. Define the
$\mathbb{C}^{\ast}$-equivariant maps $\iota:U^{\prime}\to
U,\jmath:V^{\prime}\to V,\Phi^{\prime}:U^{\prime}\to
V^{\prime},\alpha:V^{\prime}\to U,\beta:V^{\prime}\to\mathbb{C}^{n}$ by
$\iota:=\check{\alpha}\circ\jmath^{\prime}|_{U}^{\prime}$,
$\jmath:=\check{\jmath}\circ\jmath^{\prime}$, $\Phi^{\prime}={\rm
Id}_{U^{\prime}}$, $\alpha=\check{\alpha}\circ\jmath^{\prime}$ and
$\beta=(z_{1},...,n_{n})$. As $\jmath^{\prime}:V^{\prime}\to\check{V}$ is an
étale open neighborhood of $\check{v}$ in $\check{V}$ , there exists
$u^{\prime}\in V^{\prime}$ with $\jmath^{\prime}(u^{\prime})=\check{v}$, and
$z_{b}(u^{\prime}=0$ for $b=1,...,n$ as $\check{z}_{b}(\check{v})=0$, so
$u^{\prime}\in U$ with
$\iota(u^{\prime})=\check{\alpha}\circ\jmath^{\prime}(u^{\prime})=\check{\alpha}(\check{v})=i(x)$.
Also $\iota,\jmath,\alpha\times\beta$ are étale as
$\check{alpha}\times\check{\beta},\check{\jmath},\jmath^{\prime}$ are. The
same computations as in the proof of [21, Prop 2.23] show that in a
neighborhood of $u^{\prime}$, which can be taken to be
$\mathbb{C}^{\ast}$-invariant, we have $\Phi\circ\iota=\jmath\circ\Phi$, so
this making $u^{\prime}$ smaller and still $\mathbb{C}^{\ast}$-invariant it
holds on $U^{\prime}$. The equations
$\alpha\circ\Phi^{\prime}=\iota,\beta\circ\Phi^{\prime}=0$ are immediate, and
$g\circ\jmath=f\circ\alpha+(2\sum_{b=1}^{r}z_{2b-1}z_{2b}+\sum_{b=2r+1}^{n}z_{b}^{2})\circ\beta$
follows from (5.8). $\Box$

###  5.2. Smooth torus-equivariant morphism of d-critical scheme

######  Lemma 5.2.

Consider a smooth morphism of d-critical scheme $\phi:(X,s)\to(Y,t)$ of
relative dimension $d$. For $x\in X$, we can then choose two
$\mathbb{C}^{\ast}$-equivariant critical charts $(R,U,f,i)$ and $(S,V,g,j)$ of
$X$ and $Y$ containing $x$ and $y:=f(x)$ with a smooth map $\Phi:U\to V$ of
relative dimension $d$ such that $f=g\circ\Phi$ and $\Phi\circ i=j\circ\phi$.

Proof: This is a $\mathbb{C}^{\ast}$-equivariant version of the result used in
the proof of [7, Prop 4.5]: as before, we check that we can do each step in a
$\mathbb{C}^{\ast}$-equivariant way. As in the proof of [21, Prop 2.43],
because the $\mathbb{C}^{\ast}$ action on $Y$ is étale-locally linearizable,
we can find an étale neighborhood $S^{\prime}$ of $y$ with a closed
$\mathbb{C}^{\ast}$-equivariant embedding into a smooth algebraic space with
$\mathbb{C}^{\ast}$-action $V^{\prime}$, with $\dim(V^{\prime})=\dim(T_{y}Y)$.
According to the proof of the proposition 2.4, the
$\mathbb{C}^{\ast}$-equivariant smooth map $\phi^{-1}(S^{\prime})\to
S^{\prime}$ is étale locally a $\mathbb{C}^{\ast}$-equivariant affine
fibration, hence we can find a $\mathbb{C}^{\ast}$-equivariant étale
neighborhood $S$ of $y$, $\mathbb{C}^{\ast}$-equivariant embeddings
$i:R:=\phi^{-1}(S)\to U$, $j:S\to V$ into smooth algebraic spaces with
$\mathbb{C}^{\ast}$-action $U,V$ with dimension $\dim(U)=\dim(T_{x}X)$,
$\dim(V)=\dim(T_{y}Y)$, with a $\mathbb{C}^{\ast}$-equivariant smooth map of
relative dimension $d$ $\Phi:U\to V$ such that $\Phi\circ i=j\circ\phi|_{R}$.
Hence as in the proof of [21, Prop 2.43], up to shrinking $S,V$ (and then
$R,U$) to a $\mathbb{C}^{\ast}$-invariant open subspace we can find a
$\mathbb{C}^{\ast}$ invariant regular function $g:V\to\mathbb{C}$ such that
$\iota_{S,V}(t|_{S})=j^{-1}(g)+I_{S,V}^{2}$, and $S={\rm Crit}(g)$, i.e.
$(S,V,g,j)$ is a $\mathbb{C}^{\ast}$-equivariant critical chart of $Y$ near
$y$. Defining the $\mathbb{C}^{\ast}$-equivariant regular function
$f:=g\circ\Phi:U\to\mathbb{C}$, because $s=\phi^{\star}(t)$, one has
$\iota_{R,U}(s|_{R})=i^{-1}(f)+I_{R,U}^{2}$. Also $\Phi\circ
i=j\circ\phi|_{R}$ and $\Phi,\phi$ are smooth of relative dimension $d$, then
$i(R)=\Phi^{-1}(j(S))$ in a neighborhood of $i(x)$: taking the union of all
these neighborhood, one can shrink $U$ (and then $R$) to
$\mathbb{C}^{\ast}$-invariant open neighborhood satisfying
$i(R)=\Phi^{-1}(j(S))$. Because $\Phi$ is smooth, one has:

(5.9) $\displaystyle{\rm Crit}(f)=\Phi^{-1}({\rm
Crit}(g))=\Phi^{-1}(j(S))=i(R)$

i.e. $(R,U,f,i)$ is a $\mathbb{C}^{\ast}$-equivariant critical chart near $x$.
We have then build $\mathbb{C}^{\ast}$-equivariant critical charts $(R,U,f,i)$
and $(S,V,g,j)$ of $X$ and $Y$ near $x$ and $y$ with a smooth map $\Phi:U\to
V$ of relative dimension $d$ such that $f=g\circ\Phi$ and $\Phi\circ
i=j\circ\phi$. $\Box$

## References

  * [1] Jarod Alper, Jack Hall, and David Rydh, _The étale local structure of algebraic stacks_ , arXiv: Algebraic Geometry (2019), arXiv:1912.06162v2.
  * [2] by same author, _A Luna étale slice theorem for algebraic stacks_ , Annals of mathematics 191 (2020), arXiv:1504.06467v3.
  * [3] Jarod Alper, Daniel Halpern-Leistner, and Jochen Heinloth, _Existence of moduli spaces for algebraic stacks_ , Inventiones mathematicae (2018), arXiv:1812.01128.
  * [4] Noah Arbesfeld, _K-theoretic Donaldson-Thomas theory and the Hilbert scheme of points on a surface_ , Algebraic Geometry (2021), 587–625, arXiv:1905.04567v2.
  * [5] K. Behrend and B. Fantechi, _Symmetric obstruction theories and Hilbert schemes of points on threefolds_ , 2008, arXiv:0512556v1.
  * [6] Kai Behrend, _Donaldson-Thomas type invariants via microlocal geometry_ , Ann. of Math. (2) 170 (2009), no. 3, 1307–1338, arXiv:math/0507523.
  * [7] O. Ben-Bassat, Christopher Brav, Vittoria Bussi, and D. Joyce, _A ‘Darboux theorem’ for shifted symplectic structures on derived Artin stacks, with applications_ , Geometry & Topology 19 (2015), arXiv:1312.0090v2.
  * [8] Andrzej Białynicki-Birula, _Some theorems on actions of algebraic groups_ , Annals of Mathematics 98 (1973), 480–497.
  * [9] Tom Braden, _Hyperbolic localization of intersection cohomology_ , Transformation Groups 8 (2002), 209–216, arXiv:math/0202251v5.
  * [10] Christopher Brav, Vittoria Bussi, Delphine Dupont, Dominic Joyce, and Balázs Szendrői, _Symmetries and stabilization for sheaves of vanishing cycles_ , Journal of singularities 11 (2015), arXiv:1211.3259v4.
  * [11] Christopher Brav, Vittoria Bussi, and D. Joyce, _A ’Darboux theorem’ for derived schemes with shifted symplectic structure_ , J. Amer. Math. Soc. 32 (2019), arXiv:1305.6302v4.
  * [12] Christopher Brav and Tobias Dyckerhoff, _Relative calabi–yau structures ii: shifted lagrangians in the moduli of objects_ , Selecta Mathematica 27 (2018), arXiv:1812.11913.
  * [13] Vittoria Bussi, Dominic Joyce, and Sven Meinhardt, _On motivic vanishing cycles of critical loci_ , Journal of Algebraic Geometry (2019), arXiv:1305.6428v3.
  * [14] B. Davison and Sven Meinhardt, _Cohomological Donaldson–Thomas theory of a quiver with potential and quantum enveloping algebras_ , Inventiones mathematicae (2016), arXiv:1311.7172.
  * [15] Pierre Descombes, _Cohomological dt invariants from localization_ , Journal of the London Mathematical Society 106 (2021), arXiv:2106.02518.
  * [16] Vladimir Drinfeld, _On algebraic spaces with an action of $G_{m}$_, arXiv preprint arXiv:1308.2604 (2013), arXiv:1308.2604.
  * [17] Tom Graber and Rahul Pandharipande, _Localization of virtual classes_ , Inventiones mathematicae 135 (1997), 487–518, arXiv:alg-geom/9708001v2.
  * [18] Amer Iqbal, Can Kozcaz, and Cumrun Vafa, _The refined topological vertex_ , JHEP 10 (2009), 069, arXiv:hep-th/0701156.
  * [19] Joachim Jelisiejew and Łukasz Sienkiewicz, _Białynicki-Birula decomposition for reductive groups_ , Journal de Mathématiques Pures et Appliquées 131 (2019), 290–325, arXiv:1805.11558v3.
  * [20] Yunfeng Jiang, _The moduli space of stable coherent sheaves via non-archimedean geometry_ , arXiv preprint (2017), arXiv:1703.00497.
  * [21] D. Joyce, _A classical model for derived critical loci_ , arXiv: Algebraic Geometry (2013), arXiv:1304.4508.
  * [22] M. Kontsevich and Y. Soibelman, _Cohomological Hall algebra, exponential Hodge structures and motivic Donaldson-Thomas invariants_ , arXiv: Algebraic Geometry (2010), arXiv:1006.2706v2.
  * [23] Maxim Kontsevich and Yan Soibelman, _Stability structures, motivic Donaldson-Thomas invariants and cluster transformations_ , arXiv:0811.2435.
  * [24] David B. Massey, _The Sebastiani–Thom Isomorphism in the Derived Category_ , Compositio Mathematica 125 (2001), arXiv:math/9908101v1.
  * [25] David B. Massey, _Natural commuting of vanishing cycles and the Verdier dual_ , Pacific Journal of Mathematics 284 (2016), 431–437, arXiv:0908.2799v1.
  * [26] Sergey Mozgovoy and Markus Reineke, _On the noncommutative Donaldson-Thomas invariants arising from brane tilings_ , Adv. Math. 223 (2008), no. 5, 1521–1544, arXiv:0809.0117.
  * [27] Hiraku Nakajima, _Lectures on perverse sheaves on instanton moduli spaces_ , arXiv: Representation Theory (2016), arXiv:1604.06316v3.
  * [28] Nikita Nekrasov and Andrei Okounkov, _Membranes and sheaves_ , Algebraic Geometry 3 (2016), arXiv:1404.2323.
  * [29] T. Pantev, B. Toën, M. Vaquié, and G. Vezzosi, _Shifted symplectic structures_ , Publ.math.IHES 117 (2013), arXiv:1111.3209v4.
  * [30] Miroslav Rapcak, Yan Soibelman, Yaping Yang, and Gufang Zhao, _Cohomological Hall algebras, vertex algebras and instantons_ , Commun. Math. Phys. 376 (2019), arXiv:1810.10402.
  * [31] Timo Richarz, _Spaces with $G_{m}$-action, hyperbolic localization and nearby cycles_, arXiv: Algebraic Geometry (2016), arXiv:1611.01669v3.
  * [32] Morihiko Saito, _Mixed hodge modules_ , 1990.
  * [33] Hideyasu Sumihiro et al., _Equivariant completion_ , Journal of Mathematics of Kyoto University 14 (1974), no. 1, 1–28.
  * [34] Balazs Szendroi, _Cohomological Donaldson-Thomas theory_ , (2015), arXiv:1503.07349.
  * [35] Richard P Thomas, _A holomorphic Casson invariant for Calabi-Yau 3-folds, and bundles on K3 fibrations_ , arXiv preprint math/9806111 (1998).
  * [36] Yukinobu Toda, _Gopakumar–vafa invariants and wall-crossing_ , Journal of Differential Geometry (2017), arXiv:1710.01843.