Weakly Lensed Gravitational Waves: Probing Cosmic Structures with Wave-Optics Features

,

Gravitational lensing, the effect of gravitational fields on the propagation of signals through the universe, predicts a plethora of observable effects [1,2].Many gravitational lensing phenomena have been observed using light and other electromagnetic (EM) signals, leading to a wide-range of applications in astrophysics, cosmology and fundamental physics.Now, the advent of gravitational-wave (GW) astronomy provides the prospect of observing novel lensing phenomena.
Lensed GW signals stand out as highly complementary to EM observations.GWs can be detected at very high redshift and are free from many of the systematic uncertainties present in EM probes.For this reason, GWs have been proposed as an alternative tool to test the cosmological model, e.g.studying the cross-correlation of GW observations and galaxy surveys [3][4][5].Moreover, the low frequency and phase-coherence of observable GW sources make them ideal ground to probe the wave-optics (WO) regime [6].WO encompasses phenomena such as diffraction and interference, which imprint frequencydependent signatures on GW waveforms.They can hence be used to identify a signal as lensed [7][8][9][10][11][12] and even infer the lens properties accurately [13][14][15][16], at least in the strong-lensing regime.In contrast, identifying lensing in the geometric-optics (GO) limit, i.e. the high-frequency limit, requires associating multiple images from the same event, a method prone to false alarm [17], or identifying subtle waveform differences [18], which requires sources with large mass ratios [19,20].
WO signatures require lenses in a restricted mass range, set by the frequency range of the observed signal.For observable GWs, this limitation implies that only relatively light structures can be detected.This dramatically reduces the probability of detecting WO features, at least for strongly-lensed signals in which the most likely lenses are massive galaxies [21].This led to a pessimistic prospect to detect WO imprints in strongly lensed signals, e.g. by LISA [22].
WO features can also be searched for in weakly-lensed signals, which do not require a close alignment of sourcelens and observer and thus have a higher probability of occurring.It was estimated that LISA may detect WO effects at ∼ 50× the strong-lensing impact parameter, corresponding to O(1%) of massive BH binaries [23].While the above study was based purely on a mismatch analysis between lensed and unlensed waveforms, a more detailed estimate (accounting for waveform and lens parameter correlations) yields comparable values [14].The trade-off between strong and weak lensing is that of rare and dramatic versus frequent but subtle signatures.
The lens distribution can be probed in a rather transparent way via weakly-lensed GWs that contain WO features.Reference [24] showed that the frequencydependent amplification factor is determined by the shear of the Fermat potential at distances from the source given by the Fresnel radius ∝ 1/ √ f .In this way, weakly-lensed GWs probe entire regions of a lens as the source inspirals towards merger.In contrast, signals in the GO regime are sensitive to a very small portion of the lens plane where the image forms.Thus, WO provides a unique opportunity to test the structure of gravitational lenses.
The purpose of this paper is two-fold: first, we address the problem of computing WO signatures.Focusing on the single-image regime, we develop a framework to compute lensed waveforms of arbitrary lenses, efficiently and accurately.Second, we use these tools to explore the phenomenology of weakly lensed gravitational waves, the possibility of inferring structural features of lenses and the prospect of detection.The paper is organized as follows.Section II summarizes the WO regime of gravitational lensing and introduces two frameworks for the single-image regime: a general numerical computation and an expansion suitable for the weak-lensing (WL) limit.In Section III we discuss WO phenomenology, using a Green's function approach to analyze symmetric lenses, before addressing the case of a lens with substruc-ture.Section IV discusses the prospects of observation by future detectors.Possible future applications of our formalism are presented in Sec.V. We conclude by discussing our results in Sec.VI.

II. GRAVITATIONAL LENSING IN THE SINGLE-IMAGE REGIME
In this section, we develop a framework to compute diffraction effects in the single-image regime.We will first summarize the WO formalism in the frequency and time domain (Sec.II A).We will then present a method to numerically evaluate single-image signals in the time domain (Sec.II B).Finally, we will develop a perturbative WL expansion in the time domain (Sec.II C).

A. Wave optics formalism
In the frequency domain, the effect of lensing is characterized by a multiplicative factor F (f ), called amplification factor : where h0 (f ) and h(f ) are respectively the Fourier transforms of the unlensed and lensed strain amplitude.The frequency-domain amplification factor is obtained as (see Ref. [1] for a derivation assuming the weak-field limit, the thin-lens approximation and a static configuration).The integration is over the lens plane, with the coordinates rescaled by an arbitrary dimensionful scale ξ 0 (e.g. a characteristic scale of the lens), so x is dimensionless.The impact parameter y is rescaled by η 0 ≡ D S ξ 0 /D L , where D S , D L are the angular diameter distances to the lens and the source, respectively.
Here we introduced the dimensionless frequency which is given in terms of a redshifted effective lens mass: The factor d eff ≡ D L D LS (1+z L )D S also depends on the angular diameter distance between the observer and the source D LS .For a point lens, M Lz is equal to the total mass of the lens times (1 + z L ) if ξ 0 is set to the Einstein radius.However, this is not true for extended lenses (e.g.Eq. ( 46) below).
The integral in Eq. ( 2) depends on the Fermat potential : which is a dimensionless version of the time delay.Here ψ(x) is the lensing potential, which depends on the matter distribution projected on the lens plane and whose derivative gives the deflection angle.In particular, it is obtained as the solution of ∇ 2 x ψ(x) = 2Σ(ξ 0 x)/Σ cr , with ∇ 2 x being the 2D Laplacian, Σ(ξ 0 x) the projected matter density of the lens, and Σ cr ≡ (4πG(1 + z L )d eff ) −1 the critical density.We shift the Fermat potential by a constant ϕ m (y), defined in such a way as to make the minimum time delay equal to zero.
An important case of WO lensing is the GO limit, corresponding to the w → ∞ limit of the diffraction integral Eq. (2): Here the index J labels the GO images, located at stationary points of the Fermat potential x J such that ϕ ,i (x J , y) = 0, where a comma subscript indicates derivative with respect to lens-plane coordinates.The magnification µ −1 ≡ det (ϕ ,ij (x J )) and the time delay ϕ J ≡ ϕ(x J , y) are evaluated on the image positions.The Morse phase is n J = 0, π/2 or π depending on whether x J corresponds to a minimum, saddle point or maximum of ϕ, respectively.In the single-image regime, GO is simply a rescaling of the waveform, F (w) = |µ|.
We will now compute the amplification factor in time domain.We define the time-domain signal as the Fourier transform of iF (w)/w: where δ(x) is the Dirac-delta function.The expression above reduces the computation of I(τ ) to a onedimensional integral over contours of equal time delay ϕ(x, y) = τ , see [25].The amplification factor (2) follows from Fourier-transforming back to the frequency domain.
By choosing coordinates that follow the contours, the equation above reduces to where the coordinate t ≡ ϕ(x, y) and s is the arc-length distance along contours of equal time delay.The summation is over distinct contours with same time delay.We give a detailed derivation of this expression in App. A. We will consider the Green's function, defined as (also the Fourier transform of the amplification factor).
The time-domain lensed waveform is given as a convolution of the unlensed waveform h 0 with Green's function where t ≡ 4GM Lz τ .The GO image in G(τ ) appears as a singular contribution, stemming form the discontinuity of I at τ I (a Dirac delta function in the single-image regime).We will hereafter use the term Green's function when referring to the regular part, defined as When necessary, we will refer to G(τ ) as the full Green's function.
Figure 1 shows the procedure to compute WO predictions in the single-image regime.The lens is a Singular Isothermal Sphere (SIS), see Appendix B. Panels show the countours of constant Fermat potential for an example lens, the time-domain integral (8), Green's function (11) and the amplification factor (2). Colored points in I(τ ) and G(τ ) correspond to the contours in ϕ(x).Most discussions on WO lensing have focused on the amplification factor F (w).However, Green's function offers complementary insights into WO phenomena and their relation to lens properties, making certain features particularly transparent.Unlike F (w), G(τ ) is real-valued and thus easier to display.

B. Non-perturbative single-image framework
In this work, we focus on the single-image regime, in which only one GO image forms.The WL limit emerges as a particular case, in which deflections are small.While the WL regime is amenable to a perturbative treatment, which we will develop in Sec.II C, in this section we will first present the full framework needed to compute the amplification factor, in the single-image regime, without additional approximations.
The starting point for the full computation of the timedomain amplification factor will be Eq.( 7), but expressed in polar coordinates where x m is the location of the minimum time delay, i.e. ϕ(x m ) = 0.With this change of coordinates we get Our main assumption to solve this integral will be that we are in the single-image regime, so that the global minimum is the only critical point of the Fermat potential.Furthermore, if ∂ R ϕ ̸ = 0, we can invert  13).Green's function G(τ ) is then computed as the derivative of I(τ ), see Eqs. ( 9) and (11).Finally, the amplification factor F (w) is the (inverse) Fourier transform of the full Green's function G(τ ).The sharp peak in Green's function is associated with the center of the lens, which features a cusp in the SIS.In our case, the lens is located at (0, 0) and the corresponding contour is represented with a dotted line.
to obtain R(θ, τ ).Once this solution is found, we can plug it back and compute the integral as In practice, what we will do is to solve the system of differential equations In this way, the curve will be sampled with the precision needed to achieve a given tolerance in I.The system is integrated from θ = 0 to 2π and with initial conditions I(θ = 0, τ ) = 0 and R(θ = 0, τ ) chosen such that ϕ(R(0, τ ), 0) = τ .The previous derivation relied on the fact that ∂ R ϕ ̸ = 0, which is always the case for axisymmetric lenses when there are no critical points.Even though we will not need it in this work, the previous framework only needs to be slightly modified if this is not the case and ∂ R ϕ = 0.The main change to be made is that, instead of parameterizing the curve (14) as R(θ, τ ), one must use a parametric representation R(σ, τ ) and θ(σ, τ ).In this case, one should also keep track of the values of R and θ and finish the integration once the contour closes.

C. Perturbative weak lensing expansion
We can understand how to set up a perturbative calculation in WL in the following way.Let us consider the image to be at x m in the lens plane, as in Fig. 2, and let us assume that ψ(x) grows less than the quadratic part of ϕ(x, y) at large x.Then, at sufficiently large xs from the image, the contours of constant ϕ(x, y) = τ approach circles centered at y and are weakly influenced by ψ(x).One can then take into account the effect of ψ(x) perturbatively.On the other hand, at radii comparable with the distance |y − x m | or smaller, the contours are still weakly affected by ψ(x), but cannot be parametrized at the lowest order as circles centered in y.Indeed, the correct parametrization here is with ellipses centered at x m .This can be understood as the GO limit for I(τ ) since regions of small time delay correspond to the highfrequency limit for F (w) (see [25,26]).The two calculations for small and intermediate/large time delays can then be matched in an intermediate region.

Large time delays
The first step of the perturbative approach is to split the Fermat potential into a lens contribution and a "free" part, ϕ(x) = ϕ 0 (x) − ψ(x), where the free piece still contains non-perturbative information about the lens in the minimum time delay ϕ m .After plugging this result into (7), we then expand in powers of the lensing potential ψ In the second line, δ (n) stands for the derivative of the Dirac delta with respect to its argument.Without loss of generality, the impact parameter y can be taken to be parallel to the x 1 axis and with magnitude y.We can use polar coordinates again, centered at the minimum of the time delay of the free case (y, 0), to evaluate the integrals in (19).In general, for a generic function f (x), we can simplify the integral as where τ ′ ≡ τ + ϕ m and Θ is the Heaviside step function.Using this result, we can finally write the linear approximation as where This formula already captures all the essential diffraction features of the amplification factor with a very good accuracy that improves as y increases.In the next subsection, we will study the region of small time delays, −ϕ m > τ ≥ 0, where the linear formalism cannot be applied anymore, but analytic results from the GO expansion are available.We now discuss the frequency-domain version of our approximation for WL signals.This will serve to illustrate how the WOF appears in the amplification factor.However, for applications we will evaluate lensed signals starting from the time domain.Also, in the frequency domain, WL effects are given by an expansion in powers of ψ(x).We will show that at leading order in the lensing potential and for large y, the amplification factor is obtained by Fourier-transforming the signal from the large-time delay region.To obtain this result, let us write F (w) ≃ F (0) (w) + F (1) (w), where F (0) (w) and F (1) (w) are the Fourier transforms of I (0) (τ ) and I (1) (τ ), respectively.Then, by using the expressions in Eq. (22) and by performing the dτ integration using the delta function, we have Notice that ϕ m and ϕ 0 (x, y) still depend on the lensing potential.However, we are interested in keeping only leading order terms in ψ.Hence, we can expand the exponents of Eqs.(24) and (25) in powers of ψ m ≡ ψ(x m ), truncating at linear order.We also make use of the lensing equation at leading order: x m ≃ y + ∇ y ψ(y).This gives ψ m ≃ ψ(y).Then, expanding F (0) (w) up to first order in ψ(y) and adding the contribution from F (1) (w) leads to One can check that this expression correctly captures the WL features.Moreover, in the limit of large y the GO result is approximately recovered.Indeed, in this limit the location of the image approaches x m ≃ y.Expanding the integrand in Eq. ( 26) around this point and performing the Gaussian integral we obtain, at leading order in w ≫ 1, F (w) ≃ 1 + ∇ 2 x ψ(y)/2 ≃ |µ|, as expected (see next subsection and the GO expressions in App.A).Equation ( 26) can be applied to simple lenses to obtain analytic expressions in the WL regime.In App.D we present the result for the SIS lens.An expression similar to Eq. ( 26) is also derived in Ref. [24]. 1 1 From our understanding, this reference subtracts ψ(xm) instead of ψ(y) in Eq. (26).In this way, their F (w) grows at high w and

Small time delays
For τ approaching the minimum time delay, the formalism above cannot be straightforwardly applied.The main effect of the lensing potential in this region is to shift the minimum (from the lens equation, at leading order in ψ one has x m ≃ y + ψ(y)).At subleading order we also have a deformation of the contours contributing to a change in the magnification of the image.
The expansion of the contours near the minimum time delay leads to the GO expansion, which corresponds with the high-frequency limit of F (w).One can systematically obtain this expansion without making assumptions about the size of ψ.Following [15,25], we have in the time domain Here, µ is the magnification factor of the image, while ∆ (1) and ∆ (2) are the first two beyond geometric optics (bGO) corrections.We give the explicit expression for these coefficients in App. A. Higher-order terms in τ can be obtained in a similar fashion.
In the frequency domain, Eq. ( 27) becomes Notice that all these expressions only require knowledge of the image location x m in order to be evaluated.Moreover, higher orders in the bGO expansion decay with higher powers of w and are therefore subleading at high frequencies.On top of the bGO terms associated with the image, other locations in the lens plane can contribute at subleading orders in 1/w.For instance, this is the case for points where ψ(x) develops cusps (typically at the lens' center).We will elaborate on this point later, see also [26,27] for more details.
In order to connect with the expansion of Sec.II C 1 we need to pick a time delay τ match where to match the two expressions.In practice, in the large-τ expansion, it is convenient to use τ ′ = 0 (or τ = −ϕ m ) as the matching point and it is usually enough to keep only the leading order term in (27).As shown in App.E, it is possible to achieve O(1%) accuracy for y > 2 by setting ∆ (1) = ∆ (2) = 0 and interpolating between τ = 0 and τ = −ϕ m computed using (22) with n = 1.

III. ANALYSIS OF LENSING DIFFRACTIVE FEATURES
The framework developed in the previous section allows us to compute wave-optics features (WOFs) in the does not reduce to the GO limit.Therefore such an expression does not reproduce the WOF features we discuss in the next sections.
single-image regime.Let us now discuss the phenomenology of WOFs, their dependence on the lens properties and the prospect of individually identifying and characterizing sub-lenses.We will first start with the analysis of symmetric lenses (Sec.III A) before addressing models of composite ones (Sec.III B) and their signatures (III C).

A. Symmetric lenses
Let us start our discussion by considering simple, symmetric lenses (where ψ only depends on x ≡ |x|).First, we will introduce the symmetric lens models (a detailed description of these lenses and their phenomenology is given in Ref. [15]).Then, we will discuss the WOF in the time domain and their dependence on the lens parameters.Since GW analyses are often performed in the frequency domain, we will also discuss the WOF as a function of w.
Our discussion of symmetric lenses will focus on the well-known SIS and a one-parameter generalization, the Cored Isothermal Sphere (CIS).The SIS is characterized by a central cusp with diverging density, Σ(x) ∝ x −1 .In the GO limit, the SIS can have one or two images depending on whether the impact parameter is outside or inside the caustic y cr = 1, respectively.The CIS has finite density with Σ(x) ∝ x 2 + x 2 c −1/2 , where x c ≡ r c /ξ 0 is the projected size of the core.Similarly to the SIS, multiple images form for sources within the caustic y rc (x c ) ≤ 1 (smaller than for SIS, multiple images also require x c < 1/2).More details about both lens models are given in App.B.
Let us now describe the WO phenomenology of these lenses in the single-image regime.Figures 3 and 4 show, respectively, the predictions for an SIS at varying y and for a CIS at fixed y but with different x c s.For concreteness, we will first discuss Green's function G(τ ), which makes the discussion especially transparent, and comment on the amplification factor F (w) afterwards.We will describe the overall structlensing gravitational waves o3ure of the WOF, then the role of lens parameters and discuss how they could be measured from GW observations.
Single-image WOFs begin at τ I , which corresponds to the minimum of the time delay (i.e. the type I image, we set τ I = 0 by convention).The GO image appears as a delta function in the full Green's function, Eq. ( 11), while the WO piece G(τ ) features a discontinuity, associated to the bGO correction at the position of the image (originating from the coefficient ∆ (1) in Eq. ( 27)).G(τ ) is initially positive and increases with τ , as the contours approach the center of the lens.The location of the lens is associated to a peak in Green's function.At slightly higher τ , G(τ ) becomes negative and asymptotes towards zero as τ → ∞.
For symmetric lenses with a cusp (e.g.SIS) the peak of the WOF is located at τ C = ϕ(x = 0, y).The peak is due to the high curvature of the constant-ϕ contours.The curvature, and hence the height of the peak, depends on the lens profile as well as the impact parameter: denser lenses and lower impact parameter produce taller peaks.
In non-differentiable lenses the peak is singular.For instance, in the linear weak-lensing approximation for the SIS lens, Green's function is found to have a logarithmic divergence (see App. D for a derivation) In the frequency domain, the peak in G(τ ) is directly related to the damped oscillations seen in |F (w)|: sharper peaks have more pronounced features that decay more slowly with w, and can thus be observed at higher frequencies.The singular contribution is related, in the frequency domain at large w, with subleading terms in the bGO expansion originating from regions of the lens plane close to the center of the lend.See discussion in App.D around Eq. (D7) and Refs.[15,27] for more details.
Let us now discuss how the lens parameters affect the WOF, separating the peak and broad shape.WO predic-tions are independent of M Lz when expressed in terms of y, τ and w.However, M Lz can be inferred when restoring the units for t or f , given an observed waveform.We will discuss this at the end of this subsection, together with prospects for lens parameter recovery.
The impact parameter y controls the position and amplitude of the WOF peak.This is shown in Fig. 3 for an SIS at y ∈ (1.5, 12), always in the single-image regime.We find the scalings The time-delay scaling is exact for the SIS, but Green's functions dependence is only valid in the WL limit, y ≫ 1.Because the peak in G(τ ) involves short timescales, in the frequency domain it corresponds to high frequencies, i.e. the damped oscillatory pattern after the maximum of F (w).
The broadband shape of the WOF is also sensitive to the impact parameter.This is determined by the behaviour of G(τ ) over large time delays and is therefore captured by low w features of F (w).We can thus use the position and height of the first peak, w 0 and F (w 0 ) respectively, to characterize the broadband WOF.We find the scalings The scaling of w 0 is approximate, but the scaling of F (w 0 ) is very accurate.The likely cause is that the height is dominated by the transition between the behaviour near the image (with a well-defined scaling with y) and the asymptotics of the lensing potential (independent of y).See Eq. (D7) for an analytic estimate for the scalings in the frequency domain.
Measuring a WOF may allow one to infer the magnification of the image, which is not directly observable in the GO limit.This is because y determines µ, as can be seen from either the height of I(τ → 0) or the asymptotic value of |F (w → ∞)|.For the SIS in the single-image regime (this also holds approximately true for lenses with the same large-x behaviour, like the CIS).Therefore, Eqs. ( 30), ( 31) imply the following scalings The correlation between the µ and the WOF properties opens the possibility of mitigating the uncertainty due to WL in standard sirens.We will comment on this potential application in Sec.V B.
We explore the role of the lens compactness and shape by considering a CIS with variable core size x c . Figure 4 shows results at fixed y = 3, but varying x c between 0 (the SIS limit) and 1 (a sub-critical lens, unable to form multiple images even for y → 0).The main effect of x c is on the amplitude and shape of the WOF peak: smaller cores produce narrower and taller peaks that persist at higher frequencies.Larger cores also shift the position of the peak slightly towards lower τ C : in this case the peak is associated to the edge of the core region, where contour curvature is maximum, rather than the lens' center, where the contours are much smoother.In the sub-critical lens case (x c = 1) the peak in G(τ ) is barely recognizable.The broadband WOF is still apparent in F (w) by the onset of diffraction, which is caused by the overall transition rather than the peak.
Let us briefly discuss the prospect of constraining lens parameters from the observation of weakly-lensed GWs.Such an inference is possible in principle, at least if we assume a lens model.If we assume an SIS, we can infer the lens parameters from the WOF peak position and height via Eq.(30).The degeneracy between M Lz and y can be broken because the peak's position τ C = t c /4GM Lz depends on the effective lens mass, while its amplitude G(τ C ) depends only on y.Converting the projected mass M Lz into the halo mass, M vir , requires knowledge of ξ 0 , which depends on the redshift of the source and the lens.While z S can be constrained by the amplitude of the signal, z L ∈ (0, z S ) is generally unknown and only a lower bound on M vir can be derived (corresponding to the largest ξ 0 ).Nonetheless, assuming a halo mass function enables a probabilistic inference of M vir , which is sharply peaked around the minimum possible value (see Ref. [15] Sec.V A for details).Additional leverage on the lens parameters can be obtained from the broadband feature.
More general lens models will make parameter inference from the WOF more challenging.As we saw in the case of the CIS (Fig. 4), the lens' internal parameters affect the height of the peak, complicating the distinction between y and M Lz outlined above.These additional parameters can be constrained by the shape of the peak in the WOF and the broadband feature.Nonetheless, degeneracies with M Lz and y will affect the precision (see Sec. IV in Ref. [15] for examples in strong lensing) and will lead to biases if the wrong lens model is assumed.Because the WOF depends on the entire lens, it is possible to reconstruct ψ(x) given I(τ ) under several assumptions (weak lensing, symmetric lens and known y).We discuss this possibility in Sec.V A.

B. Modelling composite lenses
Let us now address how GW observations may probe a lens with a non-symmetric profile with an internal structure.We will consider a matter distribution composed of N sub objects with a common projected profile, Σ sub , and mass We consider equal-mass sublenses for simplicity, but our expressions are easy to generalize to a mass distribution.The centers of each sub-lens, x i , will be drawn from a distribution P (x i ).Together with N sub , the functions P (x i ) and Σ sub determine the statistical properties of the composite lens.We stress that this composite lens model is not intended to be a realistic realization of a halo, but it provides insights on the WOF produced by substructures.
The average surface density for the composite lens is a convolution of the distribution P (x) with the sub-halo profile Σ sub Because of linearity, an analogue expression can be derived for the lensing potential: 36) has some obvious limits: if Σ sub is a delta function then the average profile is given by P (x) and vice-versa.However, considering two extended functions gives non-trivial profiles in general.It is in principle possible to derive an expression for the variance, higher order statistics and correlations at different points x 1 , x 2 , to further characterize the convergence towards the average lens as N sub → ∞.
For the composite lens, we will assume a distribution of sublenses that follows the SIS projected density profile, and model each of the sublenses as SISs where x is written in polar coordinates x and φ.The truncation at x max ensures a finite average profile (36) for the above x, x ′ dependence.Writing x ′ = ξx and using the above definitions, we obtain the average projected profile and the lensing potential: These expressions can be solved numerically.Let us now explore different realizations of this composite lens and their average limit.

C. Signatures of composite lenses
As we saw, WOFs are characterized by their broadband modulations and the peak's position/delay, amplitude and shape.Each of these characteristics depends on the lens parameters.An important property of a WOF is that it is approximately linear in the projected density Σ(x).Hence, Green's function G(τ ) (or I(τ )) of a composite lens is well described as the sum of the WOFs associated to each of the sublenses, appropriately timeshifted, at least in the WL limit.Therefore, WOF peaks in GW data could be used to identify individual sublenses, infer their spatial distribution and constrain their properties.We will use the composite lens, Eq. ( 35), and its average profile, (36), to investigate how the number of sub-lenses and their distribution leave a characteristic imprint.
We will start exploring the effect of the number of sublenses.Figure 5 shows Σ(x), I(τ ), G(τ ) and F (w) for realizations of the composite lens (35) with N sub = 2, 4, 100 sublenses with the same total mass and at fixed y = 3. Predictions for a single SIS and the averaged lensing profile, Eq. ( 36), are shown for comparison.Information about the sublenses is most transparent in Green's function: each object forms a distinct peak, whose height and position are determined by its mass and separation from the GO image in the lens plane (following the trend seen for symmetric lenses in Fig. 3).When the number of sublenses is large, N sub = 100, G(τ ) follows the averaged profile closely, with some stochastic "cuspiness" added.In the frequency domain, the existence of multiple peaks becomes an interference pattern in the damped oscillations.For N sub = 100 the superposition is mostly incoherent, reproducing the prediction of the average profile.
The broadband shape of the WOF (as characterized by the position and height of the first peak) is independent of the number of sublenses.The relative difference in F (w 0 ) between the models shown in Fig. 5 is ∼ 10 −3 between the SIS and the average lens and ∼ 10 −4 between the average and composite lens N sub .This homogeneity is in stark contrast with the patterns of peaks/damped oscillations for different values of N sub .This consistency may reflect the broadband WOF only depending on the difference between the GO region and the asymptotic behaviour of the lensing potential, as argued above.
Let us now explore the effects of the mass distribution on the WOF. Figure 6 shows Σ(x) and G(τ ) for composite lenses made of four equal-mass SIS profiles, equally spaced along a straight line.This type of sublens distribution is a crude approximation to filamentary structures in the cosmic web.In this case, the position of the four peaks satisfies τ n ∼ 1 2 (y−αn cos θ) 2 +(αn sin θ) 2 , where θ is the angle between the line of sublenses and the optical axis and α controls the spacing.Therefore, the separation of WOF peaks in the time domain could be used to infer aspects of the spatial distribution.This shows that, despite all the information being compressed into a single dimension (e.g.G(τ ), F (w)), it might be possible to reconstruct the morphology of a 2-dimensional (projected) structure under certain circumstances.We will discuss this potential application in Sec.V C.
Just as for individual lenses, identifying separate peaks in G(τ ) can provide information on substructure mass and relative positions, cf.Eq. 30 (although assumptions about the lens profile might be necessary, cf.Sec.III A).In the limit in which the linear approximation holds, G(τ ) only depends on the projected matter distribution of the lens, with the amplitude determined by the overall impact parameter.A high-quality observation of a GW signal thus offers additional constraints on the inner structure of the lens.Identifying the WOFs allows the reconstruction of substructures, although certain degeneracies persist.First, since WO lensing encodes 1dimensional information (through I(τ ), G(τ ) or F (w)) it is not possible to recover the 2-dimensional projected distribution Σ(x).This is obvious from the linear WL limit, in which the position and height of the peak only depend on the offset between the sub-lens and the GO image (for fixed sub-lens profile).Second, lack of knowledge of the lens redshift prevents us from accurately determining the virial mass of the lens (or sublenses), as converting M Lz (observed) into a virial mass requires knowing d eff and thus z S (constrained from the signal amplitude) and z L (unconstrained).As already discussed, this degeneracy can be partially broken using probabilistic information, see Sec.VA in Ref. [15].

IV. OBSERVATIONAL PROSPECTS
Let us now discuss the prospects of observing WOFs.We will first derive the maximum impact parameter at which WOFs can be detected for a given GW source (Sec.IV A).We will then include information about the halo abundances to estimate the probability of observation (Sec.IV B) and the prospect of constraining the halo mass function (Sec.IV C).Our estimate of the probabilities includes only isolated halos: we conclude this section by discussing subhalos and their expected imprints (Sec.IV D).

A. Critical impact parameter
Assessing whether WOFs are detectable requires accounting for the details of GW sources and the instrument's sensitivity.We will assume sources to be equalmass ratio, non-spinning compact binary coalescence and describe them with the IMRPhenomD [28] model in the PyCBC package [29].We will focus on LISA [30,31] and the Einstein Telescope (ET) [32,33].For each instrument, we will include the effect of sky inclination and polarization averaging over the antenna pattern functions [34].Our results will then reflect typical sources: neither optimally aligned nor close to a blind spot of the instrument.We further consider the detector to be static.This is a good approximation since there is a single image and the SNR is very concentrated around the merger (see Sec. IVA in Ref. [26] or Ref. [35]).Finally, we will only obtain results for single detectors: a detector network will in general improve the prospects of detection by improving SNR and sky coverage.
We will assess the detectability of WOFs based on the mismatch between the lensed and unlensed waveforms.For two generic waveforms, h 1 and h 2 , the mismatch is defined as: Here we introduced the noise-weighted inner product for two signals h 1 (t) and h 2 (t) with Fourier transforms h1 (f ) and h2 (f ): where S n (f ) represents the sky-averaged one-sided detector power spectral density.In terms of this product, the signal-to-noise ratio is SNR ≡ (h|h).
According to the Lindblom criterion [36], two waveforms are considered indistinguishable if the condition (δh|δh) < 1 is satisfied, where δh ≡ h 1 − h 2 .When considering signals with comparable SNR, as for a weakly lensed and an unlensed signal, this criterion requires M × SNR 2 < 1.In general, the converse of the Lindblom criterion is not true.Other factors, such as correlations between waveform parameters and bias in the recovered source parameters, may restrict the detectability.Parameter degeneracies can be accounted for using the Fisher information matrix [37].However, Ref. [15] showed that the Fisher matrix can overestimate the precision in lensing parameters due to the breakdown of the linear signal approximation: this issue depends on the lens model and parameters and requires a case-bycase inspection.Detectability in the Fisher matrix approach is usually defined through the standard deviation of the marginalized posterior, leading to results that depend strongly on whether M Lz or y is used (see discussion below).Ultimately, Lindblom and Fisher analyses are answering different questions about the information gained from the signal.Due to computational simplicity and the arguments discussed above, we will employ the flipped Lindblom criterion, considering WOFs to be detectable if M × SNR 2 > 1. Henceforth, M is the lensedto-unlensed mismatch, unless stated otherwise.At the end of the section we will compare our results to other methods.
Weak lensing on a single image can be detected at y ≫ 1 via WOFs.We stress that in GO this is not possible because a time delay and amplitude magnification on a single image are degenerate with source properties (luminosity distance, coalescence time).The condition M(y cr ) × SNR 2 = 1 determines the critical impact parameter, y cr , which characterizes the minimum level of alignment between a lens and source required for detectability.We will omit the explicit dependence of y cr on the source and lens properties.We will assume all lenses to be described by the SIS profile, regardless of their mass.
The critical impact parameter depends on both the lens' and the binary's masses, as shown in Fig. 7. Here, the virial mass, M vir , associated with M Lz is shown on the upper scale, assuming a lens redshift of z L = z S /2 (the relationship between both masses is discussed in App.B).The mass range covered by each curve depends on the frequency range spanned by the waveform and the corresponding detector response.Typically, since the merger frequency of a BBH dominates the SNR and scales as 1/M BBH , the onset of WO is expected at larger lens masses for heavier BBH.This is evident when comparing LISA and ET systems in the Figure.However,  Figure 8 shows the mismatch as a function of y for fixed M Lz and source properties.The mismatch increases below y < y cr , but it saturates at low enough y.Interestingly, the mismatch at y ≪ y cr is lower for the cases with larger y cr : the signals with the more clear signatures (higher M) are not the ones for which the detection probability is most enhanced (high y crit ).The Figure also shows the effect on both the amplitude and the phase of the GW signal.The two contributions oscillate as a function of the impact parameter, as when approaching the WO regime.
The behaviour of the y cr curves presented above can be understood using the approximate analytical form of the WOF for the SIS lens given by Eq. (D7), valid at large w.Here, the WOF arises from the lens center as δF c = i/(wy 3 )e iwϕc .Since we are considering large y values, this is a small correction to the amplification factor, and the critical curves can be derived by expanding the mismatch M of Eq. ( 41) in this quantity.Moreover, we assume the mismatch to be dominated by a single frequency, f ⋆ .This is either the BBH merger frequency at the detector f max (that we take to be double the ISCO frequency, assuming negligible final spin, , see e.g.[38]), or the frequency at the detector's sensitivity curve minimum, f det .
At sufficiently small w ⋆ ≡ 8πGM Lz f ⋆ , the approximation used for δF c breaks down, and one has to resort to the full WO result, Eq. (D6).However, in the regime of small w we expect the lensing effect to be negligible, with F (w) ≃ 1, and the WL critical impact parameter value to drop down rapidly.This condition can be implemented, roughly, by cutting-off the curves at the onset of WO.By inspecting Eqs.(D6) and (D7), this happens when w ⋆ y 2 /2 ≃ mπ, where m is a small integer.This corresponds to a peak value of the curves y max cr ≃ SNR/(2 3/2 mπ), that depends only on the SNR of the signal.
Following these prescriptions, the critical curves can be approximated as follows: where the curve is truncated at the maximum, M max Lz ≃ 32π 3 m 3 /(w ⋆ SNR 2 ), and We find that setting m = 2 when f max < f det and m = 1 otherwise returns the closest match with the numerical results for LISA.Moreover, the agreement is improved if the f max is taken to be the actual merger frequency instead of using the ISCO approximation (for numerical fits cf.Eq. ( 29) in [39]).For ET the noise curves are flatter, making the single-frequency approximation less accurate.Notice that M Lz depends on the virial mass and distance of the lens through Eq. (46).We stress that different mismatch thresholds, i.e.M(y cr ) × SNR 2 = Λ 2 , lead to re-scalings and shifts of the critical curves.In particular, the peak's position is M max Lz ∝ 1/Λ 2 and its value as y cr ∝ 1/Λ.These scalings can be used to extrapolate our results to more stringent detection criteria (Λ > 1).
Let us now compare our results with previous analyses that addressed the detectability of WO effects by LISA by different methods.Reference [14,40] employed a Fisher-matrix analysis, including source parameters.The authors define the critical impact parameter in terms of ∆θ/θ, the ratio of the marginalized posterior width to the fiducial value of the lens parameter θ.Specifically, they require M Lz at the critical impact parameter to be measurable at 1-sigma (the relative error is order one).This criterion can account for potential degeneracy in the lenssource parameters.Table II in Ref. [14] shows results of y cr for an SIS based on θ = M Lz or θ = y in three cases: the two estimates differ by a factor ∼ 2, with ∆y/y giving the larger y cr .Our analysis gives slightly larger values, with y cr = 58, 48 and 26, which are 28%, 29% and 33% larger than the results in Ref. [14] for the same source properties.

B. WO optical depth
Here we forecast the probability for a GW signal to carry a detectable WOF signature.We will focus on isolated matter halos, described by the SIS profile and with a mass distribution characterized by a halo mass function.We will further assume that halos at any given redshift are distributed homogeneously in space.Lensing probabilities are then described by Poisson statistics Here λ is the optical depth, which we define below.k is the number of lenses contributing a detectable signature to the signal: for WOFs, cases with k > 1 would produce signatures similar to those shown in Figs. 5 and  6.The probability of having any number of detectable imprints is simply P(k ≥ 1) = 1 − e −λ ≃ λ, where the approximation holds for λ ≪ 1.
The lensing probability is given by the optical depth λ evaluated at the source's redshift.The total optical depth is given by an integral over halo masses λ = d log(M vir ) dλ d log(Mvir) , where the differential optical depth per logarithmic virial mass is Here n L is the number density of lenses at fixed M vir , χ L is the lens' comoving distance.The integrand consists of two factors.The first is the lensing angular cross section, namely πθ 2 cr , representing the projected area where WOF is detectable.For a symmetric lens θ cr = y cr ξ 0 /D L , where the critical impact parameter is the one discussed in the previous section.The second factor in Eq. ( 45) is the halo mass function, dn L (M vir , z L )/d log M vir , i.e. the number density of lenses per logarithmic virial mass unit.We adopt the Tinker halo mass function [41] extrapolating it to low halo masses, i.e.M vir < 10 10 M ⊙ (other halo mass functions produce similar results).We note that our dn/d log(M vir ) only accounts for isolated halos: subhalos within the larger structure may further contribute to the probabilities presented here.We will comment on subhalo signatures in Sec.IV D.
The optical depth depends on two definitions of the mass: M vir in dn/d log(M vir ) and M Lz , Eq. ( 4), in y cr .They are related by: , (see Eq. (B4) for further details).Note that both quantities differ by several orders of magnitude for typical halos.The above relation depends on the source and lens' redshift.z S can be inferred by the amplitude of the signal, which scales as ∝ D L (z S ) −1 , assuming a cosmology.As we mentioned multiple times, while z L is in general unknown, it is possible to place some restrictions on M vir , given an observed M Lz , the quantity that lensing is sensitive to.Let us give more details about this procedure.Assuming a halo mass function, one can assign a probability distribution to P (M vir |M Lz , z S ).Because d eff is bounded from above, it defines a minimum value M min vir (at fixed M Lz ).The probability P (M vir |M Lz , z S ) is sharply peaked around M min vir although larger values are possible, see Sec.V-A in [15] for details.
Using the definition of the effective lens mass, Eq. ( 4), and dχ , the differential optical depth can be recasted as: We will consider two types of lensing observation: • λ WO for WOF detection, with y cr computed as described in Sec.IV A.
Note that λ WO depends on z S as well as the source properties (M BBH ), while λ SL only depends on z S .Note that the strong-lensing criterion is based solely on the formation of multiple images for the SIS model.Requiring that the secondary image is brighter than a detection threshold would reduce λ SL .However, lensing situations with y ≳ 1 could still produce strong-WO phenomena even in the single-image regime [7].The probability of these situations is approximately represented by λ SL .Figure 9 shows the differential optical depth as a function of the virial mass of the lens, for various redshifts and total BBH masses (source-frame), within LISA's reach.The higher peak is reached by the 10 6 M ⊙ curve at z S = 5, which corresponds to the largest critical impact parameter and SNR among those considered.We observe that larger source redshifts do not necessarily reduce the differential optical depth.The SNR is reduced by a factor (1 + z S )D S , while the number of intervening halos and their typical Einstein radii increase with z S .In addition, as the signal gets redshifted, i.e.M D BBH = (1 + z S )M BBH , the merger can be displaced to a frequency band with lower/higher sensitivity.The final behaviour depends on M BBH .For instance, for M BBH = 10 7 M ⊙ , the closer source curve (z S = 5) encompasses the further one (z S = 10), while the opposite is true for M BBH = 10 5 M ⊙ .The optical depth for strong lensing is independent on the BBH mass and increases with redshift as expected.We plot the total optical depth λ in Fig. 10 for typical sources detectable by LISA (upper panel) and ET (lower panel).Physically, this quantity represents the average number of lenses per source entering the Poisson distribution, Eq. (44).
We find that that λ(z S ) displays a nontrivial dependence with M BBH : it can increase or decrease with z S , as a result of the competing effects discussed above (SNR, lens projection and frequency shift).In the 10 2 and 10 3 M ⊙ curves for ET we can see the effect of the signal falling out of the detector's horizon at high z S .
The prospect of observing WOFs is optimistic: our mismatch analysis shows that ∼ 20% of LISA binaries with M BBH ∼ 10 6 M ⊙ would carry observable WOF signatures.WOF-detection by LISA is substantially more likely than strong lensing, except for the heavier sources at very high redshift (M BBH = 10 7 M ⊙ , z S > 7).While WOF detection will be dominated by events where a single lens contributes (k = 1 in Eq. 44), the probability of detecting WOF signatures from multiple lenses ∼ λ 2 is non-negligible for M BBH ∼ 10 6 M ⊙ , even when considering only isolated halos.For ET, strong lensing is more likely.Still, WOFs can plausibly be observed given the large number of expected sources, with rates ∼ 10 4 − 10 5 yr −1 [42,43].
We conclude by recalling that, while the flipped Lindblom tends to overestimate the critical impact parameter, two assumptions make our estimate conservative.First, y cr was computed from averaged noise curves.This accounts for typical sources, but we expect it to underestimate λ WO , as λ ∼ y 2 cr ∼ SNR 2 and ⟨SNR 2 ⟩ > ⟨SNR⟩ 2 , i.e. events with better source/detector alignment at fixed z S will compensate for the fainter ones.This information (e.g. based on the posteriors for the sky localization, binary inclination, etc...) could be folded into the estimates discussed in Secs.IV C and V D. Second, the halo mass function only accounts for isolated halos.Light subhalos (M ≲ 10 9 M ⊙ ) are more abundant than their isolated counterpart, and will further contribute to the optical depth.They may also be distinguishable through characteristic WOF signatures, see Sec.IV D. Now, we compare our estimates with previous results in the literature.The estimates in Ref. [40] are less optimistic due to the differences in the critical impact parameter definition, as discussed at the end of the previous section.Ref. [44] obtained probabilities almost 4 orders of magnitude lower.In this case, the difference can be traced back to different modelling of the source and the lens.They use the lowest order post-Newtonian expansion of the waveform, truncated at the innermost stable circular orbit, which covers only the inspiral stage of the coalescence.As discussed above, the peak of the critical curves is associated with the large SNR contribution coming from the merger.Neglecting the merger can severely underestimate the height of the peak in y cr and its position.Finally, they work with a Navarro-Frenk-White lens profile, which is shallower than the SIS towards the center and consequently produces a less prominent WOF.Along these lines, Ref. [45] finds larger detection probabilities for cuspier profiles.This analysis is also based on inspiral-only waveforms and thus leads to low detection prospects.The impact of the lens profile on the detection prospects will be addressed in future work.

C. Total rates: probing the halo mass function
We can now estimate the total number of events with observable WOF and use that information to constrain the halo mass function on the mass range where WO effects are observable.We will focus on LISA, for which the optical depth was found to be significant.
The results depend on the total number of detectable signals, their redshifts and mass distributions.For LISA, a recent analysis placed bounds on the detection rate compatible with results from pulsar timing arrays [46].We will use the results for their "agnostic" [47] and "astro-informed" [48] models (obtained from their Fig. 2) for sources with z S ≤ 5 and M BBH ≥ 10 6 M ⊙ .We will estimate the rate of lensed events (WL, SL) as λ L (z * S , M BBH ) • Ṅdet (M BBH )T obs , integrating over the binary masses.We evaluate the optical depth z * S = 5, the upper end of their interval.This is conservative for weaklensing, see Fig. 10 (λ WO increases with z S for lighter sources M BBH = 10 5 M ⊙ , which are not included).The choice is optimistic for strong lensing/multiple images, as the optical depth grows with z S .
Table I shows the rates of detected events, WO distortion and multiple images under the above assumptions.Each row corresponds to a confidence interval within the agnostic/astro-motivated models.The detection rates vary several orders of magnitude, especially in the agnostic model.Note that the above numbers for WOF detection neglect sources at high redshift (z S > 5) and with light masses (M BHB ≤ 10 5 M ⊙ ), whose inclusion would increase the detection rate.Even so, multiplying by the observation time T obs ∼ 5 yr gives reasonable observation prospects in all but the most pessimistic cases.
Detecting WOF enables novel probes of the halo mass function.In the simplest cases, constraints rely on the ratio of observed ṄWO / Ṅdet to the values predicted in Eq. (47).To perform a quantitative estimation of the sensitivity, let us consider a constant rescaling of the halo mass function where the fiducial case f H = 1 corresponds to the case discussed in Sec.IV B. This represents a phenomenologi-cal parametrization in the range of M vir for which a detector is sensitive to WOF, cf.Fig. 9.
Let us now derive limits on f H given the observation of k det events, of which k WO carry WOF signatures.We will not consider the possibility of events with multiple WOFs, with probability O(λ 2 WO ).Our theoretical model depends on N det = r Ṅ (0) det T obs and N WO = f H Ṅ (0) WO T obs , where T obs is the survey time and the free parameters f H , r encode variation with respect to fiducial values Ṅ (0) X .The likelihood of the data (k wo , k det ) given the model (f H , r) is where P(k|λ) is the Poisson distribution, Eq. ( 44).We obtain the 1-dimensional confidence levels in f by evaluating the posterior on a grid in r, f and integrating r.This implicitly assumes a wide prior, i.e.Π(r, f ) ≃ 1 around the peak of the likelihood.For k WO > 0, the limits are consistent with estimating the posterior by sam- where k WO , k det are taken from a Poisson distributions with rates N WO , N det .
Table I shows the marginalized 90% confidence interval on f H that can be achieved by a T obs = 5 yr LISA mission, given the rates shown in Table I.We have assumed that the number of events k det , k WO is ⌊N where ⌊x⌋ is the floor function.The sensitivity ranges from sub-percent in the most optimistic case to 1-2 orders of magnitude upper bounds in the most pessimistic cases.Note that cases with k WO = 0 allow upper limits on f H , as long as k det ≤ 1.In particular, the Astro 25% c.l. scenario has a single detection, leading to a O(10 2 ) limit.This is a very simple estimate of prospective constraints on the halo mass function.In Sec.V D we will discuss how more information on dn/dM vir can be obtained from WOF observations.

D. Subhalos vs isolated halos
The rates presented above are conservative, as the halo mass function accounts only for isolated halos, excluding structures that have incorporated into more massive halos, i.e. subhalos (see discussion in Sec. 4 of Ref. [49]).While some subhalos will be disrupted, those that survive may contribute to the probability of detecting WO effects, increasing the rates presented in Sec.IV B.
Besides increasing the probability of WOF detection, subhalos differ from isolated halos in several regards.Even when comparing objects of equal mass, we expect the following effects in subhalos.I: LISA rates and prospective constraints on the halo mass function.The first three columns show rates of total mergers, WO detections and multiple-image events with MBBH ≥ 10 6 M⊙, zS ≤ 5. Rows correspond t the agnostic and astro-informed models in Ref. [46], at different confidence levels (Agnostic 5% c.l. is not shown, as it is several orders of magnitude lower).
The last two columns show the lower and upper limits on the amplitude of the halo mass function, Eq. ( 48).
3. In some cases, multiple strong-lensing images are produced by the main halo.Then, each image may contain WOF of different subhalos.
In addition, subhalos will differ from isolated halos (e.g. in their shape) due to their assembly history and interactions with the main halo and other subhalos [50].These evolutionary features will affect the probability of subhalo lensing and may have observational imprints.
Massive halos contain a large number of lighter subhalos [50].A signal propagating through a galactic-scale halo can have an enhanced probability of encountering multiple low-mass subhalos that can produce a WOF, increasing the detection rates computed for isolated halos in Sec.IV B. The subhalo lensing probabilities can be modelled using a Poisson distribution with a rate NWO,s , which depends on the point where the image forms in the lens plane and the details of the subhalo mass and spatial distribution [50][51][52].Lens configurations for which NWO,s ≳ 1 are likely to produce rich WOF with multiple peaks, as discussed in Sec.III C.
Another effect of substructures is the distortion of the lensing potential by the main halo.The leading order corrections, convergence and shear, have been shown to enhance the diffraction pattern appreciably [7] and can allow lighter microlenses to contribute to diffraction effects [8] (although see [12]).They may thus increase the prospects of detecting the WOF, particularly sources whose image forms in the inner part of the halo.The imprint of convergence and shear may also serve to distinguish subhalos from isolated halos and place them within their host halo.
For closely aligned systems, the main halo will split the source into multiple images.The associated probability is given by the strong lensing rate (λ SL , for y cr = 1) in Sec.IV B. In this case, each image may contain information from nearby subhalos, which will also be affected by convergence and shear from the main halo.Identifying multiple GW events as strongly lensed images from the same source (e.g. by overlapping sky-localization and in-trinsic parameters) would provide a unique opportunity to constrain the properties of subhalos.
Finally, we note that WOFs could also be imprinted by close-by systems, such as subhalos of our own galactic halo.The critical curves presented in Fig. 7 do not depend on the lens position and retain their validity when z L , d eff → 0. However, the physical impact parameter shrinks -the Einstein radius goes as √ d eff -and the conversion between M Lz and M vir gets offset.For example, inverting Eq. ( 46) the interval M Lz = 10 −1 − 10 3 M ⊙ is mapped into M vir ≃ 10 6 − 10 9 M ⊙ (at z S = 1 and D L = 200 kpc).We stress that the geometry of the system does not suppress the lensing probability.Indeed, despite the contraction of the physical impact parameter, the quantity χ 2 L θ 2 cr in Eq. ( 45) is invariant.If a WOF is observed, sky localization information could be used to determine the probability of local versus cosmological origin.

V. POTENTIAL APPLICATIONS
We will now discuss some possible uses of WOF to constrain properties of the large-scale structure.Our presentation will be qualitative and schematic.More detailed analyses are left for future work.

A. Lens reconstruction
The simple expressions behind the perturbative weaklensing framework (Sec.II C) open the possibility of systematically reconstructing lens features.As already explained in the previous sections, a full reconstruction of the 2-dimensional lensing potential is impossible from 1-dimensional data from a single source (G(τ ), I(τ ) or F (w)).
Nonetheless, assuming that the lensing potential is symmetric, ψ(x) → ψ(x), a formal relation between the time-domain integral and the lensing potential can be obtained from the leading-order term in Eq. (23a).First, we change integration variable from φ to x and obtain with a kernel where φ is a function of x, y and τ (according to Eqs. (23b), (23c)), and x 0 , x π are the values at which φ(x, τ, y) = 0, π, respectively.By knowing I(τ ), Eq. ( 50) gives the lensing potential as the solution of an integral equation.One can then obtain the projected mass density Σ(x) e.g.Eqs.(11) in Ref. [15].
While potentially interesting, it is not clear how useful the above expression may be in practice.Besides assuming linearity and symmetry of ψ, Eq. ( 50) requires knowing the impact parameter y.If this value is not constrained (e.g. from the WOF peak), one possibility is to perform the reconstruction for different values of y and consider the most plausible reconstructed lens, according to some prior, e.g. from theory or simulations.Ultimately, the reconstruction will be limited by how well I(τ ) can be inferred from real data.

B. GW delensing
The correlation between the WOF and the GO magnification opens the possibility of inferring |µ|.Eqs.(33) and (34) give approximate relationships for symmetric SIS models.
This information would enable a re-calibration of the intrinsic luminosity of the source, a major source of uncertainty for standard sirens at high redshift [53].Such procedures are known as "delensing".Proposed methods for GW delensing usually rely on EM follow-up to characterize the lensing potential in the direction of the source [54,55].Because GW sources are poorly localized, EM follow-up methods need to cover a large portion of the sky, which can become very costly.
The main limitation is that magnification is dominated by galactic scale lenses, too heavy to produce an observable WOF by themselves.Hence, WO-based delensing would rely on detecting substructure within the main lens (cf.Sec.IV D).This may limit the applicability to a subset of sources.However, delensing based on the WOF would not require additional observations and could be attempted without any costly follow-up.This is similar to delensing of the cosmic microwave background, which can be performed at the level of the observed maps [56].

C. LSS morphology
The large-scale structure (LSS) of the universe displays a rich pattern in the distribuiton of dark and baryonic matter.Different universe regions can be broadly classified by their morphology, determined by the number of independent directions that are expanding vs contracting.This gives rise to 4 categories: voids (3 expanding directions), sheets/walls (2 expanding, 1 contracting), filaments (1 expanding, 2 contracting) and halos (3 contracting) [57][58][59].
In Fig. 6 we explored an idealized representation of a filamentary structure.Our model, a 1-dimensional chain of equal-mass and equally-spaced lenses, left a characteristic imprint in Green's function, with a clear dependence on the angle between the chain and the location of the source.Weakly-lensed GWs have the potential to distinguish between these patterns, allowing not only to identify a filament but also to reconstruct its angle, and perhaps even other properties (e.g.mass and relative spacing between sub-lenses).
Realistic realizations of the LSS are vastly more complex.Nonetheless, information about the lens morphology will be present in the WOF, e.g. through the statistics of the peak distribution in G(τ ).This information is projected from a 3D distribution into the 1D Green's function, and thus a complete reconstruction is not possible (see Sec. V A).Nonetheless, it might be possible to obtain information on the lens morphology, e.g. in high-SNR observations where multiple peaks can be clearly located.More likely, morphology reconstruction will, at best, assign a probability to each category given an observation.
Lensed GWs may offer a complementary means to identify light halos and subhalos, distinguish between both and constrain their properties [15,44,45] and abundances.Isolated light halos at cosmological distances can be observed by future detectors thanks to the large critical impact parameter (Fig. 7).In Sec.IV C we showed how LISA can obtain constraints on the halo mass function from detection of WOFs or their absence.In the most optimistic case, a constant rescaling on the scale of interest f H , Eq. ( 48) can be constrained to sub-percent level, while a single unlensed event yields an O(10 3 ) limit.As argued in Sec.IV B, including information on the source's parameters that affect the SNR (sky localization, inclination) is likely to produce more robust limits.Information on the impact parameter posterior can also be incorporated; see Ref. [77].
A more detailed analysis can also constrain the virial mass dependence, i.e. f H (M vir ), using the fact that the source mass M BBH determines the range of M vir that can be probed (cf.Fig. 9).However, as the rates are dominated by the optimal mass (∼ 10 6 M ⊙ for LISA), constraints on the mass dependence will be far less stringent than the overall amplitude.It may also be possible to constrain the redshift dependence: the rates for different M BBH evolve differently with z S (cf.Fig. 10) and z S is known from the luminosity distance (assuming a cosmology).While important degeneracies are expected in a generic f H (M vir , z L ), it may be possible to test models with suppression of light halos (e.g.warm dark matter [78]) or a minimum low-mass cutoff (e.g.ultra-light dark matter [79,80]).
Additional information on the lens may be obtained via the WOF: the relevant observable quantities are the source's redshift (from the signal's amplitude, assuming a cosmology) and the lens' effective mass M Lz (see discussion in Sec.III A).Because the lens redshift is unknown, M Lz only provides a lower bound on the virial mass of the halo.However, assuming a halo mass function allows one to define a probability distribution for z L and M vir , which is peaked around the minimum M vir (see Sec. VA in Ref. [15] for details).The inferred values from all events with WOF would then improve on the constraints from detection counts.This method bears analogy to inferring the properties of galactic-scale lenses from strongly lensed signals using the distribution of time-delays between multiple images [81].
Lensed GWs can also probe subhalos and distinguish them from isolated halos.As argued in Sec.IV D, the expected signature of substructures is the detection of multiple peaks in the WOF (from nearby subhalos) and/or the presence of convergence/shear caused by the main halo.Subhalos may have different abundances and properties from isolated halos of the same mass (e.g.due to tidal stripping, shock heating and other interactions [50]).Being able to distinguish both populations separately may offer valuable clues about the assembly history of the large-scale structure.

VI. CONCLUSIONS
We have investigated the phenomenology of gravitational lensing in the single-image, wave optics (WO) regime, with an outlook on their potential to probe cosmic structures and prospects for observation by future GW detectors.Our results converge to the well-known limits of geometric optics (GO) in the large source fre-quency/lens mass limit.Large angular separation between source and lens corresponds to weak-lensing (WL), where WO corrections are subtle but potentially observable.
We presented two methods to solve the diffraction integral in the time-domain, adapted to the single-image regime but accounting for WO effects (Sec.II).Both approaches yield accurate results in their domain of validity.First, we present an algorithm able to explore any singleimage configuration, its accuracy limited only by numerical errors.This method is valid even in strong-lensing: it can be used to explore the outer regions of a caustic, where the GO magnification diverges and WO effects persist at high frequency.Second, we develop a perturbative expansion on the lens potential.This method is faster and converges very rapidly to the full solution in the WL limit, at large impact parameter y.The leading-order expansion is linear in the lensing potential/projected lens density, making the study of composite lenses straightforward.Both methods are fast enough for applications such as parameter estimation.
Using these algorithms, we analyze the phenomenology of Wave-Optics Features (WOFs), the WO imprint on lensed GWs (Sec.III).This discussion is particularly clear when using Green's function G(τ ).The most salient aspects of the WOF are its peak and broadband distribution.The peak forms at the center of the lens: its associated time delay and height contain information about the lens location (relative to the GO image) and the lens mass (Fig. 3), while its shape is related to the lens projected density (Fig. 4).The broadband profile is related to the large-scale properties of the lens, such as its total mass and spatial extent.In the frequency domain, the broadband feature corresponds to the first maximum of F (w), while the peak appears as a damped oscillatory pattern at higher frequencies.Our analysis also applies to composite lenses: we study the superposition of N sub equal-mass sublenses with an SIS profile.Each sublens produces a peak, whose associated time delay and amplitude are set by its distance to the GO image in the lens plane.We computed the average profile for this distribution and showed that G(τ ) and F (w) converge to it in the limit of large N sub , although finite N sub is associated to stochasticity in Green's function.
We then address the prospects of detecting WOFs (Sec.IV).Future detectors can potentially observe WO effects at impact parameter one to two orders of magnitude larger than the Einstein angle.Assuming a halo mass function (i.e.extrapolating to the mass-scales producing WO) and uncorrelated spatial distribution, our results can be directly translated into detection probabilities (via the optical depth).The prospect of detection is optimistic for LISA, with rates ∼ 20% for mergers with the optimal mass range.The rates suggest that detecting WO effects is plausible, although subject to uncertainties on the LISA detection rate for MBHB mergers: the number of detections ranges from large to unlikely in the scenarios we have considered.Detection proba-bilities are lower for ground detectors due to the reduced SNR and higher frequency, as WOFs require lighter halos with smaller Einstein radii.Finally, we note that these rates refer only to isolated halos: subhalos assembled into larger structures will increase these rates, and might be distinguishable from isolated halos in some cases.
We can summarize our main findings as follows • Gravitational lensing imprints WOF: frequencydependent modulations that can be observed in the waveform without the need of a counterpart, higher-mode emission or association of multiple events.The Green's function offers a particularly transparent way to analyze these features.
• Features of the lens (effective mass, spatial distribution, inner structure) translate cleanly into features of the Green's function (broadband shape, location and height of peaks), thanks to (approximate) linearity in the WL regime.Distinguishing these features would allow constraints on the properties of individual sublenses and their relative positions.
• Our framework explains clearly how a macroscopic lens arises effectively from a superposition of smaller objects.In the frequency domain a large number of sub-lenses add-up incoherently, suppressing the signatures at high frequencies, as in the average smooth lens (Fig. 5).
• Events with high SNR enable detections of WOFs at impact parameters much larger than the Einstein radius.For low-frequency detectors such as LISA, sensitive to heavier halos, this translates into promising prospects for observations, with probabilities well beyond those of strong lensing.
• Observation of WOFs (or lack thereof) constrains the amplitude of the halo mass function in the range 10 5 M ⊙ ≲ M vir ≲ 10 8 M ⊙ , with precision between percent-level to order-of-magnitude upperlimits (cf.Table I).This information will enable constraints on halos that are both poorly constrained (from simulations and observations) and sensitive to the properties of dark matter.
The transparency offered by WO effects to probe largescale structures and the prospect for detection suggest several applications (Sec.V).GW data may allow a reconstruction of the lensing potential under the assumption of a symmetric lens.In some cases, identifying WOFs on a GW signal might be used to infer the magnification of the signal, mitigating a major uncertainty for standard sirens.Novel probes of large-scale structure could be developed: for instance, identifying several peaks in the WOF may serve to constrain the morphology of the gravitational lens.Finally, constraints on the abundance of subgalactic halos could be improved significantly thanks to the information on the lens (projected mass, impact parameter, etc...) obtained from the WOF.
In addition, we envision future directions regarding the computational framework, lens modelling and data analysis.
Our computational frameworks are flexible, accurate and efficient.While we have focused mostly on weak lensing, our methods can be readily applied to single-image strong lensing, e.g. a source very close to a caustic on the side in which a single GO image forms.This regime has WO features extending to very high frequencies, and could be used to probe the phenomenology of strongly lensed GWs without the challenges of including multiple GO images [25,26].Another interesting extension is microlensing by extended structures, i.e. considering the effects of lens sub-structure on a macro-image.Ultimately, our algorithms will be integrated in the "Gravitational Lensing of Waves" (GLoW) code for public use by the scientific community. 2nother important extension is improved lens modelling.Our analysis largely relied on SIS.This choice was motivated by both simplicity and this lens' particular stance in the single-image regime: the SIS is the "cuspiest" profile that does not form multiple images for arbitrarily large y.It leads to the sharpest possible peak in the WOF, which becomes smoother in the presence of a central core.Other, well-motivated, symmetric profiles such as NFW [82] typically have shallower inner cusps.It is interesting to consider these well-motivated distributions, as well as profiles motivated by dark matter theories [15].Beyond symmetric lenses, future models should include elliptic matter distributions, external convergence and shear and realistic realizations of substructure.Besides more complex lens models, addressing observational prospects will require reliable halo mass functions for M vir ≲ 10 8 M ⊙ halos and detailed predictions for source rates (e.g.Refs.[83][84][85][86][87][88][89]).
Important open questions remain on the signal analysis.A limitation of our analysis is the use of mismatch as a detectability criterion, which is optimistic (e.g.ignoring parameter degeneracies).Future analyses should rely on Fisher matrix or Bayesian sampling.In addition, it is necessary to address the issue of false alarm triggers, i.e. detector noise mimicking a WOF signature.These mundane effects will affect the detection prospects, although less significantly than uncertainties on the source rate, cf.Table I.A standing open question is how to optimally identify and analyze WOFs from GW signals.A non-parametric reconstruction of the Green's function from strain data could provide a partial reconstruction of the lensing potential and hence the projected mass (Eq.( 50)).These insights could be combined with methods for lens-agnostic analyses [90].Given the chance of signals being affected by multiple WOFs (but with an unknown number), it will be necessary to devise a framework for data analysis that does not assume a fixed number of lenses, e.g.along the lines of reversible-jump samplers [91].Regardless of identifying WOF, it is important to prevent unaccounted-for signatures to bias the remaining parameter estimation, e.g. in the LISA global fit [92], or misinterpreting those residuals as new physics (e.g.violation of Einstein's GR).
There is a promising future for GW lensing at the intersection between the WO and WL regimes.While subtle, WOF may be common enough to offer a window into cosmic structures and their properties, and will likely lead to applications beyond the ones outlined here.LISA stands out as a particularly promising probe, thanks to the combination of large SNR and low-frequency sensitivity.If deployed, other proposed space detectors will yield even more impressive results, thanks to larger SNR [93,94], a lower frequency band [95] or both.These observations will provide novel means to probe low-mass halos that are both notoriously elusive and a prime test-bed for dark matter models.WOF may thus become a powerful probe of large-scale structure and fundamental physics thanks to the next generation of GW detectors on the ground and in space.
GO coefficients µ, ∆ (1) and ∆ (2) appearing in Eqs. ( 28), (27) for axi-symemetric lenses.These quantities are given by Here we defined a ≡ and all quantities are evaluated at the location of the image, x = x m .As far as we know, the coefficient ∆ (2) was not given in the literature before.

Appendix B: Symmetric lens models
Here we provide some details of the symmetric lens models considered in the text (SIS, CIS) and the relation between the effective lens mass and the virial mass of the halo.

Singular Isothermal Sphere
A commonly used approximation to the density of a halo is given by the SIS profile where σ v is the velocity dispersion of the halo.For this lens, a convenient choice for the arbitrary scale ξ 0 in Eq. ( 4) is The lensing potential ψ(x) associated to ρ(r) then becomes very simple ψ(x) = x.
In the GO limit, the SIS can have one or two images depending on whether the impact parameter is outside or inside the caustic y cr = 1, respectively.For y < y cr , two images form (a minimum and a saddle, labelled respectively by (+) and (−)) with magnifications µ ± = 1/y ±1, time delays ϕ ± = ∓y − 1/2 and Morse phases n + = 0 (minimum), n − = 1/2 (saddle).Only the image corresponding to the minimum survives for y > y cr .The minimum time delay is given by ϕ m (y) = −y − 1/2.

Cored Isothermal Sphere
We will also consider a CIS, a variant of the SIS in which the presence of a central core smoothes the density profile [96,97]: where ρ 0 is the central density and r c is the core radius.The surface density is Σ(ξ) = πρ 0 r 2 c / ξ 2 + r 2 c .Choosing a normalization scale ξ 0 = 2πρ 0 r 2 c /Σ cr , gives the lensing potential where x c ≡ r c /ξ 0 .Similarly to the SIS, multiple images form for sources within the caustic y rc (x c ) ≤ 1, which is smaller than for SIS.An additional requirement for multiple images is x c < 1/2, so the lens' central density Σ cr .The eventual additional GO image is associated to the maximum of the Fermat potential, and forms close to the center of the lens.In the SIS limit x c → 0, the GO magnification vanishes and this image is replaced by the cusp feature.This lens and the properties of the central image are discussed in detail in Ref. [15] Sec.IIIC.

Relation between MLz and Mvir
For extended lenses, the redshifted lens mass defined in Eq. ( 4) does not typically coincide with the physical mass of the halo.The two quantities can differ even by a few orders of magnitude.As a definition for the physical mass we consider the virial mass M vir , defined as the mass up to the virial radius r vir , i.e.M vir ≡ 4π rvir 0 drr 2 ρ(r) (see Ref. [15] for the full expressions for SIS and CIS lenses).On the other hand, the virial radius is defined as ρ(r vir (z L )) ≡ ∆ c ρ c (z L ), with ∆ c (z L > 1) ≃ 18π 2 and ρ c being the critical density at redshift z L , ρ c = 3H(z L )/(8πG).Due to these relations, the virial mass is a function of the lens redshift [98].
In the case of the SIS, with ξ 0 as in the previous subsections, M Lz and the virial mass are related by: Similar expressions for the CIS profile can be adapted from Sec. III-C-1 in Ref. [15].and where Π and K are, respectively, the complete elliptic integrals of the third and first kind, see e.g.Ref. [99].

Frequency domain result
At this point, one can move to the Fourier transform of I(τ ), so to obtain the amplification factor F (w). Starting from Eq. ( 2), we could not find a closed-form expression for F (w), but we were able to reduce it to a single angular integral.Using polar coordinates again, we can solve the radial integral, obtaining F (w) = e iw(y 2 /2−ϕm) 1 + Other representations for F (w) are also available.See, for instance, [26] for a series representation.
Appendix D: Weak lensing expansion for the SIS In this appendix, we will apply the weak-lensing expansion of Sec.II C to the case of the SIS profile.Thanks to the simplicity of this model, we can obtain explicit expressions and characterize their WOFs.= 2 (a + y) a 2 E (q) + 2 (a − y) a 2 K (q) .(D1) Here we defined q ≡ 4ay/(a + y) 2 , a ≡ √ 2t, and t = τ + ϕ m while the functions K(q) and E(q) are the complete elliptic integrals of the first and second kind, respectively.
The functions K(q) and E(q) are real for q ≤ 1.Notice also that the parameter q as a function of t spans from 0 to 1.The latter value is attained at t = y 2 /2, while zero is reached asymptotically for large t.We notice that the function (D1) is always positive.More specifically, it decays for large τ and has a peak for time delays corresponding to the center of the lens (i.e.around t = y 2 /2).
Although there is no analytic expression for the location of the maximum, we can see that it is located around q = 1 (here E(q) has a maximum).This is physically reasonable since it corresponds to the features due to the center of the lens.We can actually find a good approximation for the maximum by expanding I (1) around t = y 2 /2.To do so, we first write q = 1 − ϵ and expand for small ϵ: Setting the derivative with respect to ϵ of this expression to zero yields an equation for the maximum, that can be solved.Such value ε does not depend on y, and is approximately ε ≃ 2.1 × 10 −3 .Having found the approximate maximum q = 1 − ε, we can translate to t.
With these values we can evaluate I (1) (τ ) at the peak, I (1) (τ ) ≃ 4.22/y.We can actually have an analytic form for the peak, with a very good approximation of the tail at large τ .To get it, we expand the elliptic integrals around q = 1, and obtain This expression well approximates the full I (1) but becomes unreliable at small t, where it diverges.In the latter region, one can use the GO expansion instead.It is also interesting to evaluate Green's function at the peak.
From the definition of G(τ ), Eq. ( 9) and by expanding Eq. (D3) around the peak, we obtain We can also analyse the asymptotic behaviour for large τ .This limit maps to the small-frequency limit since we are considering large time delays from the image.At leading order for large τ , the radius is approximated by x ≃ √ 2t and is independent on the angle φ: see Eqs. (23b) and (23c).Then, taking the τ derivative in Eq. (23a) inside the integral we obtain Here ′ stands for the derivative with respect to x.The result above applies to all symmetric lenses, particularly the SIS, where ψ ′ ( √ 2t) = 1.We conclude that the falloff of I (1) (τ ) as a function of τ is related to the asymptotic properties of the lensing potential at large radii.In particular, more compact lenses have faster falloffs (e.g. for a point lens ψ(x) ∝ log x, the decay is I (1) (τ ) ∝ τ −1 ).

Frequency domain result
As we discussed in the previous paragraphs, the center of the lens can be responsible for a peak in Green's function in the time domain, even in the absence of an image.This is indeed the case for the SIS.In the section, using the frequency-domain approximation derived in Eq. ( 26), we obtain the analogous feature in the amplification factor.For the SIS, it is straightforward to integrate Eq. ( 26) directly: F (w) ≃ 1+iwy − πw 2 e iz−i π 4 2zJ 1 (z)+(2iz −1)J 0 (z) , (D6) where z ≡ wy 2 /4 and J ν (z) are the Bessel functions of the first kind.We can better appreciate the effect of the center by taking the high-w limit of the expression just obtained: where |µ| ≃ 1+2/y at large y and the time delay of the center is ϕ c ≡ y 2 /2.We can recognize the second term, going as ∝ w −1 , as the approximate bGO correction from the image (it does not contain phases/time delays with respect to the magnification term).On the other hand, the third contribution containing ϕ c originates from the center of the lens.In the SIS case, the dependence on w of these two contributions scales in the same way, both in w and y.Moreover, we notice that the bGO contribution is suppressed by a factor of 8.

FIG. 1 :
FIG.1: Computation of the amplification factor in the single-image regime for an SIS with impact parameter y = 3.Each contour of constant Fermat potential, ϕ(x) = τ , contributes to a point in I(τ ), see Eq. (13).Green's function G(τ ) is then computed as the derivative of I(τ ), see Eqs. (9) and(11).Finally, the amplification factor F (w) is the (inverse) Fourier transform of the full Green's function G(τ ).The sharp peak in Green's function is associated with the center of the lens, which features a cusp in the SIS.In our case, the lens is located at (0, 0) and the corresponding contour is represented with a dotted line.

FIG. 2 :
FIG. 2: Diagram for computation I(τ ) in the weak-lensing approximation.In the large-time delays region, contour of approximately equal τ are centered at y, have radius 2(τ + ϕm) and angle φ with respect to the x1 axis.The lens is indicated by the blue dot, at the origin.The red star corresponds to the source, appearing at xm.

FIG. 3 :FIG. 4 :
FIG.3: Role of the impact parameter y on an SIS lens.The columns show the projected density, time-domain integral, Green's function and amplification factor.Green's function omits the GO contribution and has been rescaled by y 3 so the WOFs can be appreciated at large impact parameter.Dashed lines show the ϕ(x, y) = τC contours passing through the lens' center and the associated values of τ .

FIG. 5 :
FIG. 5: Convergence to the homogeneous lens when the number of sublenses is increased.The columns show the projected density, time-domain integral, Green's function and amplification factor for 2-100 SIS.Note that substructures in a realistic halo would appear far more scattered.

FIG. 6 :
FIG.6: Predictions for a chain of 4 sub-lenses as a function of its orientation θ with respect to the optical axis aligned with the source (orange star).Only projected density and Green's function are shown.The sublenses appear as a series of peaks in G(τ ), whose distribution over τ depends on the orientation.A chain perpendicular to the optical axis (lower right) produces closely overlapping peaks, which have a distinct shape from the SIS prediction.The image position (red cross) at different orientations differ marginally.

FIG. 7 :
FIG.7: Critical impact parameter curves plotted against the lens mass MLz (bottom axis) and the virial mass Mvir (top axis), assuming zL = zS/2.The orange color scheme corresponds to lensed MBBH at zS observable by LISA, while lighter BBHs at zS = 1 observable by ET are in blue.

3 MMFIG. 8 :
FIG.8: Lensed-to-unlensed (SIS) mismatch as a function of the impact parameter.The dashed and dotted lines correspond to modifications in the amplitude and phase in the lensed waveform, respectively.

FIG. 9 :
FIG.9: Prospects of Wave-Optics Features observation by LISA.Lines show the differential optical depth per logarithmic virial mass.We considered equal-mass binaries at different total source-frame masses (colors) and different redshifts (line styles).The green lines show the strong-lensing (i.e.multiple images) probability.

FIG. 10 :
FIG.10: Lensing optical depth, λ, against the redshift of the source for LISA (top) and ET (bottom), for different source masses (colors).Dashed lines show the strong-lensing optical depth (independent of the source mass).

2 GH 3 ×
10 6 M ⊙ (1 + z L ) + y cos θ) , (C7) and where f and g are the auxiliary functions for the Fresnel integrals.Using the conventions of[100], they can be written in terms of the Fresnel sine S and cosine C as Let us apply the time-domain approximation to the SIS, ψ = x.Let us focus on the single-image region, y > 1. Starting from Eq. (23a), it is possible to have a closed-form expression for I (1) :I (1) (τ ) = + y)E (q)]