Gravitational wave lensing as a probe of halo properties and dark matter

,

Gravitational lensing is a consolidated astronomical probe across vastly different scales.From weak lensing of the cosmic microwave background to microlensing by exoplanets orbiting Milky Way stars, its rich phenomenology offers valuable insights about the matter distribution in the Universe [1].Often, gravitational lensing data translates into powerful tests of fundamental physics, such as the properties of dark matter (DM).Given the breadth of DM theories still viable [2], increasing the pool of available observations is of paramount importance.
Applications of gravitational lensing have relied mainly on sources observed via electromagnetic (EM) radiation.Rapid progress of gravitational wave (GW) astronomy may soon open a new window into the Universe with the advent of GW lensing. 1 Lensing of GWs is highly complementary to observations of lensed EM sources due to several key properties, such as low emission frequency, phase coherence, negligible absorption and accurate source models.
The lower frequencies of GWs enable the observation of wave-optics (WO) lensing phenomena caused by diffraction.The shortest observed GWs have wavelengths (∼ 1000 km) orders of magnitude larger than the longest EM waves able to penetrate the atmosphere (∼ 10 m).The rich WO phenomenology includes frequency-dependent patterns in the GW signal, which carry information about the lens' mass and its density distribution.WO amplification contains a wealth of information that is lost in Geometric Optics (GO), the high-frequency limit, characterised by only three parameters per image (magnification, time delay and Morse phase), which depend on the local properties of the lensing potential at the location of the images.In contrast, WO effects probe a large portion of the lens, potentially allowing for a more detailed lens reconstruction [3][4][5][6].
The phase coherence of GWs allows one to measure the interference between different GO images in the strong-lensing regime and more easily observe fainter images.For most EM sources the interference patterns are washed out due to temporal or spatial incoherence [7,8].Observing these effects on gravitationally-lensed EM sources would require a very abundant population of light, compact lenses and can be limited by the source's physical size [9,10].Coherence makes GW detectors sensitive to the field's amplitude (rather than its intensity) which decays as 1/D S , rather than 1/D2 S (here D S is the luminosity distance from the source): this scaling will allow the next generation of ground detectors to observe every stellar-mass binary black hole merger in the Universe [11].Detecting the field amplitude has the additional advantage of making fainter images comparably easier to detect, as they are dimmed by a factor √ µ I (instead of µ I , larger when µ I < 1).Here µ I is the lensing magnification factor.The negligible absorption of GWs allows the observation of images through dense or opaque regions, which would either block or outshine EM signals.Strong lensing by smooth matter distributions predicts the formation of faint images near the center of the lens [12].Combined with the comparably smaller demagnification suffered by GWs, the observation of these images may allow GWs to probe the centers of galactic and dark matter halos.We therefore use "GW lensing" instead of "Gravitational lensing of GWs", as repeating "gravitational" is both redundant and wordy.
This type of detection is challenging for EM sources: only two doubly-lensed quasars have been observed with central images where the lens is a single galaxy, with magnification ratios µ H /µ brightest = 0.004, 0.007 (see [13], Sec.4.4), where µ H is the magnification of the central image.The central regions of lenses are where differences between dark matter models become more apparent and may even allow the observation of super-massive black holes [14][15][16][17][18][19][20][21][22].A central image in a cluster-scale lens has been used to constrain the mass of a central black hole [23].GWs can also be modelled accurately.Analytic and numerical methods enable accurate waveform predictions in terms of the source parameters (e.g.masses and spins) [24][25][26].In contrast, EM sources can rarely be described from first principles.Lensing of EM radiation thus requires observing the time variation of these sources [27,28], or using objects that can be calibrated through empirical relations, such as Type Ia supernovae [29].The existence of accurate, well-understood models provides an additional handle to test GW lensing effects in general systems, without additional assumptions.In addition, GWs allow an exquisite timing, making time-delay measurements more precise and robust than for EM systems, such as quasars [30] and supernovae [7,31].
Thanks to these properties, GW lensing offers a synergy to probe the matter distribution in the Universe.Detection of lensed GWs by LISA could be used to accurately reconstruct the lens parameters [4].While strongly-lensed LISA sources with WO effects are unlikely [32], WO effects can be detected for sizeable impact parameters, substantially boosting the chance of detection [4,[33][34][35].These determinations benefit from accurate waveforms, with additional information from higher harmonics improving the lens reconstruction [36].These studies focused on symmetric and singular lenses: the point lens and/or singular isothermal sphere (SIS). 2  Work on a simple extension of the SIS [51] and symmetric but non-singular lenses (Navarro-Frenk-White, NFW) [52,53] suggests that the mismatch between waveforms may allow future observations to distinguish among different lens models.However, a more detailed analysis of the lens reconstruction, e.g.including degeneracies between parameters, needs to be performed.
GW lensing may offer constraints on DM models complementary to other gravitational probes [2,54].A prime target has been compact objects such as primordial black holes and compact and/or light DM halos [40,44,53,[55][56][57][58][59][60].GW may be able to probe the properties of DM further.Long-range interactions of DM particles may allow them to form compact substructures [61][62][63][64].Other DM scenarios have distinct predictions on small scales, such as a lower abundance of sub-structure or central cores in DM halos.Such is the case of ultra-light boson fields (also known as ULDM) [65][66][67][68] or self-interacting DM (SIDM) [69,70], which also address discrepancies between ΛCDM and observations on small scales [71].Probing these features with GWs is complementary to lensed EM sources.
No clear detection of lensed GWs has been reported yet [80,86].However, Ref. [87] found an intriguing stronglensing candidate.Ultimately, ruling out the presence of lensed events requires assumptions about the highredshift merger rate [88].Regardless of its current status, lensing of GWs is bound to become a reality as the number of signals grows with cumulative observing time and detector upgrades [37,[89][90][91][92][93][94].Estimates indicate event rates ∼ 1/yr for LIGO A+ and ∼ 50/yr for 3G detectors [94], with a caveat of false alarm events [95] and large uncertainties on the merger rates at high redshift.In addition to its potential to probe cosmic structures, GW lensing may bias the inferred source population [96] as well as distance measurements [97,98].
The purpose of this work is to address the capacity of strongly-lensed GWs to constrain the properties of galactic and DM halos by exploiting their properties complementary to EM observations.We review gravitational lensing in the WO regime in Sec.II.In Sec.III we study the phenomenology of two symmetric lenses that generalise the SIS by varying its slope and introducing an inner core.Section IV presents a forecast on the detectability of lens features by LISA and LIGO using a Fisherinformation matrix approach.As a potential application, Sec.V explores the prospect of recovering the lens' redshift and its virial mass, as well as probing the parameter space of two dark matter scenarios (self-interacting and ultra-light).We summarise and discuss our results in Sec.VI.Appendices contain further details about lensing, the Fisher forecast calculation and its validation.

A. Summary and guide for busy readers
Our main results encompass the following applications of GW lensing: • Detecting central images: GWs can probe faint images forming near the center of strong gravitational lenses, for which WO corrections are the largest.This is explored in detail in Sec.III for lenses with variable slope and an inner core.
• Probing lens features: In Sec.IV we present a forecast of the sensitivity of lensed GWs to the lens' parameters, including the core size and density slope.Precise measurements are possible thanks to multiple images and WO effects.
• Testing large-scale structure: In Sec.V A we develop a probabilistic framework to constrain the lens' redshift and virial mass, given observations of projected quantities, an expansion history and halo mass function.
• Constraining dark matter: we show the capacity of lensed GWs to set stringent limits on dark matter theories that predict the formation of cores.Self-interacting dark matter and ultra-light dark matter are discussed in Secs.V B and V C.
Our analysis highlights the differences and synergies between lensing of EM and GW sources.
Our conclusions (Sec.VI) summarize our analyses, highlight our main results and discuss implications.Readers are encouraged to skip content they are familiar with.We have provided abundant cross-references and figures describing our results and analysis.Readers interested in a specific lens model can read only the relevant subsections.

II. GRAVITATIONAL WAVE LENSING
In this Section, we will review the equations governing gravitational lensing in the wave-optics (WO) regime, Sec.II A. We will recap analytic expansions valid in the high-frequency limit in Sec.II B. The low-frequency expansion and numerical methods are discussed in Appendices A 1 and A 2, respectively.

A. Equations & definitions
The WO regime of gravitational lensing is characterised by the amplification factor F (f ), defined as the ratio between the lensed and unlensed waveforms in the frequency domain: where f is the frequency. 3At the leading order in frequency (f much higher than the typical background curvature), the polarization tensor of the wave is just parallelly transported along the null geodesics of the wave [99].Therefore, the effect of lensing on the polarization is negligible, and the polarization tensor can be regarded as a constant, so it drops out in the definition of F (f ).
Polarization-dependent corrections on the GW phase are suppressed by the curvature R and frequency by a factor ∼ R/f 2 [100,101].Here and in the following we work in units where c = 1.Let us consider a lens at redshift z L and define an effective distance Here D L , D S and D LS are the angular diameter distances to the lens, the source and between the lens and source.We will work within the thin-lens approximation, projecting the mass of the lens onto the lens plane [102].Moreover, we indicate positions on the lens plane with ξ and the impact parameter on the source plane with η, where both are two-dimensional vectors.The observerlens-source configuration is summarised in Fig. 1.
The expression for F (f ) can be obtained from the simplified Fesnel-Kirchhoff integral [12] (a static configuration has been assumed in deriving this result).Here t d (ξ, η) is the time-delay function of the lens, given by The overall phase φm (η) is chosen such that the minimum time delay is zero, meaning that the first component of the lensed signal to be received (a type I image in the GO limit, cf.Sec.II B) arrives at t d = 0.The lensing potential ψ(ξ) is determined by the projected mass distribution of the lens.In particular, given a density ρ(r), the projected mass density Σ(ξ) is obtained by integrating along the direction z perpendicular to the lens plane and ψ(ξ) is the solution of the equation where ∇ 2 ξ is the 2-dimensional Laplacian and G is Newton's constant.It is convenient to recast the diffraction integral of Eq. ( 3) in terms of dimensionless quantities.To this end, we introduce two scales, ξ 0 and η 0 ≡ D S ξ 0 /D L , that will be specified depending on the lens model.Then, we define the dimensionless quantities This allows one to construct dimensionless versions of the time delay ( 4) and of the lensing potential (6) as follows The rescaled lensing potential is given by the dimensionless version of Eq. ( 6), ∇ 2 x ψ(x) = 2κ(x), in terms of the convergence where the critical density is Σ cr ≡ (4πG The potential ψ(x) can be then obtained using the Green's function method Here and in the following, log z indicates the natural logarithm.In Eq. (11), one needs to impose the proper boundary conditions by adding solutions of the homogeneous Laplace equation to the right-hand side.In the context of the geometric-optics approximation discussed in the following subsection, it is also useful to define the reduced deflection angle The dimensionless version of the time delay, also known as Fermat potential, ϕ(x, y), takes the simple form From here on, we will suppress in our formulas the minimum of the Fermat potential ϕ m (y) and always assume that it is added to make the minimum arrival time equal to zero.When necessary, we will introduce it back.All these definitions, when applied to the diffraction integral (3), lead to the following expression Here we introduced the dimensionless frequency where the redshifted effective lens mass is given by In the point-lens case, M Lz corresponds to the total (redshifted) mass (1 + z L )M , i.e. setting the scale in Eq. ( 7) to the Einstein radius However, for extended lenses, M Lz does not correspond to the total lens mass and may differ by several orders of magnitude from the virial mass, cf.Eq. (A10). 4We will discuss the relationship between M Lz and the total halo mass for extended lenses in Sec.III.The role of the (unknown) lens redshift will be examined in Sec.V (cf.Fig. 21).
Once computed, F (w) (for a given lens and impact parameter) can be compared at different source frequencies and lens masses.In typical situations, the WO lensing phenomena can be broadly characterised in three distinct regimes depending on the dimensionless frequency • Perturbative: for w ≪ 1 the signal undergoes a small amplification F (w) ∼ 1 + w a , with a determined by the asymptotic properties of ψ(x) (Sec.A 1).This limit reflects how a wave is not affected by objects smaller than its wavelength.
• "Intermediate" wave-optics: at finite w the amplification receives contributions from large portions of the lens plane (Sec.A 2). Around w ∼ 1 we typically observe the onset of amplification, i.e. the first peak in F (w).
• Geometric optics: for w ≫ 1 the amplification factor is dominated by GO images (stationary points in the lens plane, where ∇ x ϕ(x, y) = 0).Each image gives a copy of the signal characterised by a magnification, time delay and phase factor (Sec. II B).
The three regimes are shown in Fig. 2 for a typical LVK source.Note that the accuracy of the perturbative and GO limits depends on the lens parameters and source's position, a dependence that will be discussed later.Specializing in GW detectors on the ground and in space, the typical order of magnitude for the dimensionless frequency is The onset of magnification corresponds to lenses in the range of massive stellar objects and intermediate-mass black holes (M Lz ∼ 10 − 1000 M ⊙ ) for typical LVK sources, while for LISA it corresponds to sub-halos and massive black holes M Lz ∼ 10 5 − 10 7 M ⊙ .Proposed lower-frequency detectors such as µARES in the µHz range would push the onset of magnification to even higher masses M Lz ∼ 10 10 M ⊙ [103], while for pulsar timing arrays in the nHz band it would reach ∼ 10 13 M ⊙ .

B. Geometric optics & high-w expansion
In the high-frequency limit, only the neighbourhoods of extrema of the Fermat potential (13) contribute to the amplification factor (14).Each extremum is associated with an image J, located at a position x J in the image plane where the lens equation is satisfied.The geometric-optics regime emerges from a quadratic expansion of the Fermat potential around each image, so the diffraction integral can be performed analytically.
The GO amplification factor ( 14) receives contributions from each image J where the magnification is evaluated on the image position x J (the second equality above applies to the specific case of axially-symmetric lenses).In the above expressions, ϕ J is the time delay (in units of 4GM Lz , Eq. ( 13)) of the J-th image and ϕ ,ij ≡ ∂ i ∂ j ϕ is its Hessian matrix.The Morse Phase [4,12] depends on the type of image as Minima, saddle points and maxima of the time delay function are also known as type I, II and III images, respectively.
Beyond geometric optics (bGO) corrections can be obtained as a series expansion in 1/w.We now present the expressions derived in Ref. [5] (see also Ref. [104]), focusing on axially-symmetric lenses.The amplification at order 1/w reads where the real number ∆ J characterizes the bGO correction for each image and is given by J /x J )/2.Equation (24) shows that the leading order GO result is a good approximation provided that ∆ J /w ≪ 1 for all images.
On top of the bGO corrections, originating from the locations of the images, subleading terms in 1/w also appear from non-analytic points in the Fermat potential (e.g.cusps), without the presence of a corresponding image [5,104]. 5In general, these effects need to be accounted for at intermediate frequencies.
Let us summarise the results obtained in [5] regarding contributions from cusps at high w (see also [104]).We can focus on a particular lens model, the gSIS, that we will discuss at length in this work (see Tab. I).Following the same logic as in the stationary-phase approximation (that leads to the GO and bGO amplification factors), one can isolate the contribution of a cusp by expanding ϕ(x, y) around the location of the singular point and then integrating in a small neighbourhood around it (taking a small interval is justified when w ≫ 1).In this case, the integral is not simply approximated by a Gaussian since the location of the cusp is not a stationary point of the Fermat potential.Let us focus on the gSIS lens for k ≥ 1, which has a cusp at the origin in the lens plane.We indicate with F c (w) the contribution to the amplification factor from the cusp.For small impact parameter and high w, F c (w) is found to be where ϕ c ≡ y 2 /2 − ϕ m is the time delay of the center of the lens and Γ(z) is the Gamma function.Notice that in the SIS case, k = 1, the cusp contribution becomes F c (w) ≃ i e iwϕc /w while for k > 1 Eq.( 26) decays faster than 1/w.We will indicate with rGO the approximation for F (w) that includes both bGO terms and F c (w) at high w.

III. MODELS OF LENS FEATURES
We now address the GW lensing phenomenology of the extended lens models.We focus on two 1parameter extensions of the singular isothermal sphere (SIS, briefly reviewed in Appendix A 3).The generalised SIS (gSIS) allows a variable slope of the density distribution (Sec.III A).The cored isothermal sphere (CIS) introduces a central core with finite density (Sec.III B).In both cases, we will present the relation between the halo mass and M Lz , the GO structure of the lens as well as bGO and WO features.In Sec.III C we address parameter reconstruction, detectability of central images and the mismatch with GWs lensed by an SIS.
The lens models are summarised in Table I. Figure 3 shows the matter profile, lensing potential and lens equation for a characteristic example within each lens model, highlighting some of their differences.For the gSIS the narrow and wide cases are qualitatively different and are shown separately.

A. Power law generalised SIS
We now consider a generalization of the SIS lens with a generic power-law profile.This model, which we refer to as the gSIS model, is discussed in Refs.[12,51,105].We will present the lens characteristics, its GO/bGO properties and WO signatures.For an illustrative comparison with the SIS and CIS lenses see Fig. 3.

Mass profile and scales
The density profile of the gSIS model is given by where ρ 0 is a typical value for the density, r ⋆ is a radial scale and k represents the slope of the halo: k = 1 recovers the SIS, and larger/smaller values correspond to steeper/shallower profiles, respectively.Taking the range 0 < k < 2 ensures both finite central densities and that the lensing potential grows more slowly than the quadratic part of the Fermat potential.From Eq. ( 5), the projected mass density of the lens is where for convenience we defined the constant . We make the following choice for the scale ξ 0 : Then, using the relations of Sec.II A we obtain that the lensing potential is For all the allowed values of k, the total mass of the lens is divergent, as in the case of the SIS lens.Therefore, a cutoff in the radius is needed in order to define a virialized mass.In this case, M vir is obtained as where, as for the SIS, ρ vir ≡ 200ρ c .This expression allows us to connect the effective lens mass M Lz with M vir .To do so, we can use Eq. ( 31) to obtain ρ 0 in terms of M vir , and then replace this quantity in the expression for M Lz .We can write this relation as follows where the dimensionless coefficient γ k and the mass scale M 0 are defined as This generalises the SIS result, Eq. (A11).Notice that our choice for the normalisation scale ξ 0 and the effective lens mass M Lz , Eqs. ( 29) and (32), depend on the slope k.When comparing lensing results for different ks, one has to keep in mind that these scales are not kept fixed.

GO structure & bGO corrections
The steepness parameter k leads to two distinct regimes: 1.A broad matter profile (0 < k < 1) leads to the formation of a third GO image.This central image is the closest to the lens center.It is associated with the maximum of the Fermat potential and has finite magnification.
2. A narrow lens profile (1 < k < 2) has only two GO images.However, both exist for arbitrarily large impact parameters, with the type II image becoming very faint for large y.This is similar to a point lens and is due to the compactness of the lens.
Note that a broad matter profile is associated with a steep lens potential, while a narrow matter profile is associated  with a shallower lens potential; compare Eqs.(27) and (30).The SIS is the limiting case between the two.The image positions, magnifications and bGO corrections for the broad/narrow lenses are shown in Figs. 4 and 5, respectively.
Let us see these differences in detail in the geometricoptics limit, first focusing on the critical curves (i.e. points where the magnification Eq. ( 22) diverges, see also Ref. [12]).Inspecting Eq. (22), one sees that one family of solutions, the so-called tangential curves, are given by the solutions to the equation α(x)/x = 1.For the gSIS lens, using Eq. ( 12), we have α(x) = x 1−k .Thus, for general k the tangential curves are given by x tc = 1.This is also the case for the SIS lens, so we have no difference for what concerns the Einstein's ring.On the other hand, the other family of solutions, the radial critical curves, are given by the equation dα(x)/dx = 1.They have a more interesting dependence on k.The solution is indeed and exists only for k < 1 (otherwise Eq. (34) leads to a complex solution).The SIS corresponds to the limiting case k = 1, where this curve approaches the origin and vanishes.In the case k < 1, the expression Eq. ( 34) together with the lens equation ( 20) leads to a caustic in the image plane Equation (35) gives us the region where multiple images form if k < 1.In the SIS case, this reduces to y rc = 1 as expected.For k < 1 instead, the value of y rc is always less than 1.
At this point, we can consider the magnification of the different images.The magnification for a given image x I is From this relation, we see that the critical curves are infinitely magnified as expected.Moreover, images close to the centre are highly demagnified since µ I scales as x 2k I .Central images vanish when x I → 0. However, even in this limit, they can be observed via bGO corrections.
The bGO corrections (parametrised by the coefficient  30) with k = 1/2).Top: Image positions associated with saddle point (blue), maximum (orange) and minimum (green) of the Fermat potential, with the projected density (shaded region).The saddle and maximum images exist only within the caustic (red dashed vertical line, Eq. ( 35)).The shaded area represents Σ as a function of x.Its height is rescaled to arbitrary units.Middle: Magnification for the GW amplitude for the different images from Eq. (36).The central image magnification goes to zero at low y and diverges near the caustic.Bottom: Beyond geometric optic correction times the magnification for different images, obtained from Eq. (37) and Eq.(36).The central image has a large contribution that overcomes demagnification at low y.
∆ given in Eq. ( 25)) for a generic image reads Interestingly, for a central image this correction can be- come very large given that ∆ I ∼ 1/x 2−k I .In order to assess whether bGO corrections from this image overcome the other two (brighter) terms one should look at the combination ∆ I |µ I | 1/2 since, as we previously saw, the magnification tends to zero in the same limit.This combination scales as , hence is large at small enough x I when k < 1 (which is the same situation where there is a central image in the first place).Moreover, the bGO correction is suppressed at high frequencies.Therefore, we expect this central image to be important in the wave-optics regime w ∼ 1, at least in a certain range of values for the impact parameter.
We will now explore these results explicitly for some values of k where the lens equation can be solved analyt- FIG. 6.Amplification factor of a gSIS for different values of the slope k, for a fixed value of the impact parameter y = 0.3 (in the strong-lensing regime).The WO curves are obtained with the numerical method of Sec.A 2 and are shown with solid lines.The GO result is shown in dashed lines (on the left-hand side GO is shown only for one curve).On the right-hand side, the plots are in linear scale.Here, the envelope of the WO (GO) |F (w)| is shown as thick solid (dashed) line.Top: case k < 1 where three images form within the caustic and give additional modulations in w (see right panel).Bottom: case k ≥ 1 where only two images form.In this case, no large modulations appear, and the curves quickly approach the GO approximation.

ically.
a. Broad profiles 0 < k < 1.We consider the specific case k = 1/2, where analytic solutions for the image positions exist.Although significantly broader than the SIS, this case has the same qualitative features found for all the models with 0 < k < 1.The geometric-optics results for k = 1/2 are plotted in Fig. 4. Let us discuss this case in more detail.The lens equation Eq. ( 20) is solved for x 6 Here, with an abuse of notation, we consider x I as the image position along the direction of y.Thus, it differs from the radial position x as it can also take negative signs.
where, x L , x H and x S are respectively the minimum, maximum and saddle points of the Fermat potential.The maximum and the saddle merge at the caustic of Eq. ( 35), i.e. at y rc = 1/4, or x rc = 1/4, and for larger impact parameters only the minimum image remains.The magnification, time delays and beyond geometric optics corrections can all be easily derived from these solutions in closed form.For instance, the magnification for the image x I is obtained by plugging the solutions Eq. ( 40) inside Eq. (36).
As already noted for general k, the central image x H has a vanishing magnification for y = 0.In the case at hand moreover, its magnification remains small relative to the others, unless we are very close to the caustic.However, its beyond geometric optic correction ∆ is much larger than for the other two images.This can be seen by computing ∆ from Eq. (37).Clearly, the central image has an enhanced value of ∆ because of the term in the denominator ∼ x 3/2 I .As in the case of µ, we also have a divergence at the caustic.We will come back to this issue after discussing the narrow lens.
b. Narrow profiles 1 < k < 2. Analytic solutions exist for a narrow lens with k = 3/2.In this case, the lens equation reduces to a cubic equation in √ x with a simple closed-form solution: where z ± ≡ ∓2y 3 + 27(1 ± 1 ∓ 4y 3 /27).The minimum x L corresponds to z + while the saddle x S corresponds to z − .The minimum image starts at x L = 1 for y = 0 and moves to larger values of x as y increases.The saddle starts instead at x S = −1 and moves towards zero.Contrarily to the broad case, the two images exist for all values of y, and there is no caustic.We summarise the GO properties of this case in Fig. 5.The fainter image (saddle point) has a lower magnification than the minimum away from y = 0.This can be seen by inspecting the expression for the magnification (36).That expression indeed shows that the saddle gets very demagnified as we increase the impact parameter.
Even though the saddle-point image becomes undetectable at large y, bGO corrections may offer a window to probe compact lenses at high impact parameters.For narrow gSIS at y ≳ 1, the second image approaches the center of the lens and receives large bGO corrections due to the curvature of the lensing potential (its derivatives become large), and so ∆ S ≫ ∆ L .While for y ≳ 1 the second image is already very faint, its associated bGO correction remains sizeable.Equations (36) and (37) evaluated on x S → 0 yield where the last relationship uses the position of the saddle point |x S | ∼ y −1/(k−1) for large y (this holds only for 1 < k < 2, as x S does not exist for y > y rc otherwise).While Eq. ( 42) tends to zero for large y, it remains sizeable for y ≳ 1.This is seen in the bottom panel of Fig. 5, where |µ S ||∆ S | ∼ O(0.1) up to y = 4, with the 1/y 2 behaviour appearing at larger impact parameter.
As the lensing probability is proportional to y 2 , the bGO feature opens the possibility of detecting steep lenses in high SNR, low frequency (in terms of w) GW events via bGO corrections.

Wave optics features
Here we discuss the phenomenology of the gSIS lens and the effect of the slope k on the full WO amplification.
a. Imprints of the lens slope k: Let us now discuss the effect of the gSIS slope k on the amplification factor.Figure 6 shows the full WO predictions in the stronglensing regime (fixed y = 0.3) for both broad (k < 1, top panel) and SIS/narrow (k ≥ 1, bottom panel) lenses.For easier display, each plot is divided between small w in log scale (left) and high w in linear scale (right), where also the amplitude is shown.Some of the features observed in the Figure depend on the relation between the impact parameter y and the caustic y cr , which is a function of the slope via Eq.(35).For k = 0.65, 0.8, 0.95 and 1 the caustic is located at y cr = 0.369, 0.535, 0.811 and 1, respectively.
At low frequencies and fixed y, the effect of k is to provide a different power-law behaviour in w, with F (w) − 1 ∝ w k/2 , Eq. (A4).One can heuristically interpret these results as a wave not being disturbed by objects smaller than its wavelength: broader lenses (small k) converge more slowly to the free propagation case as w ≪ 1, while the narrower lenses converge faster.In the narrowest possible gSIS (k → 2), the convergence is as fast as a point lens F − 1 ∝ w, cf.Eq. (A2).The slope k also affects the amplitude of the term F (w) − 1, with smaller k (broader lenses) producing larger |F (w) − 1|.The larger amplitude is carried over also to the regime w ≳ 1.
The high-frequency behaviour depends on the case under consideration.SIS and narrow lenses (k > 1) display a beating pattern in |F (w)|, caused by the interference between the two GO images.Narrow lenses k > 1 always form two images, and hence this pattern persists even for y ≫ 1 although with a small amplitude modulation.At intermediate frequencies there is a subtle amplitude modulation.This modulation is caused by the contribution from the center of the lens (a cusp in the Fermat potential), Eq. ( 26), and decays roughly as w −k 2−k .For the values k ≥ 1, shown in Fig. 6, the amplitude of the modulation is very small, and can only be appreciated for the SIS case (notice the slight difference between the WO and GO curves for the SIS case in Fig. 6 at around w ∼ 20).These cases converge to GO at a relatively small w (convergence to rGO is much faster).
In broad lenses, the envelope of the beating pattern persists at arbitrarily high frequencies.This envelope modulation is caused by the central image, and its amplitude is determined by the magnification µ H .The associated sizeable bGO contribution, discussed after Eq. ( 37), can be appreciated as a w-dependence of the modulation at intermediate frequencies (this is most appreciable in the k = 0.8 case in Fig. 6, where the amplitude of the modulation is largest at around w ∼ 20 and then decreases, converging to a constant at larger frequencies).The width of the modulation in frequency space is given by the difference in the time delay between the maximum and the saddle point, ∆w ∼ 1/(ϕ H − ϕ S ).As the impact parameter approaches the caustic, the amplitude and width of the modulation grow substantially (see the k = 0.65 case in Fig. 6).This regime is associated with a much slower convergence to GO.Some of these features are qualitatively similar to the cored isothermal sphere, as we will see below.

B. Cored isothermal sphere
We will now generalise the SIS by introducing a core with finite density [106,107], a feature shared by several DM scenarios that we will discuss in Sec.V.We will present the cored isothermal sphere (CIS) lens characteristics, its GO/bGO properties and WO signatures.

Mass profile and scales
The CIS is characterised by a core size r c and finite central density ρ 0 : This axisymmetric lens model reduces to the SIS for large radii r ≫ r c while providing a finite density at its centre.
The projected mass density is obtained as If we choose the normalisation scale ξ 0 to be then the lensing potential takes the form Here we have introduced a rescaled core radius x c ≡ r c /ξ 0 and added an unobservable constant factor x c to the lensing potential.
The virial radius for this lens is Contrary to the case of the gSIS lens, for the CIS the virial mass is just weakly affected by the new feature (the core radius in this case).Intuitively, for cores much smaller than the virial scale, r c ≪ r vir , the mass is dominated by the density away from the core.Hence, M vir must have the approximate form of Eq. (A10).To be more precise, the virial mass for the density profile ( 43) is Then, in the limit r vir ≫ r c , the second term on the right-hand side gives a constant negative contribution to the mass, ∼ 2π 2 r 3 c ρ 0 .From Eq. ( 47), we can notice that this limit is equivalent to the limit ρ 0 ≫ ρ vir .At first order in such limit, we can obtain ρ 0 as a function of M vir .Following the procedure used for the other lenses, we can then relate the effective lens mass M Lz to the virial mass M vir , (49) where we defined γ c ≡ 2(1 + z L ) 2 and the mass scale M 0 is given by Eq. ( 33).The second term in the square bracket of Eq. ( 49) represents a small correction from the result for the SIS (A11) due to the mass removed in the core.

GO structure & bGO features
Before discussing WO effects, let us examine the GO structure of the CIS.Depending on the core size and the impact parameter there can be one or three real solutions to the lens equation ( 20), similarly to the broad gSIS (k < 1, Sec.III A).As we will see below, multiple images can form only if x c < 1/2, i.e. if the density at the core is super-critical (Σ > Σ cr ).The GO images, their magnification and bGO corrections are shown in Fig. 7 as a function of the impact parameter for a CIS with x c = 0.15.
The critical curves are determined by the CIS deflection angle The tangential critical curve solves α(x)/x = 1 and is located at i.e. the core reduces the enclosed mass, pushing the Einstein ring inward (if the standard normalisation is kept).The tangential curve x tc is associated with a caustic at y tc = 0 when projected in the image plane via Eq.( 50).
The degeneracy of the caustic to a point is an artefact of the lens symmetry, any perturbation of the lens will spread the caustic [12].Another radial critical curve exists under the condition dα(x)/dx = 1, or which is associated with a caustic in the image plane Multiple images form only if with the two additional images corresponding to a saddle point and a maximum of the Fermat potential.The  53)).Middle: Magnification for the GW amplitude for the different images.The central image magnification varies between a finite minimum (∼ 0.429 for xc = 0.15) and the divergent value near the caustic.Bottom: Beyond geometric optics correction times the magnification for different images, obtained from Eq. ( 37) and Eq. ( 36).The central image has a large contribution that can overcome the geometric-optics magnification of the other two images for w ∼ 1.
known SIS result y tc → ±1 is recovered when x c → 0. In that case, the image associated with the maximum of the Fermat potential remains at the central cusp x H = 0 and is infinitely demagnified, µ H → 0.
The main GO effect of the core is to allow the central image to have a finite magnification µ 0 above a certain threshold.In the limit |x| ≪ x c , the central image asso-ciated to the local maximum lies at Note that this expression is valid for |y| ≪ 1/2 − x c , in addition to the strong-lensing conditions y < y rc , x c < 1/2.The magnification follows from Eq. ( 22) The first term on the right-hand side can be identified as a lower bound for the magnification (note that the second term is positive if the lens is super-critical, Eq. ( 54)).
The magnification increases monotonically with y until diverging at the caustic y → y tc , see Fig. 7.
The bGO factor |∆ H | |µ H | associated with the central image is similarly bounded from below and grows towards the caustic, where it diverges.For low impact parameters, it is smaller than the bGO factors associated with the minimum and saddle point, both diverging for y → 0. At intermediate y the bGO correction from the central image dominates over the other images.
There are several differences between the CIS and the broad gSIS (k < 1).The magnification of the central image approaches a constant µ 0 for aligned lenses y → 0, while it vanishes rapidly for the broad gSIS.In contrast, the bGO correction from the central image of the gSIS dominates over the entire strong-lensing regime, while for the CIS it only does at intermediate impact parameters.These features follow from the regularity of the CIS lens, compared to the cuspy gSIS.

Wave optics features
Let us now discuss the effect of the CIS lens parameters in the WO regime.Figure 8 shows the amplification factor at y = 0.3 for different values of the core size x c , between 0 (SIS) and 0.15, for which the source is right inside the caustic y rc = 0.3187 (Eq.53).At low w, x c has a moderate effect on the onset of magnification, with smaller cores associated with higher first peaks at low frequency.This is expected, as a larger x c decreases the enclosed mass near the lens center.
The beating pattern sets in at intermediate frequencies w ≳ 10.Its width is given by the time delays between GO images and its amplitude by the relative magnification.The primary, narrow oscillation is associated with interference between the minimum and the saddle point, ∆w N ∝ ϕ −1 S (recall ϕ L = 0).This is the only oscillation seen for the SIS at high w, although a damped oscillation on the envelope due to the contribution from the cusp (similar to the gSIS, Sec.III A and Eq.26) can be appreciated for w < 100.In contrast, the CIS has a smooth profile.
A finite core size introduces a secondary modulation in |F (w)|.This broad envelope is due to the central image, with a characteristic period of oscillation ∆w B ∝ (ϕ H − ϕ L ) −1 .Lowering the value of x c increases the width and height of this effect because the location of the caustic approaches the source (recall y is fixed).For the most extreme case shown (x c = 0.15), the magnification of the central image is sufficiently large to overtake the other two images.The central and right panels of Fig. 8 show the envelope of |F (w)| to highlight the secondary modulation (this is similar to the broad gSIS, cf.Fig. 6, top panel).
The proximity to the caustic determines the convergence to GO.This can be seen by the difference between the WO (solid) and GO (dashed) envelopes in Fig. 8.For WO and GO to agree well, a first requirement is to have negligible bGO corrections, w ≫ |µ I |∆ I for all images I, Eq. ( 24).In addition, it is necessary that w ≫ |ϕ I − ϕ J | −1 for all pairs of images with I ̸ = J.That is, the singular features in the time-domain integral are sufficiently far apart (relative to 1/w) to be resolved by the Fourier transform at a given w, cf.Sec.IV and Fig. 2 of Ref. [5].The low/intermediate frequency behaviour of the x c = 0.15 curve in Fig. 8 can be heuristically understood as the two images (maximum and saddle point) acting as a single GO feature and gradually becoming distinguishable at higher frequencies.
Let us now examine the effect of the impact parameter and the transition to the single-image regime.Figure 9 shows the amplification factor for x c = 0.05 and y between 0.1 and 2.5, with a value within and closely outside the caustic y rc = 0.581.The lowest impact parameters (y = 0.1, 0.3) show the same trends discussed above, with a narrow primary modulation and a broader secondary one due to the central image.The range of amplitudes |F (w)| for the more aligned system, due to larger magnification.The prediction for y = 0.55 ≲ y rc shows the same enhancement and broadening of the secondary modulation due to proximity to the caustic (but in the multiple-image regime) seen in Fig. 8.
The predictions are qualitatively different in the singleimage regime.At high frequencies there is only a single image, and therefore |F (w)| → |µ L | when w → ∞.For sizeable impact parameters, this limit is achieved at comparably low w.An additional feature is a weak pattern of the amplification at low frequencies 1 ≲ w ≲ 10 (e.g. at y = 2.5), which becomes more apparent for moderate impact parameters (y = 1.2).As the source approaches the caustic, y = 0.6 ≳ 0.581, this oscillation not only becomes very pronounced at low w, but it also extends to very high frequencies w ∼ O( 103 ).In this regime, one observes a primary beating pattern, but with a damped envelope (rather than a modulated one, as in the multiple-image regime).This is a WO feature, absent in GO.

C. Qualitative comparison and detectability
After discussing the lens models in the previous sections, we would like to assess whether the CIS and gSIS lens profiles can be distinguished from SIS by future GW detectors.We first discuss the information contained in the GO/bGO limit (Sec.III C 1).Then, we address the central images and bGO signatures (Sec.III C 2).Finally, we present the mismatch between the models and the SIS as a function of the lens parameters (Sec.III C 3).A quantitative analysis of the lens parameter reconstruction will be presented separately, Sec.IV.

Parameter reconstruction for w ≫ 1
Let us start our discussion by establishing how many lens parameters can be independently constrained for high-frequency sources.The GO amplification factor ( 21) is fully characterised by two continuous parameters per image: the magnification µ J and the time delay ∝ M Lz ϕ J . 7However, one can only measure differences in time delays (the overall Fermat potential is degenerate with the coalescence time) and magnification ratios (the overall amplitude is degenerate with the source distance).Therefore, in GO one can determine at best lens parameters, even in the absence of measurement uncertainties.A known consequence of this relation is that weakly-lensed signals, N images = 1, do not allow any lens parameter to be inferred.The situation is slightly better when including bGO corrections (24), as all additional ∆ I factors can in principle be measured.In that case, and one can in principle reconstruct N images additional lens parameters.Similar consideration can be applied to other 1/w corrections, such as the cusp contribution of the SIS and gSIS lenses, Eq. ( 26).
An additional consideration pertains to the precision with which the GO parameters can be reconstructed.In general, the time delays can be determined with precision 1/M Lz , while the magnification ratios uncertainty scales as 1/SNR [78].Ultimately, the quality of the different parameters affects the GO limit of the reconstruction, as we will explore in detail in Sec.IV B, IV C.
Addressing lens-parameter reconstruction in the full WO regime is not as straightforward.Non-perturbative WO corrections are given by a general function of w which cannot be captured in a finite number of parameters.This is in principle promising for reconstructing generic lensing potentials.However, important degeneracies still exist even for observations with SNR ≫ 1 (e.g. the projected matter distribution is a two-dimensional function, while the WO amplification factor is one-dimensional).In practical situations, we are mainly limited by the precision of GW observations, both in the WO and perturbative regimes.
Let us now focus on the lens models discussed above.The distinction between GO and bGO reconstruction is starkest for the narrow gSIS.As these lenses form two images but have three parameters, a reconstruction is only possible for intermediate-mass lenses, when bGO effects are relevant.This will be shown explicitly in Sec.IV C. Both broad gSIS and CIS lenses predict the existence of a third image, which qualitatively changes the parameter estimation in the high-w limit.Its detection allows us to reconstruct the three lens parameters even for very massive lenses, at least in principle, as we will see more in detail in Sec.IV B.

Properties of central images
Both CIS and broad gSIS lenses form a central image x H in the strong-lensing regime.This image is associated with the maximum of the Fermat potential, with |µ H | < 1.It is the closest to the lens' center.The central image has a large bGO correction (the largest if the lens is singular at the origin, as in the case of the gSIS, Figs. 4, 5).Since it is absent for the SIS (or narrow gSIS), its detection provides important information about the lens.
For generic values of the impact parameter (not too close to the critical curve), the central image is fainter than the other two.Depending on the lens parameters and the SNR of the unlensed GW, it might not be possible to detect it.For CIS, the magnification of the central image is given approximately by Eq. ( 56) and is controlled by x c .A similar reasoning applies to the broad gSIS lens with the parameter k < 1 (here the relevant parameter controlling the magnification of the third image is 1 − k).We first consider for what values of the lens parameters the central image has an SNR higher than a given detection threshold (which we call SNR th ).For LIGO, the typical threshold value is SNR th ≳ 8 for an event to be considered detectable.The detectability of the central image can be increased by considering sub-threshold triggers SNR < SNR th [73,74], but this increases the chance of a noise fluctuation or other event mimicking the signature [95].We take SNR = 30 for LIGO events (note that LVK events are volume-limited, hence typical SNR values for lensed events will be close to the detection threshold).Even in the case of LISA we consider the same threshold SNR th ≳ 8, although the typical SNR can be much higher than for LIGO.For LISA we take SNR = 10 3 .
Figure 10 shows the amplitude of the central image as a function of x c and k, for various values of the impact parameter.We can notice, in the CIS case, that x c ≳ 5 × 10 −3 is typically detectable with LISA signals and the dependence on y is minimal (unless we are very close to the critical y).This results from having a finite minimum magnification for the central image, as obtained in Eq. (56).On the other hand, for gSIS detectability requires deviations (1 − k) ≳ 0.2 and the de-pendence on y is more pronounced than for CIS.
As our analysis shows, central images are very sensitive to matter distribution near the center of the lens.In contrast, other images depend on the total projected mass, up to the radius where the image forms.This dependence is exact for axially-symmetric lenses in GO, as the reduced deflection angle α = M (x)/x, where M (x) is the total projected mass up to a radius x.Thus, finding central images of GWs might provide unique insights into the densest regions of galactic and DM halos [16,17,[20][21][22].An application of this idea will be the prospects for DM tests, cf.Sec.V.
We note that observing a central image might be more difficult for lensed EM sources.First, the EM signal is demagnified by |µ H | ≪ 1, which can be much smaller than the GW GO amplification |µ H |. Second, because of its central location within the halo, the EM signal might be blocked by dust or gas, or outshined by other sources.Finally, because of the high frequency of EM signals, it is highly unlikely that additional information from bGO terms can be retrieved.Thus, GWs provide a unique opportunity to probe central images and can be highly complementary to EM observations of lensed sources.See Sec.4.4 of Ref. [13] for a discussion of central-image searches with EM sources.

Mismatch in the WO regime
In this Subsection, we discuss how well one can distinguish a GW waveform lensed by either CIS or narrow gSIS (with k > 1) and a waveform lensed by an SIS.To assess the difference between the waveforms, we compute their mismatch, given a detector.For simplicity, we assume that both waveforms share the same parameters M Lz and y.The mismatch between two waveforms h 1 and h 2 is defined as Here we introduced the noise-weighted inner product.For two signals h 1 (t), h 2 (t) with Fourier transforms h1 (f ), h2 (f ), it is defined as where S n (f ) is the sky-averaged detector strain sensitivity.The signal-to-noise ratio (SNR) is defined in terms of this product as SNR 2 ≡ (h|h).The right-hand side of Eq. ( 59) is minimized over the phase and time difference between the two signals.Notice also that the mismatch is invariant under rescalings of either waveform; therefore it is automatically minimized over the luminosity distance of, say, h 1 (which is indeed just a rescaling).
A simple condition for detectability is that M ≥ 1/SNR 2 .This criterion is optimistic, as it neglects correlations between lens and source parameters that might make it more difficult to distinguish between lens models.Thus, the mismatch analysis is just a first step to addressing under which conditions two lenses can be distinguished.We will address degeneracies when performing forecasts, in Sec.IV. 8In Fig. 11 we provide few examples of the mismatches between a CIS and SIS (top row) and gSIS and SIS (bottom row) for impact parameters of y = 0.3 (left column) and y = 0.6 (right column).The mismatch is given as a function of both the lens mass (horizontal axis) and the additional lens parameter (x c or k−1).We have assumed a LISA equal-mass binary with mass M BBH = 10 6 M ⊙ at redshift z = 3.This type of source has SNR ∼ 10 3 , allowing the two waveforms to be distinguishable when the mismatch exceeds ∼ 10 −6 .See Sec.IV B for more details on the waveforms we use.
The merger frequency of the signals enters the GO regime for lens masses larger than M Lz ∼ 10 7 M ⊙ (that is, w ISCO ≫ 1, where ISCO refers to the innermost stable circular orbit): one can notice a change of behaviour of the mismatch in all the plots around this value.Increasing the impact parameter improves the mismatch slightly.This can be understood, at least in the GO regime, from the larger magnification of the third image for CIS.On the other hand, for gSIS the contribution of the cusp becomes larger.
In the case of the CIS, the mismatches are very sharp as a function of M Lz , with a milder dependence on x c .Thus, very small cores, even x c ∼ 10 −3 , can be distinguished for sufficiently large lens masses.This result will be also confirmed by the Fisher analysis in the next section.The case of the gSIS is less favourable for probing k − 1 close to zero, as the mismatch curves are flatter in M Lz .Then, the difference in slope k − 1 can be distinguished only up to ∼ 10 −2 for y = 0.3 (10 −3 for y = 0.6).

IV. PROBING LENS FEATURES WITH LIGO AND LISA
The existence of central images allows GW observations to probe the inner regions of matter halos.Together with the additional images and WO effects, they allow GWs to probe the matter distribution of gravitational lenses.In this Section, we perform a Fisher matrix forecast of the sensitivity expected from LISA and LIGO observations.Section IV A presents the frameworks and assumptions.We present the results for CIS in Sec.IV B and for gSIS in Sec.IV C. Appendices B, C, D and E provide further technical details.

A. Framework and assumptions
We will evaluate the sensitivity of different experiments to lens properties using a simplified Fisher matrix forecast [108], assuming Gaussian and stationary noise.The probability distribution of the model parameters θ i around their true values, θ i = θi + ∆θ i , can be estimated with a linearised-signal approximation (LSA).In the limit of high SNR, the leading contribution to the likelihood L only depends on the signal linearised around the fiducial values where the Fisher matrix F ij is defined as In the high-SNR limit, the standard deviations σ i for the parameters θ i are obtained from the inverse Fisher matrix as σ i = (F −1 ) (where no summation is implied).The Fisher matrix gives the unmarginalised error on a parameter, i.e. keeping all parameters fixed σ fix i = (F ii ) −1/2 .In conclusion, for high SNR the leading contribution to the σ i s is solely determined by the first derivatives of the signal.In this sense, the LSA is a valid approximation.
We consider the signal h to be given by the lensed waveform hL (f ) = F (w(f )) h0 (f ).This simplification assumes a single detector and neglects the time dependence of the antenna pattern.This is a good approximation if the signal's duration is much shorter than the timescale of the detector's motion This criterion is always satisfied for LVK signals, for which t signal can last as long as few minutes for neutronstar mergers, extending to tens of minutes for nextgeneration detectors.For heavier sources, signals can be as short as a fraction of a second.In the case of LISA, high-mass binaries can be observed months or years before coalescence.For this reason, we will therefore only consider the last month before the merger for LISA sources.The results change only minimally when considering the entire duration of the signal, as most of the SNR is located close to the merger [109].For LIGO, the information about the source sky location (from multiple detectors) can be used to correct for the antenna pattern.Therefore, we do not expect it to be a limiting factor in our analysis.Strictly speaking, the staticity requirement also sets a limit on the time delay between images, and thus on the effective lens mass Given the short-signal condition (63), violations of the above relation imply observations in the GO regime, since t IJ ≳ t motion ≫ t signal ≫ 1/f .In practice, Eq. ( 64) is necessary to ensure that the different antenna pattern is approximately constant for all GO images.This is a very stringent requirement.In practice we expect our results to be valid as long as all images can be detected and associated (this will be the case if the SNR of the different images is comparable).In this sense, gaps in observation, finite survey time [32] and false-positive image association [95] might be more stringent requirements.We will further assume an edge-on source aligned with the detector, so h0 = h+ .Because of this simplified setup, we will neglect the error on the source orientation and sky localization, whose determination requires including rotating antenna patterns or multiple detectors.Regarding the source parameters, we will marginalize over the coalescence phase ϕ 0 and luminosity distance from the source D S = (1+z S ) 2 D S .We will assume equalmass ratio, non-spinning sources, with all other parameters fixed.This follows the approach in Ref. [4], which focused on point-lens and SIS.A detailed analysis of the same lenses including detector motion, source properties and waveform models found that the results are robust against additional source parameters, and might even be more optimistic (e.g. by the introduction of higher harmonics) [36].
We will normalize our results to a fiducial SNR of the lensed source (10 2 for LIGO, 10 3 for LISA).The actual errors scale as σ i = σ i fid SNR fid /SNR.We expect some deviations from our assumptions to be approximately captured by a reduced SNR: for instance, a sub-optimal sky localization or a different source orientation.In our analysis, we will be primarily working with high-SNR signals.As we will see, however, the validity of the linearisation of Eq. ( 61) is not always granted and needs to be explicitly checked.We will come back to this point below and when interpreting the results.
We compute the amplification factor and its derivatives using three methods, appropriate for different regimes.At high frequencies, w > w high ≳ 100, we use the bGO approximation summarised in Sec.II B. In the case of the gSIS lens, on top of the bGO contributions, we add the terms associated with the cusp, summarised in Eq. (26) (see also Sec.III-A-2 of Ref. [5]).For intermediate frequencies w low < w < w high we use instead the contour method summarised in Sec.A 2. Derivatives of F (w) with respect to the lens parameters are computed using the prescription defined in Appendix A of Ref. [5].For frequencies below w low = 0.05, we extrapolate using the analytic dependence described in Sec.A 1. The transition frequencies are defined so that the errors between methods are below 1% for the amplification factor and 3% for its derivatives with respect to the lens parameters, see Fig. 24.More details about the implementation of the Fisher matrix are given in App.B.
The inversion of the Fisher matrix is very sensitive to numerical errors from the FFT truncation at high w.To test the high-M Lz limit calculation, we developed the Geometric Optics Diagonal Approximation (GODA) of F ij , as well as its extension to bGO and rGO (the latter defined as including both bGO and cusp contributions).This approximation, described in App.C, greatly speeds up the computation of the Fisher matrix by neglecting contributions that oscillate in frequency space which originate from products of amplification factors corresponding to different images F * I F J .The result is a sum over images of terms in |F I | 2 (hence diagonal in image space) times the unlensed SNR weighted by different powers of the frequency, such as (h 0 |h 0 ) and (wh 0 |wh 0 ).
At sufficiently high M Lz we expect the LSA approximation to break down.We diagnose this by evaluating the phase difference in the GO limit Here J labels the images while θ i are the lens parameters.The reason behind the LSA breakdown can be understood as follows.In the lensed waveform h L , the lensing parameters appear in the magnification and time delays.The latter, in particular, are multiplied by w in the phases of Eq. (24).Linearisation is a good approximation if the changes in the phases are small.However, this is not necessarily the case for large w, or equivalently large M Lz .
We will indicate the values of the masses at which ∆ LSA > 1 in our plots (we check for this condition in the GO regime, w > 100).We evaluate Eq. ( 65) using Eq. ( 15) to write w as a function of f and M Lz , the latter displaced 1σ from the fiducial value (this is the parameter that gives the largest contribution to ∆ LSA when varied from its fiducial value).The other lens parameters θ i are taken as fiducial.For the frequency we choose f * , the value such that 90% of the total SNR is given for f ≤ f * .
The image giving the largest contribution to ∆ LSA is the central one (having the largest time delay).Therefore, the expression we use for ∆ LSA reduces to where ∆ log(M Lz ) is the 1σ uncertainty on log(M Lz ), as obtained from the Fisher analysis.For the gSIS lens, with k > 1, ϕ H is replaced by the time delay of the lens center (the location of the lens' cusp).Above these lens masses, the breakdown of LSA indicates that the Gaussian approximation for likelihood does not hold, and the results need to be taken with caution.Nonetheless, whether the results are over or underestimated will depend on the lens model: for the CIS (3 images) the Fisher approximation is correct, while for the narrow gSIS (2 images + 1 cusp) it leads to an overestimation of the sensitivity at high M Lz .This result is consistent with the lack of available GO information (relative time delays, magnification ratios) to constrain all the lens parameters, as discussed in Sec.III C 1. We will comment on the case of each lens in the corresponding Sections, with more details on Appendices D, E.
Finally, we note that the LSA validity condition (66) depends on the SNR of the signal (through the uncertainty ∆ log(M Lz ), that scales as ∼ 1/SNR).To evaluate this condition in our results, we use the fiducial SNR of the given experiment (10 3 for LISA and 10 2 for LIGO).Events with lower SNR will violate the LSA at lower lens masses.

B. Cored Isothermal Sphere
Let us now consider the capacity of different experiments to measure the lens parameters in the CIS lens.We will focus on the strong-lensing regime, with multiple GO images.We will present results for different source masses, impact parameters and fiducial slopes.We show the marginalised 68% c.l. posterior in the Fisher-matrix approximation as a function of the effective lens mass M Lz , which can be obtained without recomputing F (w).Our fiducial core size is x c = 0.01.
We start by considering LISA sources.We explore equal-mass, non-spinning binaries at z = 3, with total mass 10 4 − 10 8 M ⊙ (source frame), using IMRPhenomD waveform [110,111].Figure 12 shows the lensed and unlensed strains along with the expected LISA sensitivity for M Lz = 10 7 M ⊙ .The figure highlights the last month of data before the merger, the only portion of the  signal we included in the analysis, cf.Eq. ( 63).We verified that the constraints vary negligibly when including longer data, except for the lightest sources 10 4 M ⊙ for which the merger occurs off the LISA band.The second panel shows the lensed SNR for those sources as a function of M Lz .The results of the marginalised errors will be rescaled to a fiducial SNR of 10 3 .This factors out the dependence on the distance exactly and sky localization and source inclination approximately.
Figure 13 shows the expected 68% 1D marginalised posteriors for the lens parameters log(M Lz ), y and x c .As just explained, the results are rescaled to a fiducial lensed SNR of 10 3 : the precision of a given source needs to be rescaled by 10 3 /SNR, e.g. with the lensed SNR given in the right panel of Fig 12 .At fixed SNR and for highmass lenses, the precision of all parameters converges to a constant value, of order ∼ 1/SNR and independent of the source mass.The saturation of the sensitivity in the GO limit will be explained below in terms of parameter degeneracies (Figs.16,17). 9ighter sources are more effective at probing lower lens masses, at fixed SNR and when the signal is dominated by the merger.This is a consequence of the higher GW frequencies, emitted up to ∼ w ISCO ∝ M Lz /M BBH so the onset of magnification w ≳ 1 corresponds to lower M Lz .In this case, the lens parameters can be independently determined for M Lz ≳ 0.1 M BBH (10 3 /SNR).For the lightest sources shown (M BBH = 10 4 , 10 5 M ⊙ ) the SNR is not dominated by the merger, but by the waveform portion in the "sweet spot" of the LISA noise curve (cf.Fig. 12).In these cases, the higher frequencies of the source do not play a role, which explains their similar  shape in Fig. 13.Note that the effect of different SNRs needs to be included when considering specific sources.For instance, a 10 6 M ⊙ binary has a lensed SNR ∼ 40 times larger than that of a 10 4 M ⊙ source.In what follows we will show results for lensed 10 6 M ⊙ LISA sources, as they are the most favourable case.
Varying the impact parameter changes the precision in the high M Lz limit, with higher precision for larger y, for fixed SNR.This is shown in Fig. 14 for y = 0.2, 0.3, 0.5, 0.6, sufficiently away from the caustic to avoid numerical problems when computing derivatives in the amplification factor, cf.App.B. The increase in precision when y is increased is larger for log(M Lz ) and x c .This is explained by the lens mass being probed by time delays, which are larger for less aligned configurations.Similarly, the imprint of the core is more obvious for larger y, as both the amplitude of the central image and its bGO corrections increase as the source approaches the caustic, cf.Fig. 7 (for our fiducial model y rc ≃ 0.805, Eq. ( 53)).This gain in sensitivity is partially offset when including the actual SNR, which is higher when the source and lens are aligned, at low y.
Varying the core size at fixed y has a small impact on the constraints unless x c becomes large enough for y to start approaching the caustic y rc .Figure 15 shows marginalised 1D posteriors for x c = 10 −4 , 10 −3 , 0.01 and 0.1, for fixed SNR and y = 0.3.In these cases the caustic Eq. ( 53) is located at y rc = 0.98, 0.937,  0.805, 0.427, respectively, with a slow convergence to SIS (x c → 0, y rc → 1) due to the √ x c term in Eq. (53).All forecast curves show qualitatively similar dependence on the lens mass, except for the larger cores x c = 0.1, for which the impact parameter approaches the caustic and the constraints are slightly enhanced.The lens parameters log(M Lz ), y, x c follow very similar trends regarding the dependence with source and lens mass (Fig. 13), as well as with the impact parameter (Fig. 14) and core size (Fig. 15).In contrast, the FIG.18. Same as Figs. 13 and 14, but for advanced LIGO Sources (see text).In the case of LIGO, the LSA condition ( 66) is typically violated as soon as we enter the GO regime given the larger uncertainties compared to the LISA case.Nonetheless, as we motivate in the main text, for the CIS lens we expect our results to be still reliable, as there are enough GO parameters to reconstruct all the lens parameters.
source parameters (log(D S ), ϕ 0 ) can always be determined, although their precision suffers at intermediate masses (cf.lower panels in Fig. 14) due to degeneracies with the lens parameters.
Let us now address the correlations between different parameters.Our previous discussion (Figs.13,14,15) focused on the 1D marginalised posteriors, which quantify the ability to constrain the parameters separately.All these plots showed how the 1D posteriors saturate at high M Lz at fixed SNR.This can be explained by the amount of GO parameters available to reconstruct the lens (Sec.III C 1): two relative time delays, whose precision improves as ∼ 1/M Lz and two magnification ratios, whose precision remains ∼ 1/SNR [78].Thus, one combination of log(M Lz ), y and x c remains poorly constrained (for fixed SNR but increasing M Lz ).This trend can be observed explicitly in Fig. 16, representing the 2D marginalised posteriors in the Fisher-matrix approximation.The bottom-left panels show the posteriors shrinking as M Lz increases, for moderate lens masses.For high M Lz (top right panels) the ellipses become narrower but otherwise span the same range of lens parameters.
To further understand the degeneracies it is instructive to look at the principal components of the Fisher matrix, i.e. its eigenvectors ⃗ u (k) and their associated eigenvalues Each ⃗ u (k) represents a combination of parameters that can be constrained independently (i.e. the eigenvectors are normal ⃗ u . The left panel of Fig. 17 shows the eigenvalues and eigenvectors of F ij , ranked by precision, as a function of the lens mass.We note two distinctive regimes for light (M Lz ≲ 10 6 M ⊙ ) and heavy (M Lz ≳ 10 8 M ⊙ ) lenses, with an intermediate trend for M Lz = 10 7 M ⊙ .For lighter lenses (10 4 −10 6 M ⊙ ) the two most precise combinations are dominated by the source parameters log(D S ), ϕ 0 , whose precision saturates for light lenses.The remaining three combinations are mainly combinations of lens parameters, whose precision decreases substantially as M Lz → 0, explaining the loss of sensitivity seen for light lenses in Fig. 13, and expected from the convergence to F (w) → 1 at low w.
For heavy lenses, M Lz ≳ 10 8 M ⊙ , Fig. 17 shows how the precision of the two best-measured combinations increases with M Lz , while the other three combinations remain independent of the lens mass.1D marginalised posteriors are a combination of all eigenvalues and hence dominated by the least-precise ones, which saturate for M Lz /M BBH ≳ 10 2 .The higher precision in the first two eigenvalues for larger M Lz appears in parameter correlations as a thinning of the posteriors, seen in Fig. 17.The large degeneracy would allow for substantial improvement in the constraints if one or more of the lensing parameters can be measured independently, for instance if the lens can be located via an EM counterpart or crosscorrelating with optical catalogues [112].
Let us now consider the capacity of ground detectors to characterize CIS lenses, focusing on the advanced-LIGO sensitivity curve at design values [113].We will consider equal-mass non-spinning binaries with total mass M BBH = 3, 15, 30, 60 and 120 solar masses observed respectively for 1 h, 6 min, 1 min, 1 min and 10 s before merger.Reasonable variations of the observing time do not affect the final results.We will show results up to M Lz < 10 6 M ⊙ .For the lens parameters we consider (e.g.away from caustics), this is sufficient to reach the GO saturation described above.This choice of masses also justifies considering a static detector, as typical time delays between GO images scale as ∼ 4GM Lz = 19.4s M Lz /10 5 M ⊙ .As the time delays remain much smaller than changes in detector orientation, on the scale of 1 day = 8.6 • 10 4 s, our results could be extended to slightly higher lens masses.
Figure 18 summarises the results for equal-mass nonspinning sources at z = 0.1 as detected by one LIGO detector.In this case, the results are normalized to a fiducial lensed SNR of 100, which is plausible for lensed sources at low redshift (see the top-right panel in Fig. 18).The trends are qualitatively similar to LISA detections, shifted towards a lower range of M Lz , as expected from the scaling of w, Eq. ( 15).Note that the curves are closer together because the hierarchy of M BBH considered spans fewer orders of magnitude than in the LISA case.The uncertainties on individual parameters are larger, but this is mostly due to the choice of fiducial lensed SNR.
To test the validity of the CIS Fisher matrix, we compared the full log-likelihood (D1) to the quadratic expansion (61), focusing on the GO limit for the LISA case for simplicity (see App. D for details on how these quantities are obtained).
We sampled the full likelihood at 10 3 random points located at 2σ from the fiducial, as defined by the Fishermatrix eigenvectors.We find that in 99.7% of the cases, the full likelihood (evaluated at 2σ) is smaller than the Fisher likelihood (evaluated at 1σ).This shows that the true 1σ contours (as defined by the full likelihood) are bounded to be within the 2σ contours of the Fisher, indicating that the Fisher approximation is an adequate description, i.e. it captures the degeneracy directions and magnitude correctly.

C. Narrow generalised-SIS
Let us now consider the capacity of different detectors to measure the slope in the matter distribution using the gSIS lens.We will focus on narrow gSIS (k > 1) with two GO images, as the broad gSIS is qualitatively similar to the CIS lens (Sec.IV B).Having fewer images, its analysis is qualitatively different.Our fiducial value of k is very close to the SIS, k = 1 + 10 −4 , indistinguishable from the SIS for reasonable values of the lensed SNR (≲ 10 4 ). 10efore discussing our results, we note that the LSA breakdown leads to artificially optimistic results in the full WO calculation at high M Lz .For the narrow gSIS, only two GO images exist, not providing enough information to recover M Lz , y and k, cf.Sec.III C 1 (unlike for CIS, where the existence of three images allows a full reconstruction in GO).Thus, information from bGO corrections ∝ ∆ I /w and the cusp F c ∝ 1/w plays a critical role in the narrow gSIS.While these terms become negligible at high M Lz , the WO Fisher-matrix results predict a constant sensitivity.
The breakdown of the LSA is due to the cusp contribution.It can be traced back to the linearisation of the phase e iwϕc(k) in F c (w): if linearly expanded around k = k+∆k, it gives ∼ iw∂ k ϕ c ∆ke iwϕc( k) .Clearly, the linearisation is invalid if M Lz becomes too large and moreover, if trusted beyond this limit, this factor becomes much larger than the phase e iwϕc(k) .Therefore, the linearisation seems to overestimate the sensitivity on k.In order to confirm our intuition, we have checked that indeed the contribution of the cusp to the Fisher-matrix log-likelihood (given by Eq. ( 61)) becomes much larger than the corresponding contribution in the full likelihood (i.e.not expanding in powers of ∆θ i ), for some sample points at high lens masses.We expand on the issue in App.E with an analytical toy example and show different predictions (Fig. 26) that confirm the disproportionate effect of the cusp at high M Lz .
To circumvent the cusp problem without abandoning the Fisher-matrix analysis, we will show two results for each case: the full WO Fisher forecast and the bGODA approximation.Since the former overestimates the contribution from the cusp and the latter neglects it, the truth should lie between both.In addition, we highlight the values of M Lz where the LSA approximation, Eq. ( 65), breaks down.In each plot, we show this for the model where the breakdown happens at a lower lens mass.
Given the premise above, we can now briefly discuss the results of our forecast, obtained using a combination of WO and bGODA, for LISA.As in the CIS case, we conventionally give results normalized to a fiducial SNR = 10 3 .See the discussion in Sec.IV A on how to rescale the results for a given source.
Figure 19 shows the dependence on y, at fixed k for  65) is violated.This criterion is obtained as ∆LSA < 1 and also by requiring we are in the GO limit (wmax > 100).Dashed lines are the results of the Fisher forecast using the bGODA approximation which neglects the cusp contribution.The correct constraint should lie between WO (which overestimates the contribution from the cusp) and bGODA (which neglects it completely). .The grey area is the region where the LSA condition is violated and the dashed lines are the Fisher forecast results in the bGODA approximation (see caption of Fig. 19).The solid dark-green lines (k = 1.01) are expected to start growing at large lens masses, following the qualitative trend of the light-green curves (k = 1.3).This ultimately follows from the scaling of the cusp contribution Fc(w) as a function of w (see Eq. ( 26)).
two fiducial values of y = 0.2 and y = 0.6.Larger values of y give worse ∆y and better ∆k, ∆ log D S instead.
From this Figure, one also notices how ϕ 0 can always be recovered at high M Lz , contrary to the recovery of lens parameters.These trends follow the CIS case already described.Interestingly, the inclusion of WO effects breaks the parameter degeneracy, and the slope k can therefore be measured, both in WO and bGO.The best accuracies for the lens parameters are expected when the on-set of WO, w ∼ 10, happens around the peak of the signal SNR.Typically, this corresponds with the ISCO frequency, or w ∼ 8πGM Lz f ISCO ∼ M Lz /M BBH ∼ 10.
On the other hand, the results degrade when M Lz is increased or decreased (if we take the worst between the WO and bGODA results, which is conservative).Fig. 20 shows the dependence on the fiducial value of k, at fixed y.We can notice slight differences at the onset of WO, around M Lz ∼ 10 7 M ⊙ : larger k gives better results.For higher lens masses instead, the noticeable difference is in ∆k, whose 1D posterior worsens for larger ks (and for the other lens parameters but only in the region where the LSA fails).This trend is explained by the different scaling of the cusp corrections ∝ w −k/ (2−k) .Due to this scaling, the cusp contribution decays faster for larger ks and does not lead to a saturation of the Fisher posteriors.

V. TOWARDS TESTS OF LARGE-SCALE STRUCTURE AND DARK MATTER
We will consider the potential of GW lensing to constrain dark matter (DM) properties.First, we will discuss the prospects of inferring the virial mass of the lens, given that its redshift is unknown, Sec.V A. Then we consider two scenarios that lead to the formation of cores: modelling DM halos as a CIS allows us to translate our forecasts into a range of DM parameters accessible to LIGO and LISA.We discuss two theories proposed in the context of small-scale challenges to standard cosmology [71]: self-interacting DM (SIDM, Sec.V B) and ultralight DM (ULDM, Sec.V C).We conclude by discussing the assumption of axially-symmetric lenses, Sec.V D.

A. Lens redshift and halo mass
Let us now address the prospect of relating the measured variables (M Lz , y, x c , D S ) to the physical properties of the lens.The main obstruction is the lack of knowledge of the lens' redshift z L , which is needed to determine the virial mass of the halo (Eqs.(A11), ( 32) and ( 49)) and the physical scales, i.e. r c = x c ξ 0 (M Lz , z L , z S ) (Eqs. ( 2) and ( 16)).The dependence of physical scales (virial mass) on the lens redshift is shown on panel A (B) of Fig. 21: the maximum (minimum) value is reached at an intermediate redshift, while for z L → 0 and z L → z S the curves approach zero (infinity).Values several orders of magnitude away from the maximum (minimum) require lenses that are very close to either the source or the observer.While not impossible, this situation is unlikely.
We can now compute the probability distribution for z L and M vir , given M Lz , z S and some reasonable assumptions about the cosmology and the halo mass function.The probability distribution for M vir and z L describing the lens is dP dz L dMvir ∝ σ dV dz L dn dMvir .The angular crosssection is σ = 2πy∆y (ξ 0 /D L ) 2 , i.e. the projected area of the annulus compatible with the measured impact parameter.The volume element dV dz L = and halo mass function are both in comoving coordinates.We will use the Tinker et al. form of dn/dM vir [114], as implemented in the Colossus package [115].We will assume fixed cosmological parameters, given by Planck ΛCDM best fit [116] and focus on an SIS when concrete results are needed.Changes due to the cosmology or halo mass function do not qualitatively affect our analysis. 11e will focus on z L and M vir , assuming that M Lz , y and z S are known.The former can be inferred precisely, while the latter can be constrained from the source's luminosity distance D s assuming an expansion history (including uncertainties on these and other parameters,e.g. from parameter posteriors is straightforward).Then, the probability of a certain lens redshift is where in the second line we omit all the constant terms.
Here MLz is the measured value, and a bar means that a quantity is evaluated on it, e.g.ξ0 = ξ 0 ( MLz , z L , z S ) in Eq. ( 16).The jacobian stems from integrating the delta function over M vir and is ∝ M A similar calculation allows us to find a probability distribution for the lens' virial mass.Starting from the r.h.s. of Eq. ( 68) but integrating over z L gives (70) Here k = 1, 2 labels the two solutions of M vir (z L , MLz ), see panel A in Fig. 21.The probability distribution for the virial mass, Eq. ( 70), is shown in panel D of Fig. 21.The probability is rather peaked at the lowest value due to the divergence of the jacobian at the minimum of M vir .The probability of M vir being much larger than the minimum is very suppressed: it requires z L → 0, z S , a limit in which the cross section (∝ ξ 2 0 ) vanishes very rapidly (see Fig. 21, A, B).
This framework allows one to constrain the virial mass to within a factor of a few.For M Lz = 10 7 M ⊙ , the 95% c.l. limits are M vir ∈ [3.8, 11) 10 9 M ⊙ for z S = 1 and [2.3, 7.1) 10 9 M ⊙ at z S = 3.The distribution becomes wider for higher masses and redshifts: for z S = 10 the 95% c.l. ranges are [1.7, 17)  M Lz , z S and other parameters will lead to an uncertainty on this lower bound, as well as a broadening of the probability density shown in panel D of Fig. 21.Nonetheless, given the precision in M Lz , D S (e.g.Fig. 14), M vir 's upper limit uncertainty will be dominated by lack of knowledge of z L .
We note that the limitations discussed above may be lifted if the source/lens system can be accurately identified, e.g.via EM follow-ups [112,117].

B. Self-interacting dark matter
SIDM has been proposed to solve the core-cusp and missing satellites problems [69].In this scenario, DM particles scatter elastically with each other, leading to deviations from CDM predictions regarding the inner halo structure [70,[118][119][120][121][122].The self-interaction is described by a cross-section, which is in general velocity-dependent.For a recent summary of constraints and measurement claims see Table I of Ref. [70].
We will estimate the size of the core through the condition [123]: Here ρ(r c ) is the density evaluated at the core radius r c , ⟨σv⟩ is the velocity-averaged self-interaction crosssection, m is the mass of the DM particle and t age is the time since the formation of a given structure.Equation (71) establishes that the core forms where the DM density ρ is high enough for DM particles to interact on average at least once since the formation of the halo.We will take t age = 5 Gyr as our fiducial choice.This corresponds to the age of clusters of galaxies.It is a conservative choice, as lighter structures will have higher t age , leading to larger cores and more stringent constraints.Let us first give an order-of-magnitude estimate for the constraints on the DM cross-section ⟨σv⟩/m expected from a detection of a lensed GW event.At large-enough radii, the density of DM halos is typically described by an NFW profile, ρ NFW (r) ≡ 4ρ s (r/r s ) −1 (1 + r/r s ) −2 for r > r c .Self-interactions can instead lead to a constant density inside the core, ρ(r) ≃ ρ NFW (r c ) ≃ 4ρ s /(r c /r s ) for r < r c (with r c ≪ r s ).In this case, the core size can be characterised by the dimensionless variable xc ≡ r c /r s (analogous to x c for the CIS lens discussed previously).Using the relation (71) Here we replaced ρ s using its phenomenological relation with the virial mass (see e.g.Eqs. ( 42) and (43) of Ref. [51]).This relation can be turned into a bound for ⟨σv⟩/m.Indeed, using Eq. ( 72), a bound from GW lensing on the dimensionless core size xc ≲ O(1) for virial masses M vir ≃ 10 10 M ⊙ would constrain the averaged cross section to be ⟨σv⟩/m ≲ 5 • 10 −4 cm 2 /g.In this estimate we made the assumption of having xc ≲ O(1) at M vir ≃ 10 10 M ⊙ .This is expected to hold, as virial masses in this range typically produce lensing signals in the GO limit, where lensing parameters are obtained accurately.More sophisticated lensing forecasts, for a cored NFW profile, are needed in order to make this estimate more precise.
Assuming instead that the matter density is described by a CIS profile, we can outline how to obtain concrete bounds using our forecast results.In this case, the projected core size x c ≡ r c /ξ 0 is given by . (73) The lens redshift is chosen to minimize the projected core size (see text).Constraints on the cross section scale as ∝ ρ(rc), Eq. ( 71).Therefore, results for an NFW profile are expected to scale as the square root of the CIS ones, compare Eq. ( 72) and ( 73).The forecasted CIS results are orders of magnitude tighter than existing constraints [70].Ultimately, probes based on very small cores (low ⟨σv⟩/m) will be limited by astrophysical uncertainties.
This expression has been obtained by evaluating Eq. ( 71) on the CIS density profile, substituting ξ 0 from Eq. ( 15) in terms of the virial mass (49), where we neglect terms ∼ r c /r vir ≪ 1.Note that the dependence between the projected core size and the cross-section is quadratic, as given by the CIS dependence.Therefore, the capacity to probe very small cores grants access to very low SIDM cross-sections.This is in contrast to the NFW case, where this relation is only linear, Eq. ( 72).The different behaviour originates from the different slopes of the two profiles in the region r c ≲ r < r s .We will assume the lens' redshift that maximizes (1 + z L )d eff , for fixed z S .This choice turns the equality in Eq. ( 73) into a lower bound on the projected core size and is a pessimistic assumption in terms of using x c to constrain SIDM.We will show our results as a function of M vir , obtained from M Lz with this choice of d eff .The choice of z L results in a plausible value for the virial mass given the uncertainties, see Fig. 21.
Figure 22 shows the values of the SIDM cross-section accessible to a typical lensed LIGO and LISA source (y = 0.6), as a function of the halo mass for the CIS profile.The curves are derived by interpreting the 95% marginalised posterior (2∆x c , cf.Sec.IV B) as the minimum core size that can be probed, and using Eq.(73) to relate it to ⟨σv⟩/m.This interpretation is supported by the convergence of the posteriors as x c → 0, Fig. 15.We limit the constraints to the region where x c ≤ 1, as no multiple images occur otherwise, Eq. ( 54).Self-Interacting DM may still be probed in that regime, but a separate analysis needs to be performed.
The capacity of lensed GWs to probe SIDM is optimal at the onset of magnification w ISCO ∼ 10.While ∆x c is smaller at higher M Lz , lighter lenses have smaller ξ 0 , making small cores appear larger.At masses below the onset of magnification, the decrease in sensitivity offsets the projection effect.Interestingly, under our conservative assumptions for d eff , the optimal LIGO sensitivity corresponds to halos dominated by DM, with very low baryonic mass.For instance, for ultra-faint dwarfs, M vir ∼ 10 9 M ⊙ , the stellar-component mass is M ⋆ ∼ 10 −5 M vir (see Fig. 6 in Ref. [71]).While the uncertainty on z L implies that M vir can not be precisely determined, the probabilistic treatment discussed in Sec.V A can be used to place plausible limits.
The above estimates rely on several simplifications that need to be revised to obtain accurate constraints on the properties of DM.We have used the CIS lens model, which might fail to capture further details of the lens, such as the density profile outside the core (different from NFW).The fact that CIS grows as 1/r 2 at small r (but outside of the core) gives very strong constraints while shallower inner profiles, such as NFW below the scaling radius, will plausibly lead to milder constraints, as seen by the estimate in Eq. ( 72).Another necessary refinement is the inclusion of stellar/gas components on top of the core.This will be most important for large galactic halos, where cored profiles can arise due to baryonic effects (cf.Fig. 13 of Ref. [71]). 12This contribution may be important even for low M vir when probing low cross sections, and hence very small cores.Eventually, realistic analyses should extend the probabilistic treatment of Sec.V A and consider non-axisymmetric lens models.We briefly discuss these extensions in Sec.V D.
C. Ultra-light dark matter ULDM refers to models where DM is a bosonic particle with a very low mass m ϕ , typically in the range 10 −22 eV ≲ m ϕ ≲ 1 eV [67,68].These models exhibit interesting phenomenology on small scales, relative to the field's de Broglie wavelength λ ϕ ≃ 1.92 kpc 10 −22 eV/m ϕ ((10 km/s)/v) [66].These fields appear generically in extensions of the Standard Model of particle physics e.g. to solve the strong-CP problem (QCD axion) [128,129].In this Section, we will consider non-interacting, scalar ULDM.The shaded regions show ULDM masses accessible to lensed GW in which a core can be excluded, Eq. ( 76).We set y = 0.6 and zS = 3, MBBH = 10 6 M⊙, SNR = 10 3 for LISA, zS = 0.2, MBBH = 30 M⊙, SNR = 10 2 for LIGO.The lens redshift is chosen to minimize the projected core size (see text).The green region represents comparable limits from ultra-compact dwarf galaxies [131].
ULDM predict cored DM halos with a minimum size and a maximum central density.Here M c is the mass of the solitonic core [66], and equality holds when the solitonic core forms.The solitonic-core mass is related to the virial mass by (see Eq. ( 6) in Ref. [130], neglecting a redshift factor ∝ (ζ(z)/ζ(0)) 1/6 ). 13  Assuming that DM halos can be modelled by the cored isothermal sphere with x c ≃ r 1/2 /ξ 0 results in a minimum projected core size (76) This expression follows from combining Eqs. ( 74) and (75), expressing ξ 0 = √ 4GM Lz d eff , Eq. ( 15), also in terms of the virial mass for the CIS (49).We will assume the lens redshift is such that (1 + z L ) 3/2 d eff is maximized, for fixed z S (see discussion after Eq. ( 73) and in 13 Ultra-light scalars also predict a minimum M vir for DM halos (Eq.( 42) in Ref. [66]).Strongly-lensed GWs constrain M Lz , and thus can bound M vir , following Eq.( 70).This test is less powerful than constraining the size of the core (although less model-dependent) and we will not consider it further.
Sec. V A).This gives a lower bound for the inequality in Eq. ( 76), and hence pessimistic results.
Figure 23 shows the range of DM masses that produce cores detectable by LIGO and LISA.Masses in the shaded regions and below can potentially be excluded by observations of strongly-lensed GWs (we have taken 2σ exclusion in ∆x c , see discussion on SIDM, Sec.V B).Again, we terminate the plots at low masses, where ∆x c > 0.5, i.e. when no multiple images can form.Note that the sensitivity is well above the m ϕ ∼ 10 −22 eV value for a broad range of halo masses, potentially probing fields with λ ϕ ≳ O(1 kpc).
We expect our modelling of the core to be accurate for halos below a certain mass: cosmological simulations of ULDM show how low mass halos (M vir ≲ 10 8 M ⊙ ) flatten below a certain radius (Fig. 13 of Ref. [132]), consistently with our CIS modelling.In contrast, heavier halos develop dense and steep cores ρ ∝ (1 + (r/r c ) 2 ) −8 [130], not captured by the CIS parametrisation.More massive halos will also have a larger baryonic component contributing to the profile (Fig. 4 of Ref. [133]).Outside the core, simulations predict an NFW profile with a different slope.This difference will be most relevant for large impact parameters and less so for the central image.Other aspects of the analysis need to be refined, as already discussed in the last paragraph of Sec.V B.

D. Non-axisymmetric lens profiles
Let us now comment on how our results extend to nonaxisymmetric lenses.To keep our discussion simple, we will focus on elliptical lenses, i.e. replacing the dependence on x → x 2 1 + x 2 2 /q 2 [134].This introduces two parameters, the ellipticity q and the angle between the lens' axis and source δ = arctan y 2 /y 1 .Ellipticity can be quite substantial (i.e. the axis ratio of the ellipse can be few tens of percent) and degenerate with other deviations from spherical symmetry [13,135].Elliptical lenses can form additional GO images in certain regions of the source plane [134].
Breaking axial symmetry will introduce additional degeneracies.For example, an elliptic generalization of the SIS [134] at moderate impact parameter will form only two images.Then, in the GO limit only two combinations of M Lz , y, q and δ can be constrained, cf.Eq. (57).A CIS (or broad gSIS) in the same regime forms 3 images, allowing to constrain 4 combinations of M Lz , y, q, δ and x c but leaving one combination unconstrained.
While full reconstruction of ellipsoidal lenses is compromised in the 2-3 image GO regime, WO effects greatly improve the prospects.In the bGO approximation, an event with 3 images allows the inference of up to 7 parameters, Eq. ( 58).This is enough to fully constrain the ellipsoidal generalization of the CIS or broad gSIS (5 parameters).We expect the parameter reconstruction to improve in the full WO regime, in which the entire extent of the lens contributes to the amplification factor [6].
Information based on central images is robust to ellipticity and can lead to constraints even in the GO limit.
The existence of a central image is due to the lens' profile and is not affected by ellipticity or lens orientation.Detection of a central image (i.e. a system with 3 GO images) thus sets limits that are independent of q, δ, such as a lower bound on the lens' core size (CIS), or an upper bound on the density slope (broad gSIS), see Fig. 10.Conversely, non-detection of a central image (after modelling selection effects) would favour a singular isothermal ellipsoid and limit deviations from it.
At low impact parameter the number of images increases by two: singular lenses (e.g.singular isothermal ellipsoid) form up to 4 images, while regular nonsymmetric lenses form up to 5, with the central image being typically the faintest.These offer enough information to reconstruct all parameters within the model, even in the GO limit.Moreover, as argued above, simply counting the images should lead to an upper bound on deviations from the SIS.Selection effects need to be modeled, accounting for the chance of images not being detected (e.g.due to demagnification or detector duty cycle).This possibility complicates the identification of the central image and requires a more sophisticated analysis.

VI. CONCLUSIONS
We investigated the phenomenology of GW lensing and its potential to probe cosmological structures utilizing WO effects.In Sec.II we reviewed the equations governing GW lensing and the methods used to obtain predictions (see Ref. [5] for details).
We then addressed the phenomenology of two symmetric lens models that generalise the singular isothermal sphere (SIS): the power-law generalised SIS (gSIS, Sec.III A), in which the slope k of the lensing potential is arbitrary and the cored isothermal sphere (CIS, Sec.III B), with a finite-density core of size x c .For each lens we presented the physical scales, GO structure (images, caustics, magnifications), bGO corrections, and WO predictions, including the convergence to the SIS limit.We mostly focused on the strong-lensing regime with mild magnification, although we showed examples of WO signatures in the single-image regime (weak lensing) and for sources near the caustics (Figs. 8, 9).
A distinctive feature of these lenses is the existence of a central image associated with the maximum of the Fermat potential.Central images produce a secondary modulation in the amplification factor, encoding additional information about the lens.For sources near a caustic, the additional modulation of the amplification becomes very broad in frequency space and pushes convergence to GO towards very high dimensionless frequency w ≫ 1/min |ϕ I − ϕ J |.The central image vanishes for SIS and narrow gSIS profiles (k ≤ 1): its detection is thus a direct probe of the inner structure of halos.We discussed qualitative aspects and observational prospects (Sec.III C), including parameter estimation in the w ≫ 1 limit, the prospect of detecting central images and mismatches between GWs lensed by gSIS and CIS, relative to SIS.
We then addressed the potential for GW detectors to probe lens features (Sec.IV).We used a Fisher-matrix analysis to forecast the ability of LISA and LIGO to reconstruct the parameters of the lens (mass, impact parameter, slope/core size) and of the source (signal amplitude and initial phase), as well as their correlations.We focused on equal-mass, non-spinning sources and kept their parameters fixed (except the amplitude and overall phase).While our analysis is simplified, it produces reasonable agreement with more detailed treatments of the source and the detector [36].The large degeneracy between lens parameters makes the Fisher matrix very sensitive to numerical errors in the computation.To circumvent this issue, we used approximations to the Fisher matrix at low/high frequencies and focused on typical strong-lensing situations (away from caustic) where convergence to GO occurs at relatively low dimensionless frequencies.We also addressed the validity of the linear signal approximation (LSA) in different cases.
LISA and LIGO can identify the CIS lens parameters with precision ∼ O(1/SNR) for large w chirp ∝ M Lz /M BBH (GO limit).The precision is somewhat higher for the core size and lower for the lens mass.The marginalised constraints saturate in the GO limit because of the large degeneracies between lens parameters, although the 2D posteriors become thinner as M Lz increases.This behaviour is generic in strongly-lensed systems with a central image (three images total), allowing to recover up to four parameters in GO (two time delays with high accuracy ∼ 1/M Lz and two magnification ratios ∼ 1/SNR).The constraints degrade at lower lens masses (such that the merger occurs below the onset of magnification), due to both parameter degeneracies and diffraction suppressing the amplification factor.
Wave optics can substantially improve the constraints in regimes where GO does not provide enough information.We saw this explicitly for the narrow gSIS lens, which only produces two GO images.In this case, the most stringent constraints are found when the onset of magnification occurs close to the merger, w chirp ∝ M Lz /M BBH ∼ 1 − 10.The constraints are driven by WO and degrade for increasing M Lz , as the bGO corrections become negligible.The study of the gSIS revealed that the LSA overestimates the sensitivity at high M Lz due to the contribution from the lens' cusp.For this reason, in addition to full WO we quoted bGO-only results, which neglect the cusp contribution and are thus conservative for w chirp ≫ 1.Our Fisher matrix results are supported by the mismatch analysis (Fig. 11), which does not rely on the LSA.Given the limitations we uncovered for the Fisher analysis on lensed signals, it would be interesting to perform more sophisticated forecasts by sampling the full likelihood and better asses past results obtained in the literature (e.g.[4]).
Identifying lens features at percent and sub-percent precision (for LIGO and LISA respectively) opens up novel applications of GW lensing.While GW lensing measures projected quantities, limits on the lens mass and physical scales can be obtained under reasonable assumptions regarding the cosmology and the halo mass function (Sec.V A and Fig. 21).In particular M vir , (ξ 0 ) has a hard lower (upper) limit, and their 95% credible regions can be obtained within a factor of a few, except for very high redshifts.These limitations may be lifted if the lens redshift can be identified, e.g.via EM follow-up [112,117].
The possibility of simultaneously constraining the lens mass and core size enables tests of Dark Matter scenarios.We addressed the prospects of probing selfinteracting and ultra-light DM models (Sec.V B and V C, respectively), which predict cored halos whose size depends on the properties of DM and the halo mass.Deriving simple expressions for the projected core size in these models, we translated the sensitivity from our forecasts into the parameter space accessible to strongly lensed GWs.GWs can probe SIDM cross sections as low as 10 −4 − 10 −5 cm 2 /g, with the maximum sensitivity achieved around the onset of magnification (Fig. 22).For ULDM, the core sizes of DM halos can be tested if m ϕ ≲ 10 −17 eV (cf.Fig. 23).These bounds assume the largest possible projected core size for a given source and are thus conservative within our model.Nonetheless, the assumption of a ρ ∝ r −2 profile might give overly optimistic constraints for SIDM, compared to a shallower central profile such as NFW.
The capacity to detect very small cores (recall x c is normalized by the Einstein radius) leads to high sensitivity, particularly for lighter halos, and in some cases well beyond current limits.The capacity of GWs to constrain M Lz and x c independently may allow observations to discern DM-dominated halos (M vir ≲ 10 9 M ⊙ ), and thus whether the core is likely to represent a signature of DM properties, rather than baryonic effects.We note that turning these prospects into actual constraints requires further modelling of lens features, astrophysical effects and lensing cross sections, as well as robust detection of lensed GW signals.To illustrate the need for more realistic modelling, we discuss how our analyses, based on axisymmetric lenses, can be generalized to elliptical lenses (Sec.V D).Despite the challenges ahead, our results warrant further investigations into the potential of DM tests with strongly-lensed GWs.
Our analysis highlights several ways in which GW lensing is highly complementary to electromagnetic observations: 1.The low frequency of GW sources enables observation of WO effects in lenses with mass M Lz ≳ 10 2 , 10 5 M ⊙ (point-lens) and M vir ≳ 10 4 , 10 7 M ⊙ (extended) for LVK and LISA, respectively.WO provides information about a region of the lens plane of size ∼ w −1/2 .
2. WO effects provide additional information through the onset of magnification and mild frequencydependence at low/high dimensionless frequencies (beyond GO), as well as the contributions from cusps.These effects persist in the single-image (weak-lensing) regime at moderate w (Fig. 9).6.While lensed waveforms are sensitive to projected parameters (M Lz , y, x c ), limits on the physical scales and virial mass can be placed under reasonable assumptions, Fig. 21.In particular, M vir can be obtained within a factor of a few.

GW
7. Central images may enable novel tests of astrophysics and fundamental physics, such as probing DM through lens cores (Figs.22,23).Variations of these tests may rely on measuring the slope of the inner halo cusp or the effect of microlenses and supermassive black holes close to the central image.
Based on these prospects, we envision further extensions of our work regarding lens models, source properties, statistical analysis and towards other lensing regimes.
While we focused on simple symmetric lenses, we expect many of our conclusions to hold in more realistic situations.For instance, the central GO image probes the inner region of the lens, which may differ from the outer, averaged lens distribution.While this dependence is exact only for GO and symmetric lenses (the lens equation depends on the integrated projected mass), it should approximate more general settings.Further characterization may enable analyses based on central images complementary to those proposed for EM sources, e.g.[16,17,[20][21][22].Still, more complex lenses may mimic or obscure some of the signatures of cores and cusps.Studying other lenses [105], adding large-scale lens features (variable slope, ellipticity, external convergence & shear) and small-scale substructure (DM subhalos, stellar population) is thus relevant.Future work needs to address these modelling challenges and opportunities.
Another interesting direction is a more detailed exploration of DM scenarios.By modelling DM halos as CIS, our work has shown the promise of GW lensing under rather simplistic assumptions.As discussed in Sec.V, it will be necessary to model lenses more accurately, along the lines mentioned above.Our work can, in principle, be extended to other DM models.Warm dark matter halos have been shown to have small cores, r c /r 200 ≲ 10 −3 [136].Another promising signature is the inner halo slope ρ ∝ r −3/2 [137,138].This feature is also present in other models such as compact DM structures [63], with a sharp slope ρ ∝ r −9/4 (although see Ref. [139]) and constant density core, whose size is set by the DM velocity dispersion or annihilation cross-section [140].Ultimately, deriving constraints will require a detailed assessment of lensing probabilities, including false positives due to baryonic or complex structures.
A key advantage of GWs is the existence of accurate source models.While we considered the source parameters as fixed (except for distance and phase) and simple sources (equal mass, non-spinning binaries), further understanding the interplay between source models and lensing effects is necessary.Besides varying the source parameters, other directions include the study of different sources (e.g.Extreme-mass-ratio inspirals), and more realistic waveform models (higher harmonics, spinprecession, eccentricity, see Ref. [36]).Ultimately, full Bayesian parameter estimation of lensed signal injections will be necessary to derive accurate posteriors and address the shortcomings of the Fisher matrix.These studies will allow an understanding of both parameter degeneracies and potential systematics, e.g. if a lensed event is analyzed with an unlensed waveform.Eventually, WO amplification factors need to be included in GW searches and parameter estimations, perhaps taking advantage of machine learning algorithms, e.g.[141].
Our results also warrant investigating other regimes of GW lensing.Here we have focused on the strong-lensing regime, at mild magnifications.An interesting extension is highly magnified configurations in which the source is close to a caustic [39,142], and the maximum magnification is limited by wavelength rather than source size [143,144].This requires pushing our numerical methods to new limits, especially in computing forecasts.Another future direction is the study of line-of-sight halos [145][146][147], and explore the interplay between multi-plane lensing and wave optics [148] to enable ray-tracing in cosmological simulations.Finally, wave optics in the singleimage regime is also promising, as the frequency dependence may be detectable up to large impact parameters [4,6,33,36,149].Since the lensing probabilities scale as the square of the impact parameter, weak-lensing effects might be the first detectable signature, or even contribute as a systematic in the waveform reconstruction.
As the detection rate increases, GW lensing is bound to become a reality.It is a matter of chance whether LVK searches return confident identification of lensed signals in the coming years.However, the capacity of next-generation detectors to observe binary mergers in the entire Universe [11] will make GW lensing into a reality: even the least likely lens configurations will stand a chance to be observed.Other proposed missions will offer new opportunities, such as WO effects for galactic scale lenses from µHz space-borne detectors [103].As GW astronomy matures, wave optics effects may serve both to identify lensed sources and to study the matter distribution in the Universe.The complementarity between lensing of GW and of EM sources may enable novel probes of the structure of halos, perhaps helping us identify the nature of DM and opening new insights into astrophysics, fundamental physics, and the history and evolution of the Universe.
where q ≡ e iπ/4 z/ √ w.This expansion is well defined provided that ψ(q) grows less rapidly than the quadratic part of the Fermat potential as z → ∞.
For the lenses we are going to discuss (see Tab. I), the integrals in Eq. (A1) can be evaluated explicitly.By suppressing the overall factor e −iwϕm , for the point lens ψ(x) = log x we have where γ E is the Euler's constant.Similarly, for SIS While for gSIS we have Finally, for the CIS: We will use these expressions to compute F (w) and its derivatives at low frequencies.

Wave Optics: Numerical Methods
Numerical methods are necessary to compute the amplification factor at intermediate frequencies, bridging the gap between the low-and high-frequency expansions described above.In Ref. [5] we developed and validated two numerical methods, regularized contour flow [3] and complex deformation [150].Both yield the same results at high accuracy, but the regularized contour flow is substantially faster and will be our method of choice.We will give an overview now; further details can be found in Ref. [5].
The starting point is the Fourier transform of the integral in (14) where δ D (x) is the Dirac-delta function.The above integral is computed on contours of equal Fermat potential ϕ(x, y) = τ .As contours end on critical points, the first step is to solve the lens equation to find the initial/final conditions for the contours.The integral is then sampled on the nodes that define each contour.Then the nodes are evolved to the next value of τ , adding or removing more nodes depending on the curvature of the contour.The next value of τ is chosen adaptively and depends on the derivatives of Ĩ(τ ).
The amplification factor (14) in the frequency domain is then obtained via a fast-Fourier transform (FFT).This is computed by combining all contours, interpolated on a homogeneous grid, as defined by the frequency range of the FFT.In order to improve the computation, we split Ĩ(τ ) into a regular and a singular part.The singular part is the time-domain counterpart to the GO amplification factor (21).It can be computed analytically and consists of discontinuities (for maxima/minima of ϕ) and logarithmic divergences (for saddle points of ϕ), which if not dealt with separately contribute spurious high-frequency noise to the FFT.
This algorithm is efficient and accurate at intermediate frequencies.At high frequency, we use the bGO approximation described in Eq. (24).The transition between the two regimes, w high , is chosen such that the relative deviation |F bGO /F WO − 1| is sub-percent.At the same time, w high is chosen relatively low to avoid numerical aliasing from the FFT, usually w high ∼ O(100).The choice for w high is more delicate when computing derivatives of F (w) with respect to the lens parameters since numerical noise is higher in these cases: Appendix B discusses this issue in the context of Fisher-matrix forecasts.

Singular Isothermal Sphere
The SIS is characterised by the following density distribution where σ v is the line of sight velocity dispersion.This lens is often used to model dark matter halos.Because of the falloff as ∼ 1/r 2 of the density profile, the halo mass is formally infinite, and the model thus requires a physical cutoff for r, after which the density is imagined to drop to zero quickly.In the context of lensing, this is not an issue as long as the cutoff is much larger than the scale ξ 0 specified below.The projected mass density of this lens is (Eq.( 5)) This simple lens model falls within the class of axisymmetric lenses since the projected mass only depends on |ξ|.The lensing potential is then obtained from Eq. ( 11) as ψ(x) = σ 2 v x/(ξ 0 G Σ cr ).In our work we use a standard choice for the normalisation of the lens, ξ 0 = σ 2 v /(G Σ cr ).With this choice the lensing potential becomes particularly simple ψ(x) = x . (A9) In the GO limit we can have one or two images depending on whether the impact parameter is outside or inside the caustic y cr = 1 (see Fig. 3, Panel C).For y < y cr , the two images have magnifications µ ± = 1/y±1, time delays ϕ ± = ∓y − 1/2 and Morse phases n + = 0 (minimum), n − = 1/2 (saddle).For y > y cr only the minimum image survives.Consistently with the prescription of Eq. ( 4), the minimum time delay needs to be set to zero by adding the appropriate ϕ m (y) = −y − 1/2.
In order to properly interpret our results, we now discuss the relationship between the quantities above and the physical properties of the SIS.To define the total mass of the SIS it is necessary to truncate the mass density at a finite radius.We will follow the standard definition of considering the virial radius r vir as the radius at which the density reaches 200ρ c , with ρ c = 3H 2 0 /(8πG) being the critical density of the Universe today (H 0 = 100 h km s −1 Mpc −1 is the Hubble constant), and choose this as the truncation radius.Then, the virial mass M vir is defined as the mass within the virial radius and coincides with the mass of the halo.Using Eq. (A7) and the expression for r vir , M vir is obtained as Alternatively, the expression above can be seen as a relation between the velocity dispersion and the halo mass.
With our normalisations, the effective lens mass is instead given by M Lz = 4π 2 (1+z L ) 2 σ 4 v d eff /G (see Eq. ( 16) and our choice of ξ 0 ).Using Eq. (A10) we can then relate M Lz with the virial mass: This relation, together with Eq. ( 19), suggests that for LISA frequencies and virial masses around M vir ∼ 10 9 M ⊙ , lensing is best described by wave optics.

Appendix B: Forecast implementation and tests
In this appendix, we are going to discuss the methods and validation procedures we employed in evaluating Fisher-matrix elements.The numerical derivatives of F (w) with respect to the lens parameters θ l are particularly delicate.Indeed, as we are going to mention soon, the large degeneracies in the Fisher matrix can amplify numerical errors when computing its inverse.It is then also important to develop approximation methods at low and high frequencies in order to validate our results in these regions.
At high frequencies we compute derivatives of the amplification factor using the beyond GO expression, Eq. ( 24), which can be arranged as a sum over GO images Here uppercase Latin indices I, J, ... refer to GO images, lowercase indices l, m, ... to lens parameters and commas denote partial derivatives with respect to the lens parameters (f ,l ≡ ∂f /∂θ l ).Additionally, F I (w) is the GO amplification of the I-th image, Eq. ( 21).Note that the Morse phases n I are assumed to be constant in the derivatives.
Parametrizing the effective lens mass by its logarithm provides a convenient rescaling of the Fisher matrix.Since M Lz only enters through the dimensionless frequency w, Eq. ( 15), the associated derivative reads (B2) Here the first equality is fully general and will be used in the intermediate and low-frequency regimes.The second expression is used in the bGO regime.Similar expressions are used for the derivatives of the cusp terms appearing in the gSIS lens, Eq. (26).
At low frequencies we compute F (w) and its derivatives using the low-w expansion of Sec.A 1 (one needs to re-introduce the minimum time delay ϕ m in these expressions, as it depends on the lensing parameters and so it affects the derivatives).Moreover, we checked these results do not change when using the full WO result at low frequency (computed using Eq. ( 16) of Ref. [5]).
Derivatives of F (w) in the intermediate region are computed in the time domain via finite differences and then Fourier transformed.This process significantly worsens the numerical noise at high frequencies, in a manner analogous to the introduction of terms ∝ w in the bGO derivatives, Eq. (B1).Given the sensitivity of the inverse Fisher matrix to numerical errors, it is necessary to achieve sufficient accuracy.To this end, we use a supplementary regularization procedure, described in App.A of Ref. [5].In addition, we choose w high , the transition between WO and bGO, such that the errors between both computations are smaller than 1% for F (w) and 3% for the derivatives with respect to the lens parameters.Figure 24 shows how these tests are implemented for the CIS and gSIS.We have found that achieving our target precision is hardest for the core size and the density slope.
The large degeneracy between lens parameters makes the inversion of the Fisher matrix (i.e. the posterior) very sensitive to numerical errors.This can be seen in the condition numbers, which are typically O(10 4 ) at best and much larger at very low frequencies w ≪ 1, as the lack of magnification prevents the recovery of lens parameters (see Fig. 17 in the main text and Fig. 27 for the dependence of the eigenvalues with the lens mass in a specific example).The conditioning number also grows as M 2 Lz when w ≫ 1 due to the degeneracies in lens parameters in the geometric optics limit.These degeneracies are responsible for the saturation of the precision at high lens masses (for fixed SNR), e.g.Figs. 13, 14 and 15.As a further test of the stability of the computation, we have checked the convergence between our full calculation, GODA and bGODA (introduced in App.C) at sufficiently high lens masses (see Fig. 25 in App.C).In the gSIS case with k > 1 (only two images) the Fisher matrix for GODA is not invertible: with only two images one can extract at most two lens parameters.The inclusion of bGO terms allows instead to extract more information, although the precision falls with M Lz , since GODA is recovered in this limit.Finally, the rGODA terms (containing also the cusp contribution) allow to formally obtain constant precision at high M Lz , cf.Fig. 26 in App. C. We explain this (paradoxical) result in Sec.IV C, motivating that it cannot be trusted, as the LSA becomes invalid in this limit.

Appendix C: Geometric Optics Diagonal Approximation
We will introduce an approximation valid at high frequencies, based on the bGO Fisher derivatives (B1).Here we will give expressions for GO, but bGO (and rGO) generalizations can be obtained straightforwardly.In the GO limit the Fisher-matrix elements involving lens parameters l, m read The integrand associated with each term is I Il,Jm = µ I,l µ J,m µ I µ J + 4w 2 ϕ I,l ϕ J,m cos(∆ IJ ) + 2w ϕ I,l µ J,m µ J − µ I,l µ I ϕ J,m sin(∆ IJ ) .(C2) The oscillatory pieces are given by the phase difference between images At high frequencies, the GO Fisher matrix will in general be dominated by the diagonal terms I = J, for which the oscillatory contribution cancels.We, therefore, introduce the Geometric Optics Diagonal Approximation (GODA) for the Fisher matrix, defined by neglecting the cross-image terms as follows: The first term is proportional to the unlensed signal's SNR 2 .The second term involves the signal weighted by the dimensionless frequency (wh 0 |wh 0 ) ∝ M 2 Lz (f h 0 |f h 0 ) and typically dominates in the high-mass limit.Note that GODA is equivalent to GO in the weak-lensing regime, as there is a single image.
Figure 25 shows the convergence of the full calculation to GO and bGO for an example lensed LISA source.Good convergence is achieved above lens masses such that the dimensionless frequency corresponding to the innermost stable circular orbit (ISCO) is large, w ISCO > 10.For lower lens masses the GODA overestimates the sensitivity: it traces the shape of the error in the mass, but completely overestimates other lens parameters.This is because information about the mass is only contained in the phase/time delays, and hence lost in the limit w(ϕ I − ϕ J ) ≲ 1.Other lens parameters affect the GO magnification factors, and hence appear to be measurable in the GODA limit even when w ≪ 1, as GODA fails to account for the lack of magnification in the low-w regime.WO has an additional source of information on the lens mass, namely the onset of magnification around w ∼ 1.
the strong-lensing regime, not only to the CIS.To our knowledge, this issue was not investigated previously in the literature.It would be interesting to revisit the Fisher results in the literature and characterize the sensitivities with a more sophisticated analysis.
As we are soon going to mention, this issue becomes more problematic in the gSIS case, where the Fisher analysis yields very optimistic results that contradict the number of parameters that can be reconstructed in the GO limit, Eq. ( 57).We give a concrete analytical example of this issue in App.E. As we comment at the end, the analysis of this example also supports the conclusion that a violation of the LSA condition for the CIS should not call for concern.

Appendix E: LSA approximation for cusps
In the main text we argued that the linearisation of the likelihood used in the Fisher matrix formalism might not be always appropriate, even for moderately high SNR.In turn, the linearisation procedure can overestimate the sensitivity on parameters, yielding unrealistic, or even paradoxical, results.Specifically, this seems to be the case for the parameter k in the gSIS forecast.For a fiducial value k = 1 the Fisher matrix predicts ∆k ∝ 1/SNR in the GO limit (large lens masses M Lz ).On the other hand, in this limit the characteristic feature due to k (the signal from the cusp) quickly decays.Thus, no meaningful constraint on k should be present.Here we provide a toy example highlighting this feature: by studying the full probability distribution we can see that indeed the error on k grows for large M Lz .
Let us consider a toy lensing signal with a single generic lensing parameter θ (that could represent k in the case of the gSIS).The lensed GW waveform h L (f ) is taken to be where h 0 (f ) is the unlensed waveform, w = 8πGM Lz and M Lz is taken as fixed.This signal is indeed analogous to the cusp obtained for the gSIS in the limit of k → 1, see Eq. ( 26).One can easily generalise what follows to cases with different w dependencies.Let us write θ = θ + ∆θ, where θ is the fiducial value.We want to compute the probability distribution for ∆θ, assuming some analytic behaviour for h 0 (f ) and the noise.
Let us call ∆h ≡ h( θ) − h(θ).Then, the exponent Λ of the likelihood is defined as where we are neglecting the noise term.In this way the likelihood is written as L ∝ e −Λ .We can write the full If we were to expand the square bracket for small ∆θ, we would recover the Fisher formalism.Notice that if we expand and trust the result for w∆θ > 1, then we are overestimating Λ.In the full result, we see indeed that the square bracket remains bounded for large ∆θ.
At this point, we can understand the behaviour of the integral by taking a simplified form for h 0 (f ) and S n (f ).The lower band of the LISA noise is characterised by S n (f ) ∼ s n f −4 , where s n is a dimensionful constant that depends on the noise properties of LISA.This is a good approximation for events with large masses M BBH ≳ 10 7 M ⊙ since the merger happens before the noise curve reaches its minimum.For the waveform, we take simply h 0 (f ) = Af −7/6 e iΨ(f ) , where A is the amplitude and Ψ(f ) the phase evolution [151].We perform the integration over f from some low f 0 , where the signal enters the LISA band, and we cut the signal at the ISCO frequency, which we indicate as f 1 .
The SNR for this signal is then where we assumed f 1 ≫ f 0 .We can use the expression above to simplify Λ, and get Λ ≃ 8 SNR 2 This integral can be solved explicitly in terms of exponential integrals.However, for our purposes, we can focus on the two limits of w∆θ ≪ 1 and w∆θ ≫ 1, where w represents the dimensionless frequency contributing mainly to the integral (i.e.w at the ISCO).
For small w∆θ we can expand the square bracket, recovering a quadratic function of ∆θ: Λ ≃ SNR 2 ∆θ 2 .This corresponds with the Fisher-matrix result.Naively, it would predict ∆θ ∼ 1/SNR, even though the correction due to θ in Eq. (E1) becomes negligible at high M Lz .From this discussion it is now clear that this result is only applicable for w ISCO ∆θ < 1 and therefore breaks down at high M Lz .This is analogous to the LSA condition in Eq. ( 66) for this simplified case.Instead, in the opposite limit, we can estimate Eq. (E5) by noting that the cosine function becomes rapidly oscillating.This means its contribution to the integral is negligible: the log-likelihood then becomes independent of ∆θ: This result shows θ cannot be efficiently constrained when w ISCO ∼ SNR, since its probability distribution stops being localized around ∆θ = 0 and becomes flat, with L ∝ e −Λ ∼ 1.In the more realistic situation where, on top of the cusp contribution, there are more GO images the situation is of course more involved.Indeed since the cusp is very faint, its contribution to the likelihood might be dominated by the product with other images.The general logic is expected to apply nonetheless.
In conclusion, we expect the Fisher approximation to stop being accurate around ∆θ ∼ w −1 ISCO , with the constraint on θ becoming weaker and weaker as we increase M Lz .Eventually θ becomes unconstrained when we reach w ISCO ∼ SNR.
We can also comment on the case in which the lensing effects are not suppressed by 1/w, as in Eq. (E1), but are due to a GO image.(Notice also that bGO terms do not have the form of Eq (E1), since the phase factor would be factorised outside of the bracket.)It is straightforward to obtain the analogous of Eq. (E5) in this case (the factor w −2 disappears).Linearisation becomes unreliable when ∆θ ∼ w −1 ISCO , as in the previous case.However, the likelihood is still highly suppressed for large w ISCO , since Λ ∼ SNR 2 (the tails of the distribution differ from a Gaussian, but the associated probability remains small).Applied to the CIS case, this suggests the Fisher analysis still reliably captures the size of the sensitivities even after the LSA condition is violated.

FIG. 1 .
FIG.1.Setup for strong gravitational lensing.A lens located at (angular diameter) distance DL magnifies a source at DS with (dimensionless) impact parameter y.In the geometricoptics limit, the lens produces multiple images (stars), whose positions are given by the dimensionless coordinates xI in the lens plane.GWs allow us to observe all images, including the central image (red), whose EM counterpart is further suppressed (1 > |µ| > |µ|), and blocked or outshined by the lens matter.

12 30 + 30 MM Lz = 3 • 10 4 MFIG. 2 .
FIG.2.Left: dimensionless amplification factor for an SIS with impact parameter y = 0.7.Different regions correspond to different lens masses for a given source, cf.Eq.(15).Right: time-domain waveforms for a typical LVK source (30 + 30 M⊙) for three different effective lens masses representative of the perturbative (A), "intermediate" (B) and geometric optics (C) regimes.Each panel corresponds to a shaded region in the F (w) plot.

FIG. 3 .
FIG.3.Overview of the lens models considered: SIS (dotted), gSIS with narrow (dashed) and broad (dash-dotted) profiles and CIS (solid).Panel A: radial density profile (arbitrary units).All except the CIS have a divergent density at the center.Panel B: lensing potential from the projected density.Narrow (broad) density profiles correspond to shallow (deep) lensing potential.Panel C: lens equation.Solutions are given by the set of points where the curves intersect y = const.Extrema of the lens equation separate regions with different numbers of solutions.Magnification is related to the slope of the curves.

FIG. 4 .
FIG. 4. Geometric optics image positions and magnificationsfor broad gSIS lens (Eq.(30) with k = 1/2).Top: Image positions associated with saddle point (blue), maximum (orange) and minimum (green) of the Fermat potential, with the projected density (shaded region).The saddle and maximum images exist only within the caustic (red dashed vertical line, Eq. (35)).The shaded area represents Σ as a function of x.Its height is rescaled to arbitrary units.Middle: Magnification for the GW amplitude for the different images from Eq.(36).The central image magnification goes to zero at low y and diverges near the caustic.Bottom: Beyond geometric optic correction times the magnification for different images, obtained from Eq. (37) and Eq.(36).The central image has a large contribution that overcomes demagnification at low y.

FIG. 7 .
FIG.7.Geometric-optics quantities for a CIS with core size xc = 0.15 as a function of the impact parameter.Top: Image positions associated with saddle point (blue), maximum (orange) and minimum (green) of the Fermat potential, with the projected density (shaded region) and core size (vertical grey line).The saddle and maximum images exist only within the caustic (red dashed vertical line, Eq. (53)).Middle: Magnification for the GW amplitude for the different images.The central image magnification varies between a finite minimum (∼ 0.429 for xc = 0.15) and the divergent value near the caustic.Bottom: Beyond geometric optics correction times the magnification for different images, obtained from Eq. (37) and Eq.(36).The central image has a large contribution that can overcome the geometric-optics magnification of the other two images for w ∼ 1.

8 FIG. 10 .
FIG. 10.Detectability criterion in GO: magnification of the third image as a function of the lens parameter.Magnifications above the threshold values SNR th = 8 are detectable by GW experiments.Horizontal lines correspond to typical LIGO and LISA events, with SNR of 30 and 10 3 respectively.Left: CIS case, as a function of the core size xc.Right: gSIS case, as a function of the slope parameter k.

2 MFIG. 11 .
FIG. 11.Mismatch between CIS and SIS (top row), and gSIS and SIS (bottom row) for impact parameter y = 0.3 (left column) and y = 0.6 (right column) for a source mass MBBH = 10 6 M⊙.The mismatch is evaluated for different lens masses MLz and by varying the lens parameter of the lens (xc for CIS and k for gSIS).

4 FIG. 12 .
FIG. 12. Summary of LISA sources.Left: Fourier-domain waveforms for different source masses.Dashed lines show the unlensed waveform, solid lines are lensed by a gSIS with MLz = 10 7 M⊙, y = 0.3, xc = 10 −2 .Frequency components emitted before tmerge < 1 Mo, neglected in our analysis, are shown in lighter color.Right: corresponding lensed SNR as a function of the effective lens mass, including only the last month of the signal before merger.

FIG. 13 .
FIG.13.Constraints on the lens parameters (log(MLz), y, xc) as a function of the lens mass for LISA sources (cf.Fig.12) lensed by a cored isothermal sphere with y = 0.3, xc = 10 −2 .Results are normalized to a lensed SNR of 10 3 .Fiducial values are shown as horizontal dashed lines.Vertical lines mark the breakdown of the LSA, condition(66).

FIG. 14 .
FIG.14.Constraints on parameters of the lens (log(MLz), y, xc) and source (log(DS), ϕ0) and lensed SNR as a function of the lens mass and the impact parameter for a cored isothermal sphere with xc = 10 −2 .The source is a MBBH = 10 6 M⊙ equal-mass non-spinning binary observed by LISA.Results are normalized to a lensed SNR of 10 3 .Vertical lines mark the breakdown of the LSA condition(66).

FIG. 15 .
FIG.15.Constraints on the lens parameters (log(MLz), y, xc) as a function of the lens mass and the core size xc for a MBBH = 10 6 M⊙ LISA source lensed by a CIS (y = 0.3).Results are normalized to a lensed SNR of 10 3 .Vertical lines mark the breakdown of the LSA condition(66).

3 M 10 FIG. 17 .
FIG. 17. Left: error associated with the principal components of the Fisher matrix for the CIS lens, ordered from best to worst constrained (Fisher matrix, its inverse and eigenvectors for MLz = 10 4 , 10 7 , 10 10 M⊙ are shown in Fig.27).Right: 68% and 95% c.l. MLz − xc marginalised posteriors for the largest masses.For clarity, we show one quadrant in symlog scale.

FIG. 19 .
FIG.19.Constraints on the gSIS lens (log(MLz), y, k) and source (log(DS), ϕ0) parameters, and lensed SNR for different impact parameter (fiducial slope k = 1+10 −4 ).The source is a MBBH = 10 6 M⊙ equal-mass non-spinning LISA binary.Results are normalized to a lensed SNR of 10 3 .For values of MLz in the grey-shaded area, the LSA criterion Eq. (65) is violated.This criterion is obtained as ∆LSA < 1 and also by requiring we are in the GO limit (wmax > 100).Dashed lines are the results of the Fisher forecast using the bGODA approximation which neglects the cusp contribution.The correct constraint should lie between WO (which overestimates the contribution from the cusp) and bGODA (which neglects it completely).

FIG. 20 .
FIG.20.Constraints on the gSIS lens parameters (log(MLz), y, k) for different fiducial lens slope k values (y = 0.3).The source is a MBBH = 10 6 M⊙ LISA binary.Results are normalized to a lensed SNR of 10 3 .The grey area is the region where the LSA condition is violated and the dashed lines are the Fisher forecast results in the bGODA approximation (see caption of Fig.19).The solid dark-green lines (k = 1.01) are expected to start growing at large lens masses, following the qualitative trend of the light-green curves (k = 1.3).This ultimately follows from the scaling of the cusp contribution Fc(w) as a function of w (see Eq. (26)).

−1/ 3
vir for the SIS, cf.Eq. (A11).The resulting probability distribution is shown on panel C of Fig. 21.Higher values of M Lz are more skewed towards lower z L because heavy halos are rare at high redshift.

FIG. 22 .
FIG.22.GW lensing as a probe of SIDM for the CIS profile.The shaded regions can be probed by strongly lensed GWs.Here y = 0.6, zS = 3, MBBH = 10 6 M⊙, SNR = 10 3 for LISA and zS = 0.2, MBBH = 30 M⊙, SNR = 10 2 for LIGO.The lens redshift is chosen to minimize the projected core size (see text).Constraints on the cross section scale as ∝ ρ(rc), Eq. (71).Therefore, results for an NFW profile are expected to scale as the square root of the CIS ones, compare Eq. (72) and(73).The forecasted CIS results are orders of magnitude tighter than existing constraints[70].Ultimately, probes based on very small cores (low ⟨σv⟩/m) will be limited by astrophysical uncertainties.

3 FFIG. 24 .
FIG.24.Convergence of the amplification factor and the Fisher derivatives to bGO and rGO respectively for a CIS lens (left) and gSIS (right).The bGO/rGO values are used for high frequencies (shaded region) to avoid aliasing noise.
, xc is approximately equal to xc = 6.4 • 10 2 10 12 M ⊙ lensing is a unique probe of the inner structure of halos.Central, demagnified images are brighter for GWs than for EM sources ( |µ H | ≫ |µ H |),