Natalie B Hogg

Pierre Fleury

Julien Larena

Matteo Martinelli

Measuring line-of-sight shear with Einstein rings: a proof of concept

Keywords: gravitational lensing, strong -cosmology, theory -software, development

Line-of-sight effects in strong gravitational lensing have long been treated as a nuisance. However, it was recently proposed that the line-of-sight shear could be a cosmological observable in its own right, if it is not degenerate with lens model parameters. We firstly demonstrate that the line-of-sight shear can be accurately measured from a simple simulated strong lensing image with percent precision. We then extend our analysis to more complex simulated images and stress test the recovery of the line-of-sight shear when using deficient fitting models, finding that it escapes from degeneracies with lens model parameters, albeit at the expense of the precision. Lastly, we check the validity of the tidal approximation by simulating and fitting an image generated in the presence of many line-of-sight dark matter haloes, finding that an explicit violation of the tidal approximation does not necessarily prevent one from measuring the line-of-sight shear.

INTRODUCTION

Strong gravitational lensing is one of the most elegant probes of the late-time Universe. Massive objects such as galaxies and galaxy clusters bend the local spacetime geometry and produce multiple, distorted images of distant background sources. When the sourcelens-observer alignment is good, the image forms an Einstein ring. In an otherwise perfectly homogeneous and isotropic Universe, the morphology of that ring is a direct probe of the gravitational field generated by the strong lens.

However, on the scales probed by strong gravitational lensing, the Universe is very inhomogeneous. Therefore, we expect that inhomogeneities -typically dark matter haloes -near the line of sight (LOS) will generate weak lensing distortions that will affect the images produced by the strong lensing process. We generically refer to the effects induced by this weak lensing of strong lensing images as LOS effects. These effects have been studied for their impact on, for example, the determination of the Hubble parameter via time-delay cosmography (Gilman et al. 2020), in which case they are considered to be a nuisance effect. But when LOS effects are observable, they also provide some insight into the distribution of dark matter. For example, the effects of individual haloes of dark matter, whether in the main lens or along the LOS, can lead to the detection of dark matter structures on small scales (Vegetti & Koopmans 2009;Vegetti et al. 2012;Despali et al. 2018;Sengül et al. 2022;[START_REF] Nightingale | Scanning For Dark Matter Subhalos in Hubble Space Telescope Imaging of 54 Strong Lenses[END_REF]) and the study of their properties (Vegetti & Koopmans 2009;Vegetti et al. 2014;[START_REF] Ostdiek | [END_REF][START_REF] Zhang | Inferring subhalo effective density slopes from strong lensing observations with neural likelihoodratio estimation[END_REF]. Lineof-sight haloes and main-lens subhaloes lead to different signatures on strong lensing images, and their effects might be disentangled [START_REF] Dhanasingham | Interlopers speak out: Studying the dark universe using small-scale lensing anisotropies[END_REF], although the detectability of individual ★ E-mail: natalie.hogg@ipht.fr

haloes remains an open question (He et al. 2022). Finally, LOS effects in strong lensing may constitute a novel probe of weak lensing in their own right, which is the motivation of the present work.

Line-of-sight perturbations can be modelled in two complementary ways. One option is to add thin perturbers along the LOS in a multiplane lensing formalism (Blandford & Narayan 1986), a methodology which is particularly suitable for taking into account the presence of several subdominant yet clearly identifiable mass concentrations near the LOS. However, this approach is not suitable for modelling the collective effect of a large number of perturbers, due to an excessive number of extra parameters thereby added to the lens model.

This difficulty can be overcome by encapsulating the collective effects of LOS inhomogeneities into an effective, external convergence and shear perturbing the LOS (Kovner 1987;Bar-Kana 1996;Schneider 1997;Birrer et al. 2017;Fleury et al. 2021b). Although this neglects higher-order effects such as flexion, this parametric method is usually considered good enough on large scales. An analysis of both approaches and how to connect them is presented in Fleury et al. (2021a), where one can also find a discussion of higher-order effects.

The idea that the LOS shear can, in principle, be extracted from strong lensing images in galaxy lensing was proposed by Birrer et al. (2017Birrer et al. ( , 2018) ) as an alternative probe of cosmic shear. Based on this approach, a first attempt at correlating this LOS shear with standard weak lensing was made by Kuhn et al. (2021). However, the set-up developed in this programme relied on a non-optimal treatment of LOS effects that overlooked some degeneracies between foreground shear and the main lens modelling. It also employed a somewhat optimistic and simple model for the main lens.

In this paper, we re-analyse the possibility of measuring the LOS shear with Einstein rings generated by galaxy lensing using the minimal model introduced by Fleury et al. (2021a), in lieu of the parametrisation used by Birrer et al. (2017). This model fully accounts for degeneracies between the main lens and its LOS corrections, thereby providing an effective description that depends on only one additional complex parameter. This parameter, denoted 𝛾 LOS , is expected to be mostly independent of the main lens model. In our analysis, we take great care in studying the level of complexity required when modelling the main lens, showing that too simplistic a model can result in biases and a loss of accuracy in the reconstruction of the LOS shear.

We present the first demonstration that the LOS shear can be observed using a series of simulated strong lensing images of increasing complexity. In particular, we show that the LOS shear is generally measurable without systematic biases, and that the next-order LOS effects, while present, do not introduce further biases. These results are obtained thanks to our modifications of the lenstronomy software 1 (Birrer & Amara 2018;Birrer et al. 2021), to which we contributed a new subpackage, LineOfSight, included in the software from version 1.11.0.

The paper is organised as follows: in Section 2 we recapitulate the LOS theory and demonstrate how degeneracies between the LOS shears and other parameters are avoided in the so-called minimal model; in Section 3, we show how the LOS shear is systematically measurable; in Section 4 we explore the validity of the tidal regime. In Section 5 we present our conclusions. In Appendix A we describe our modifications to the strong lensing software lenstronomy, while Appendix B contains details of the parameter choices used to create the mock image data analysed in this work.

Throughout this work, we employ the spatially flat ΛCDM cosmological model described by the Planck 2018 data release (Aghanim et al. 2020), with 𝐻 0 = 67.4 kms -1 Mpc -1 and Ω m = 0.315.

THE MINIMAL LINE-OF-SIGHT MODEL AND ITS ADVANTAGES

In this section, we review the theoretical underpinnings of LOS effects as presented in Fleury et al. (2021a) and demonstrate how the inherent degeneracies can be absorbed in the so-called minimal model, paving the way for the LOS shear to be measured.

Tidal perturbations to a strong lens

As a ray of light passes through the inhomogeneous Universe from source to observer, its path is perturbed by the presence of massive objects. In the simplest strong lensing situation, the deflection of light is assumed to be caused by a single lens. In this case, the whole problem is encapsulated in the lens equation,

𝜷 = 𝜽 -𝜶(𝜽), ( 1
)
where 𝜷 is the angular position of the source, 𝜽 is the angular position of the observed image and 𝜶(𝜽) is the displacement angle, which is in turn determined by the gravitational potential associated with the surface density of the lens, projected along the LOS.

A more elaborate description can be made by including the effects of other nearby objects on the image. The most general approach to deal with this situation is the multi-plane lensing set-up (Blandford & Narayan 1986), in which each additional lens mass is considered to reside in its own plane along the LOS. However, the recursive nature of the multi-plane lens equation makes analytical treatments of complicated lensing situations intractable.

1 https://github.com/lenstronomy/lenstronomy. A simpler approach is to treat the inhomogeneities along the LOS as tidal perturbations which induce a shear and convergence on the image produced by the main lens. This has the advantage of a far simpler theoretical description than the multi-plane case (Kovner 1987;Seitz & Schneider 1994;Bar-Kana 1996;Schneider 1997;Fleury et al. 2021b); see McCully et al. (2014); Fleury et al. (2021a) for explicit derivations of the perturbed lens equation in the tidal regime, starting from the general multi-plane framework.

In this simplified description, the lens equation ( 1) becomes

𝜷 = A os 𝜽 -A ds 𝜶(A od 𝜽), (2)
where the amplification matrices A os , A ds , A od encapsulate all the LOS perturbations. They are commonly parameterised as

A ab = 1 -𝚪 ab , 𝚪 ab =        𝜅 ab + Re (𝛾 ab ) Im (𝛾 ab ) -𝜔 ab Im (𝛾 ab ) + 𝜔 ab 𝜅 ab -Re (𝛾 ab )        , (3)
for ab ∈ {os, ds, od}.The quantities 𝜅 ab , 𝛾 ab , 𝜔 ab represent the convergence, the complex shear, and the rotation of an image as seen from the position (a) with respect to its source at the position (b), due to the inhomogeneities between (a) and (b) -see Fig. 1. For example, if a circular source were placed in the source plane (s), and observed from the main-lens plane (d), then 𝛾 ds would be the complex ellipticity of its image due to the perturbations between (d) and (s). When all the perturbations are set to zero, 𝚪 os = 𝚪 ds = 𝚪 od = 0, we recover the standard single-plane lens equation (1). It is expected that the tidal regime, and hence Equation (2), is valid when all the perturbers besides the main lens are either a smoothly distributed mass on the LOS, which would induce convergence, or compact masses lying far from the LOS, which would induce shear; rotation is caused by successive shears along different axes. Throughout this paper, we shall assume that 𝜅, 𝛾, 𝜔 are all very small, although Equation (2) would remain valid even for large values of those parameters. This description of a perturbed strong lensing system is simple enough to be easily implemented in and rigorously tested by the current state-of-the art strong lensing software, as we will describe in Section 2 and Section 3. It is, in addition, easily extended beyond the tidal regime following the dominant-lens approach of Fleury et al. (2021a), as we will touch on in Section 4.

Absorbing degeneracies in a minimal model

In this work, we aim to show that the LOS shear can be measured from strong lensing images. However, as shown by Fleury et al. (2021a), the three shears involved in Equation (2), 𝛾 os , 𝛾 ds , 𝛾 od , cannot be measured independently, because they are degenerate with the main-lens properties and with one another. This degeneracy is a generalisation of the well-known mass-sheet degeneracy (Falco et al. 1985), and is called the source-position transformation (Schneider & Sluse 2014), since it stems from our ignorance of the exact position and morphology of the source that is lensed. For example, one may multiply Equation (2) with any matrix A and redefine the unknown source position as 𝜷 → A 𝜷; the resulting lens equation is formally identical to Equation ( 2), but with A os → AA os and A ds → AA ds .

Degeneracies may be removed by arbitrarily fixing a definition for the unknown source 𝜷. This is analogous to lifting a gauge freedom by fixing the gauge. As shown in Fleury et al. (2021a), there exists a particularly meaningful choice, leading to a minimal model for the perturbed lens equation. The first step consists in rewriting the displacement angle as the gradient of a potential,

𝜶(𝜽) = d𝜓 d𝜽 , (4)
where

𝜓(𝜽) ≡ 𝐷 ds 𝐷 od 𝐷 os ψ(𝐷 od 𝜽) , (5)
and ψ(𝒙) is twice the projected gravitational potential produced by the surface density of the main lens, Σ(𝒙),

ψ(𝒙) ≡ 4𝐺 𝑐 2 ∫ d 2 𝒚 Σ(𝒚) ln |𝒙 -𝒚| , (6)
where 𝐺 is Newton's constant and 𝑐 the speed of light. Then, multiplying Equation (2) with the combination A od A -1 ds , we obtain the minimal model for the lens equation,

β = A LOS 𝜽 - d𝜓 eff d𝜽 , (7)
with the transformed source position β ≡ A od A -1 ds 𝜷, and

A LOS ≡ A od A -1 ds A os , (8)
𝜓 eff (𝜽) ≡ 𝜓(A od 𝜽) .

Equation ( 7) effectively describes a main lens with potential 𝜓 eff and external tidal perturbations, A LOS , located in the same plane. Since Equations ( 2) and ( 7) are equivalent up to a source position transformation, the corresponding lens models would succeed or fail equally when fitting a strong lensing image. Therefore, any two distinct situations with the same A LOS and 𝜓 eff are observationally indistinguishable. This is why the (od), (os) and (ds) perturbations cannot be independently measured. However, when the perturbations are small,

|𝚪 od |, |𝚪 ds |, |𝚪 os | 1, 𝚪 LOS = 1 -A LOS is linearised as 𝚪 LOS ≈ 𝚪 od + 𝚪 os -𝚪 ds . ( 10
)
This means that the potentially measurable shear component is

𝛾 LOS ≈ 𝛾 od + 𝛾 os -𝛾 ds . (11)
We shall see in Section 2.3 that 𝛾 LOS can indeed be measured from a mock strong lensing image with a much higher precision than any one of 𝛾 od , 𝛾 os or 𝛾 ds taken independently. Another degeneracy revealed by the minimal model ( 7) is between the foreground perturbations (od) and the properties of the main lens, through the effective potential 𝜓 eff of Equation ( 9). Specifically, the 𝜓 eff = cst contours are the images of the 𝜓 = cst contours by the foreground perturbations encoded in A od . This generally makes it hard to distinguish between, for example, the effect of the ellipticity of the main lens and a foreground shear 𝛾 od . Note, however, that in many widely used elliptical lens models, such as the singular isothermal ellipsoid (SIE, Kormann et al. 1994), or its generalisation, the elliptical power law (EPL, Tessore & Metcalf 2015), the foreground shearellipticity degeneracy is only approximate, because in those models the lens ellipticity is implemented at the level of the iso-density contours, not in the iso-potential contours. That is why, in their pioneering numerical experiments, Birrer et al. (2017) could measure both 𝛾 od and 𝛾 os -𝛾 ds = 𝛾 LOS -𝛾 od . Since this only applies to very specific lens models that are not necessarily well motivated from the astrophysical point of view, such an operation should nevertheless be considered very risky; we thus recommend not to try measuring 𝛾 od in general, unless there is a compelling reason to trust the result.

We have chosen to focus the above discussion on the shear. The other components of the tidal matrix 𝚪 LOS are the convergence 𝜅 LOS and the rotation 𝜔 LOS . It is well known that the convergence and source position can be transformed in such a way that the observed lensing image is left invariant -this is precisely the mass-sheet degeneracy -which means that any attempt to measure or constrain the LOS convergence from a strong lensing image would be a pointless undertaking. For this reason, we shall keep the convergence parameters fixed to zero throughout this paper, both when generating mock images and when fitting them.

The rotation is induced by lens-lens coupling, and is generally small when the secondary lenses are individually weak, 𝜔 = O (𝛾2 ). However, it can play a significant role when the shear on the image is large -say, 10% or more. Shears this large would be unusual in real strong lensing data so when generating our mock images we restrict ourselves to smaller values, and keep the rotation parameters fixed to zero. Nevertheless, when fitting images with the minimal model, we include 𝜔 LOS as a free (and in this case nuisance) parameter. Its presence accounts for any slight non-linear effect that may be induced by larger shears, and thus ensures a better recovery of 𝛾 LOS from a given mock image.

A simple example with lenstronomy

To demonstrate the advantage of the minimal model over the full model with a concrete example, we simulate a strong lensing image using the software lenstronomy, which we modified to include the LOS formalism described above. Our modifications are now part of the public version of this code. For full details of our modifications of lenstronomy, see Appendix A.

We create a simple lens using an EPL profile, with additional shears 𝛾 od , 𝛾 os , 𝛾 ds . The ellipticity of the profile, along with its centre and the LOS shear components, are all drawn randomly. We model the source and the lens light with elliptical Sérsic light profiles. The ellipticity of the source is also generated at random. Throughout this work we simulate the noise on our mock images according to the Hubble Space Telescope Wide Field Camera 3 F160W noise settings in lenstronomy (Windhorst et al. 2011). The resulting simulated image is shown in Fig. 2 and the specific model parameters associated with this image are listed in Table 1.

Using the affine-invariant ensemble sampling Markov chain Monte Carlo (MCMC) parameter inference method available in the emcee package 2 (Foreman-Mackey et al. 2013), we then fit the simulated image shown in Fig. 2 with two different models: a model comprised 1).

of an EPL plus the full LOS model (2), i.e. the same model as was used to generate the mock image; and a model comprised of an EPL plus the minimal model ( 7). We sample all the parameters listed in Table 1 except for the lens and lens light central positions, which are kept fixed at the centre of the image. We impose uniform priors on all the parameters and run chains containing ten walkers per parameter. As discussed above, we expect that the individual shear components in the full model will be degenerate with each other, while the 𝛾 LOS combination of the minimal model will be accurately recovered. In both cases (full and minimal model), we expect the foreground shear 𝛾 od to be degenerate with the ellipticity of the main lens.

In Fig. 3, we show the one-and two-dimensional marginalised posterior distributions for the (os) and (ds) shears (left panel) and (od) shear and lens ellipticity (right panel) which result from fitting the simple simulated image with the full LOS model. Note that all the contour plots presented in this paper were produced using the chainconsumer package3 (Hinton 2016).

From these plots, we can see that while the shear parameters are correctly recovered at the level of 1𝜎, there are very strong degeneracies between the (os) and (ds) shears. This results in large uncertainties for 𝛾 os and 𝛾 ds , on the order of 10 %, i.e. comparable or larger than the expected signal. Furthermore, the expected degeneracy between the (od) shear and the lens ellipticity is clearly present. As discussed in Section 2.2, the simplicity of the EPL model artificially alleviates that degeneracy, giving the impression that the main lens's ellipticity can be measured independently of the foreground shear.

In Fig. 4, we show the one-and two-dimensional marginalised posterior distributions for the (od) shear and lens ellipticity (left panel) and LOS shear and lens ellipticity (right panel) which result from fitting the simple simulated image with the minimal model. From these plots, we can see that, as expected, the minimal LOS shear 𝛾 LOS evades any of the previous degeneracies. The input value is perfectly recovered with, in this simple example, a remarkable precision on the order of a few parts in 10 4 .

The results presented in this section demonstrate the clear advantage of using the minimal model over the full model. This advantage can Table 1. The parameters used to simulate the image shown in Fig. 2, along with the priors on the parameters sampled in the MCMC to produce the results shown in Fig. 3

𝑅 Sérsic 0.3 [0 , 10 ] 𝑛 Sérsic 1.0 [0.5, 5] Source light (Sérsic ellipse) 𝑒 1 0.104 [-0.5, 0.5] 𝑒 2 -0.038 [-0.5, 0.5] 𝑥 -0.003 [-0.1 , 0.1 ] 𝑦 0.048 [-0.1 , 0.1 ] 𝛾 od 1 -0.0027 [-0.5, 0.5] 𝛾 od 2 -0.0060 [-0.5, 0.5] LOS 𝛾 os 1 -0.0091 [-0.5, 0.5] 𝛾 os 2 -0.028 [-0.5, 0.5] 𝛾 ds 1 0.019 [-0.5, 0.5] 𝛾 ds 2 0.018 [-0.5, 0.5] 𝛾 LOS 1 -0.031 [-0.5, 0.5] Derived 𝛾 LOS 2 -0.052 [-0.5, 0.5] 𝜔 LOS 0.00033 [-0.5, 0.5]
be understood as being due to a more efficient parameterisation of the problem at hand. It is comparable to the use of 𝑆 8 ∝ 𝜎 8

√ Ω m in standard cosmic shear, which is better constrained than 𝜎 8 (see Hall 2021 for a thorough explanation as to why); or to the use of 𝜔 = Ωℎ 2 in cosmic microwave background constraints. -0 . 0 5 6

- 0 . 2 5 0 . 0 0 0 . 2 5 γ os 2 1 σ 2 σ - 0 . 2 0 . 0 0 . 2 γ ds 1 1 σ 2 σ 1 σ 2 σ - 0 . 2 0 . 0 0 . 2 γ os 1 - 0 . 2 5 0 . 0 0 0 . 2 5 γ ds 2 1 σ 2 σ - 0 . 2 5 0 . 0 0 0 . 2 5 γ os 2 1 σ 2 σ - 0 . 2 0 . 0 0 . 2 γ ds 1 1 σ 2 σ - 0 . 2 5 0 . 0 0 0 . 2 5 γ ds 2 - 0 . 0 4 - 0 . 0 2 0 . 0 0 γ od 2 1σ 2 σ - 0 . 1 0 . 0 0 . 1 0 . 2 e 1,EPL 1 σ 2 σ 1 σ 2σ - 0 . 0 4 - 0 . 0 2 0 . 0 0 γ od 1 - 0 . 1 0 . 0 0 . 1 0 . 2 e 2,EPL 1σ 2 σ - 0 . 0 4 - 0 . 0 2 0 . 0 0 γ od 2 1 σ 2 σ - 0 . 1 0 . 0 0 . 1 0 . 2 e 1,EPL 1 σ 2σ - 0 . 1 0 . 0 0 . 1 0 . 2 e 2,EPL
- 0 . 0 5 2 - 0 . 0 4 8 γ LOS 2 1σ 2 σ - 0 . 1 0 . 0 0 . 1 0 . 2 e 1,EPL 1 σ 2σ 1σ 2σ - 0 . 0 3 6 - 0 . 0 3 2 - 0 . 0 2 8 γ LOS 1 - 0 . 1 0 . 0 0 . 1 0 . 2 e 2,EPL 1σ 2 σ - 0 . 0 5 6 - 0 . 0 5 2 - 0 . 0 4 8 γ LOS 2 1 σ 2σ - 0 . 1 0 . 0 0 . 1 0 . 2 e 1,EPL 1 σ 2 σ - 0 . 1 0 . 0 0 . 1 0 . 2 e 2,EPL - 0 . 0 4 - 0 . 0 2 0 . 0 0 γ od 2 1 σ 2σ - 0 . 1 0 . 0 0 . 1 0 . 2 e 1,EPL 1 σ 2 σ 1 σ 2 σ - 0 . 0 4 - 0 . 0 2 0 . 0 0 0 . 0 2 γ od 1 - 0 . 1 0 . 0 0 . 1 0 . 2 e 2,EPL 1σ 2 σ - 0 . 0 4 - 0 . 0 2 0 . 0 0 γ od 2 1 σ 2 σ - 0 . 1 0 . 0 0 . 1 0 . 2 e 1,EPL 1 σ 2 σ - 0 . 1 0 . 0 0 . 1 0 . 2 e 2,EPL

LINE-OF-SIGHT SHEAR MEASUREMENTS FROM REALISTIC MOCK IMAGES

In this section, we show how the LOS shear is systematically measurable from a larger set of more complicated images than the simple example shown in the previous section. We also assess how the measurement of the LOS shear is affected if the images are fit with lens models that are slightly to drastically different to those used to create the images.

Methodology

We generate a catalogue of 64 mock images, depicted in Fig. 6, following a protocol which we will describe in Section 3.1.1. We (0, 0) halo

( 𝑥 h , 𝑦 h ) baryons ( 𝑥 s , 𝑦 s ) source 𝛽 Figure 5
. Schematic of our composite lens model on the celestial sphere. The baryonic component (orange) is described by a an elliptical Sérsic model whose centre is placed at (0, 0). The dark-matter halo (blue) has an elliptical NFW profile and its centre ( 𝑥 h , 𝑦 h ) is offset with respect to the baryons. This composite main lens is then perturbed by the LOS shears 𝛾 od , 𝛾 os , 𝛾 ds . The source (green) is a perturbed elliptical Sérsic placed at ( 𝑥 s , 𝑦 s ).

then fit those images with four different models which we present in Section 3.1.2, focusing on the inference of the shear 𝛾 LOS . The precision and accuracy of the recovery 𝛾 LOS for each model is assessed using statistical tools we will describe in Section 3.1.3.

Complex mock images from a composite lens model

We increase the complexity of the lens modelling in a number of ways compared to Section 2, so as to mimic real strong lensing images. We now model the lens with two constituent parts: a baryonic component modelled with an elliptical Sérsic profile, and a dark matter halo modelled with an elliptical Navarro-Frenk-White (NFW) profile (Navarro et al. 1996). Note that in both profiles, the ellipticity is implemented at the level of the potential rather than at the level of the convergence. This means that their iso-potential contours are elliptical, but not their iso-density contours. The halo is offset from the centre of the baryons and the ellipticities of each component are not constrained to be of the same magnitude or direction. The set-up is depicted in Fig. 5. We model the lens and source light with elliptical Sérsic profiles. Furthermore, we perturb the source light by adding three additional elliptical Sérsic profiles on top of the original, each with an amplitude of 1 % of the main source. Such perturbations will not be included in the fitting model, and thus aim to account for our general ignorance of the details of the source light. Finally, we shear each image using the full LOS shear model.

We simulate 64 such systems, each with every model parameter drawn at random from physically motivated distributions. The full parameter ranges and distributions are described in detail in Appendix B. The resulting simulated images are shown in Fig. 6, where their broad diversity in morphology can be clearly seen.

We quantify the quality of an image with a notion of cumulated signal-to-noise ratio. Specifically, for an image 𝑖 composed of pixels 𝑝, we define

𝑄 𝑖 ≡ ∑ 𝑝 ∈𝑖 SNR( 𝑝) Θ[SNR( 𝑝) -1] , (12)
where SNR( 𝑝) denotes the signal-to-noise ratio of pixel 𝑝. The lens light is not included in the signal because it is unrelated to the quality of the lensed image, and the noise accounts for both the background and shot noise. The Heaviside function Θ in the sum of Equation ( 12) implies that only the pixels with a signal-to-noise ratio larger than one are accounted for in the computation of the quality 𝑄 𝑖 . The advantage of this definition of quality is that it accounts for both the traditional notion of SNR (the brighter, the better) and for the number of pixels that are covered by the image (the larger and more resolved, the better).

Four fitting models

We employ four models to fit the simulated images in our catalogue, each with a different level of complexity compared to the mocks. We describe them below; see Table 2 for a summary of the parameters included in or excluded from each model.

(i) The comprehensive model is identical to the set-up used to generate the mock images, except that LOS effects are described by the minimal model ( 7) with parameters 𝛾 od , 𝛾 LOS , 𝜔 LOS ; the rotation term 𝜔 LOS is a nuisance parameter included to absorb the small non-linear effects arising from shear-shear couplings. As mentioned earlier, perturbations to the source light are not included.

(ii) The no foreground shear model is identical to the comprehensive model, except that the foreground shear 𝛾 od involved in the effective potential of Equation ( 9) is set to zero, thereby implying 𝜓 eff = 𝜓 here. The motivation for this model is twofold. On the one hand, the presence of 𝛾 od in the minimal model is precisely what makes it different from a more traditional model comprising a composite lens plus external shear, such as those used for some of the H0LiCOW lenses (Wong et al. 2020). It is thus essential to evaluate the role of this new parameter. On the other hand, we expect a degeneracy to exist between the ellipticity of the main lens and the foreground shear. If this degeneracy were perfect, then only two of the three complex parameters 𝛾 od , 𝑒 Sérsic , 𝑒 NFW would suffice to describe an image, allowing us to set one of them to zero. The performance of the no foreground shear model will thus assess the validity of this degeneracy in a composite lens scenario.

(iii) In the aligned halo model, the centres of the baryon and dark matter components are both fixed to (0, 0), thus removing the possibility for the halo to be offset from the baryons. The purpose of this model is to check if neglecting this offset leads to systematic biases in the inference of the other parameters, as was noticed by Gomer & Williams (2021). The aligned model is otherwise identical to the comprehensive one.

(iv) In the power law model, finally, the entire main lens is modelled using a single EPL profile, and with no additional NFW halo, as we used in our simple example in Section 2. This model aims to quantify the impact of neglecting the composite character of a lens on the recovery of the LOS shear.

Evaluating the precision and accuracy of the models

We now present the quantitative criteria that we use to assess the precision and accuracy of the inference of 𝛾 LOS in each model. Firstly, we note that when fitting the images in our mock catalogue, the MCMC parameter inference does not always succeed in obtaining either an upper or lower bound on one or both components of 𝛾 LOS , particularly in the less informative models. We therefore discard the cases where one or both bounds were not obtained from the quantitative comparison of the models that will follow. We refer to the remaining images in the mock catalogue as the productive images, P, denoting their number as 𝑁 P .

Furthermore, for an image 𝑖 ∈ P with a given model, the posterior distribution of 𝛾 LOS is generally non-Gaussian; in particular, the 68%confidence domain is not always symmetric, which translates into asymmetric error bars. We call 𝜎 ± 𝑎,𝑖 the corresponding upper and lower uncertainties at 68%-confidence level on the component 𝑎 ∈ {1, 2} of 𝛾 LOS for the 𝑖 th image. This implies that for each image and model, the uncertainty on the LOS shear is a priori described by four numbers, which is somewhat cumbersome. For further simplicity, we summarise them all by the geometric mean of the symmetrised uncertainty of each component of the LOS shear:

𝜎 𝑖 ≡ 𝜎 + 1,𝑖 + 𝜎 - 1,𝑖 2 𝜎 + 2,𝑖 + 𝜎 - 2,𝑖 2 . ( 13
)
The general precision of a model is then quantified from the mean uncertainty on the LOS shear over all its productive images,

σ ≡ 1 𝑁 P ∑ 𝑖 ∈ P 𝜎 𝑖 . ( 14
)
The accuracy of a model, i.e. its ability to recover the true value of 𝛾 LOS , may be evaluated with a quantity inspired by the reduced 𝜒 2 , accounting for the asymmetry of the uncertainties. Namely, for each component 𝑎 ∈ {1, 2} of the LOS shear, we define

𝜒 2 𝑎 ≡ 1 𝑁 P ∑ 𝑖 ∈ P 1 𝜎 2 𝑎,𝑖 𝛾 LOS,out 𝑎,𝑖 -𝛾 LOS,in 𝑎,𝑖 2 (15) with 𝜎 𝑎,𝑖 ≡ 𝜎 + 𝑎,𝑖 if 𝛾 LOS,out 𝑎,𝑖 ≤ 𝛾 LOS,in 𝑎,𝑖 𝜎 - 𝑎,𝑖 if 𝛾 LOS,out 𝑎,𝑖 > 𝛾 LOS,in 𝑎,𝑖 (16)
where 𝛾 LOS,in 𝑎,𝑖 and 𝛾 LOS,out

𝑎,𝑖

are respectively the true (input) value and the best fit (output) value of the component 𝑎 of the LOS shear of the 𝑖 th image. Again, for further simplicity and to avoid describing the accuracy of a model with two numbers 𝜒 2 1 , 𝜒 2 2 , we shall summarise them by their arithmetic mean,

𝜒 2 ≡ 1 2 𝜒 2 1 + 𝜒 2 2 . (17)
Lastly, we can examine the number 𝑁 >𝑛 𝜎 of outliers, i.e. the images for which the best fit output for 𝛾 LOS lies further than a certain number 𝑛 of 𝜎 from the true input value. Specifically, we consider that an image 𝑖 is an 𝑛𝜎 outlier if

𝛾 LOS,out 𝑖 -𝛾 LOS,in 𝑖 2 > 𝑛 2 𝜎 2 1,𝑖 + 𝜎 2 2,𝑖 , (18)
where 𝜎 𝑎,𝑖 was defined in Equation ( 16). This definition thus accounts for the asymmetry of the error bars.

Results

Examining the fits of a single image

We firstly present the results of the fits to a single image as a representative example of the results in the rest of this section. In Fig. 7 we show the 44 th image in our mock catalogue, along with the one-and two-dimensional marginalised posterior distributions for 𝛾 LOS 1 and 𝛾 LOS 2 obtained when fitting with all four models described above. We chose this particular image to focus on as it is of roughly median quality and the shear parameters were well-sampled in all of the models used to fit it.

From the contour plot, we can see that in both the comprehensive model and the no foreground shear model, the LOS shear is well recovered. On the other hand, in both the aligned model and in the power law model, the posteriors are fairly biased away from the expected values. This indicates that the less informative the model is, the worse the recovery of the shear is. However, the precision of the recovered shears is roughly the same in all four models, implying that the error bar on 𝛾 LOS is insensitive to the model used. This is in line with the results we presented in Section 2, where we showed how 𝛾 LOS is not degenerate with lens model parameters. We will explore this point further in the following.

We now show the result of fitting our entire simulated image catalogue with the four models described above. In Fig. 8

Results of the comprehensive model

In the top two panels of Fig. 8, we show the result of fitting all the images in our catalogue with the comprehensive model. The number of productive images in this case is 𝑁 P = 63, that is 98 % of the complete sample. We can see that in this model, the LOS shear is always very well recovered from the images.

The mean absolute precision of the recovered shear is σ = 0.01 and does not show any correlation with the amplitude of 𝛾 LOS . In comparison, the precision of the inference of 𝛾 LOS found in Section 2 when demonstrating the advantage of the minimal model was about two orders of magnitude better. The better precision obtained in that case can be attributed to the relative simplicity of the lens and source model being used.

Our measure for the accuracy of the LOS shear recovery yields 𝜒 2 = 1.0, which shows that uncertainties are well estimated for both components. As for the number of outliers, we find 𝑁 >2𝜎 = 0 points away from the true value by more than 2𝜎. This is well within the expected proportion of 2𝜎 outliers in the Gaussian case, 5 %. Furthermore, the image-quality criterion 𝑄 defined in Equation ( 12) is well correlated with the size of the error bars on each point, with the most precise constraints generally coming from the best images (the darkest coloured points).

Results of the model with no foreground shear

Next, we fit our catalogue of mock images with the model containing no foreground shear, 𝛾 od = 0. The input-output plot for this case is shown in the second two panels of Fig. 8, with the points plotted in magenta. The number of productive images in this case is one fewer than in the comprehensive model, 𝑁 P = 62, that is 97% of the catalogue.

With the exception of two outliers (which are not visible due to the scale of the plot), the two LOS shear components are generally well recovered in this case too. The precision measure is unchanged compared to the comprehensive model, with σ = 0.01. The accuracy is slightly reduced overall, 𝜒 2 = 3.4, while the number of outliers is increased, with 𝑁 >2𝜎 = 2 (3 % of the productive images), though still in line with the expectation for a normal distribution. This implies that the inclusion of foreground shear in the comprehensive model is not strictly necessary for a better fit to an image. We will discuss this point further below.

Results of the aligned model

Thirdly, we fit our catalogue of images with the aligned halo model. The input-output plot for this case is shown in the third two panels of Fig. 8, with the points plotted in green. We now see how worsening the model can degrade the fit; there is a larger number of outliers and unconstrained points, leading to the number of productive images shrinking to 𝑁 P = 56, that is 88 % of the catalogue. While the precision of the recovery of the two components of 𝛾 LOS is unaffected, σ = 0.01, the accuracy is now very bad, 𝜒 2 = 62.9, with more than two thirds of the productive images being outliers, 𝑁 >2𝜎 = 36.

Results of the power law model

Lastly, we fit our catalogue of images with the power law model. The input-output plot for this case is shown in the bottom two panels of Fig. 8, with the points plotted in blue. Of the four models used to fit the simulated images, this is clearly the worst, with a large number of unconstrained points that result from the difficulty the MCMC has when trying to fit complex images with an overly simplistic model: 𝑁 P = 28. Of the remaining points, the accuracy measure is 𝜒 2 = 92.2, with again over two thirds of the productive images leading to a measurement of 𝛾 LOS that is biased by more than 2𝜎, 𝑁 >2𝜎 = 20. The mean precision of the recovered shear components remains in line with that of the other models, σ = 0.01.

Discussion

From these results, summarised in Table 2, we can clearly see the importance of using a sufficiently feature-rich model to fit complicated images, if a reliable and precise measurement of the LOS shear is desired, with well-sampled posterior distributions that yield upper and lower bounds on the recovered shear parameters. This is emphasised in Fig. 9, where we plot the histogram of the mean error on the shear in each model (left panel) and the histogram of the difference between the input and output values of the shear (right panel). From this plot, we can see how the distribution of the size of the error is extremely similar across all the models. However, in the right-hand panel, we can see how the comprehensive model yields consistently more accurate shear measurements than the rest, with the greatest difference between input and output shear being less than 2𝜎; this is in stark contrast to the aligned halo and power law models, where MNRAS 000, 1-19 (2023) the greatest difference grows as large as 20𝜎. The accuracy of the no foreground shear model is very close to that of the comprehensive model, barring the two outliers previously mentioned above.

The MCMC sampling is also less successful if a poor model is used: when using the no foreground shear model, we found two-sided bounds on both components the LOS shear for nearly 100 % of the catalogue, whereas when using the power law model, we found twosided bounds on both the shear parameters for less than half of the catalogue.

We can further examine the correlation between the quality measure 𝑄 𝑖 defined in Equation ( 12) and the uncertainty 𝜎 𝑖 on 𝛾 LOS defined in Equation ( 13) for each image 𝑖, as plotted in the left-hand panel of Fig. 10. From this plot, we can see a strong negative correlation between these quantities, with larger average errors coming from images of worse quality (smaller 𝑄). A simple linear fit yields log 10 𝜎 ≈ -0.76 log 10 𝑄 + 1.01. This implies that for precise measurements of 𝛾 LOS to be obtained, poor quality images should be removed from the catalogue. On the other hand, in the right-hand panel of Fig. 10, we plot the difference between the input and output 𝛾 LOS in units of 𝜎 against 𝑄. We can see that there is virtually no correlation between these quantities; in other words, the accuracy of the shear recovery is not affected by the image quality. Note that we chose to plot these relationships for the comprehensive model as an example, but the trend remains the same for the no foreground shear model.

The excellent recovery of 𝛾 LOS obtained from fitting with the no foreground shear model show that the degeneracy between the foreground shear and the properties of the main lens is indeed strong enough that the composite lens images can be satisfactorily obtained by the MCMC parameter inference even when 𝛾 od is fixed to zero. This is another powerful argument in favour of the use of the minimal model and again strongly emphasises the fact that 𝛾 LOS is the only notion of shear that can evade such degeneracies and thus be accurately measured from strong lensing images.

For test purposes, we tried to increase the amplitude of the source perturbations, which we recall are not accounted for in the fitting models, so as to mimic a poorer modelling of the source light. In that case, we found that the reduced 𝜒 2 , which quantifies the accuracy of the measurement of 𝛾 LOS , increases and becomes larger than unity even for the comprehensive model. This suggests that a too simple modelling of the lens light could lead to systematic biases in the measurements of the LOS shear. Such a systematic effect could be mitigated by using complex source models when fitting the images, but this approach may not translate easily when it comes to analysing real data (Nightingale & Dye 2015;Nightingale et al. 2018).

The results presented in this section show that the LOS shear is systematically measurable from simulated images if the lens model used to fit the images is of comparable complexity to that used to generate them. They also show that the foreground shear does not need to be included for good inference if the model for the main lens is sufficiently rich.

VALIDITY OF THE TIDAL APPROXIMATION

The analysis and results presented so far were based on the lens equation ( 2), which supposes that all LOS perturbations to a strong lens -typically dark matter haloes that would lie near the optical axis -can be modelled as tidal fields, whose effect is encoded in the amplification matrices A od , A os , A ds . In this section, we assess the validity of this tidal approximation by fully simulating the effect of randomly distributed haloes on an image.

The tidal regime and beyond

A LOS perturber is said to be in the tidal regime if the deflection field that it produces can be described by a quadratic potential across the strong lensing image under consideration. Suppose for instance that the observer sees a source that is lensed by a halo ℎ in the absence of the main lens and any other perturber. The corresponding Fermat potential can be expanded at second order in the position 𝜽 in the vicinity of the LOS (𝜽 = 0),

𝜓 ℎ (𝜽) = 𝜓 ℎ (0) + 𝜶 ℎ (0) • 𝜽 + 1 2 𝜽 • A ℎ 𝜽 + . . . ( 19
)
where 𝜶 ℎ (0) is the displacement angle caused by ℎ for an image observed along the optical axis, and A ℎ is the contribution of ℎ to the amplification matrix A os . A similar expansion could be performed for the (od) and (ds) cases; see Fleury et al. (2021a) for further details on the general theory of LOS perturbations in strong lensing. We have chosen to drop the (os) subscripts or superscripts when there is no ambiguity to alleviate notation. The tidal regime for the halo ℎ is satisfied as long as Equation ( 19), together with its (od) and (ds) analogues, remain valid up to an angle on the order of the main lens's Einstein radius, 𝜃 E . This requirement is satisfied when ℎ lies sufficiently far from the LOS; as it gets closer, the higher-order terms hidden in the suspension points of Equation ( 19) may become non-negligible. Among them, the first one that we may worry about is the cubic term, O (𝜃 3 ), generally referred to as flexion because its effect is to produce arc-shaped images (Bacon et al. 2006). Flexion can be understood as gradients of A, i.e. gradients of convergence and shear; it consists of two complex numbers4

F ≡ 𝜕𝛾 = 𝜕 * 𝜅 , G ≡ 𝜕 * 𝛾 , (20)
where 𝜕 ≡ (𝜕/𝜕𝜃 1 -i𝜕/𝜕𝜃 2 )/2 is the standard complex derivative with respect to the complex counterpart 𝜃 1 + i𝜃 2 of 𝜽 = (𝜃 1 , 𝜃 2 ).

When the flexion due to LOS perturbers is properly accounted for in the lens equation, Equation ( 2) is significantly more complicated, with the addition of four extra complex parameters (Fleury et al. 2021a).

A line of sight populated with haloes

Since all the previous results we presented were obtained under the assumption that the tidal approximation is valid, we need to understand if it is a good description of the real Universe: how often do perturbers escape this regime, and do they jeopardise our ability to measure the LOS shear from strong lensing images?

In order to test this, we simulate a simple strong lens with a very large number of dark matter haloes in a volume around and along the LOS using the multi-plane lensing formalism. We compute the shear induced on the image by these LOS haloes and then fit the image using the same MCMC procedure followed in the previous sections, with a model again based on the minimal model. If the tidal approximation is valid, we expect the model to provide a good fit and an accurate recovery of the induced shear from the image. If the tidal approximation is significantly broken, we expect to see some sign of it, which could be a biased recovery of the shear, or a bad fit due to the non-negligible role of higher-order distortions.

The main lens and source light are modelled in the same way as in Section 2, except that, for the sake of simplicity, we do not include any lens light. The ellipticities of the lens and source are drawn at random, while the remaining parameters are chosen in a physically 13)] in the comprehensive model. The grey dashed line represents a linear fit, log 10 𝜎 ≈ -0.76 log 10 𝑄 + 1.01, highlighting the correlation between these quantities. Right panel: quality criterion 𝑄 versus the residuals of 𝛾 LOS in units of 𝜎 in the no foreground shear model. There is no significant correlation between these quantities. realistic way. Specifically, we set the Einstein radius of the main lens 𝜃 E = 1.5 , the slope of the elliptical power law potential 𝛾 EPL = 2.6, the apparent magnitude of the source as 23, the half-light radius of the source 𝑅 Sérsic = 0.3 and the Sérsic index 𝑛 Sérsic = 6.0. The main lens is placed at 𝑧 d = 0.5 and the source is placed at 𝑧 s = 1.5. The angular positions of the centre of the main lens and of the source are chosen so as to be observed at 𝜽 = 0 in the presence of the haloes, 𝜷 s = -𝜶 os (0) and 𝜷 d = -𝜶 od (0). Conversions between distances and redshifts are performed with astropy (Astropy Collaboration et al. 2022). As a reminder, we work in the spatially flat ΛCDM cosmology described by the Planck 2018 data (Aghanim et al. 2020), with 𝐻 0 = 67.4 km s -1 Mpc -1 and Ω m = 0.315. We generate the halo population as follows. We consider a region of space consisting of a comoving cylinder around the optical axis, with comoving radius 𝑅 = 11.6 Mpc. This number is chosen so that a point lens with mass 𝑀 = 10 12 𝑀 located halfway between the observer and the source and with a transverse comoving distance 𝑅 from the optical axis would produce an (os) shear of 10 -6 on the latter. We then randomly populate this region with a uniform comoving number density

𝑛 = ∫ 𝑀 max 𝑀 min d𝑀 𝑀 d𝑛 d ln 𝑀 , (21)
where d𝑛/d ln 𝑀 denotes the halo mass function (HMF) evaluated following Tinker et al. (2008) as implemented in the colossus package5 [START_REF] Diemer | [END_REF] at 𝑧 d = 0.5. We fix a lower mass cutoff to 𝑀 min = 10 8 𝑀 corresponding to the typical mass of galactic subhaloes in the cold dark matter scenario, and 𝑀 max = 10 15 𝑀 , which yields 𝑛 = 5.6 Mpc -3 . The mass of each halo is randomly drawn from the same HMF, shown in Fig. 11. This procedure results in a very large number of haloes, 𝑁 tot = 1.1 × 10 7 in the cylindrical region, among which many are practically too light or too far from the LOS to be relevant to the analysis. For computational efficiency, we thus remove the haloes for which a point lens with the same mass would individually produce an (os) shear with magnitude below 10 -6 on the optical axis. After this operation, we are left with 𝑁 = 3.4 × 10 4 haloes , whose positions are depicted in Fig. 12. The points in red indicate the haloes that significantly evade the tidal regime (see Section 4.3).

The haloes themselves are modelled with NFW profiles. We compute their concentration using the Diemer & Joyce (2019) model for the mass-concentration relation, which is available in colossus. We then use the cosmology subpackage of lenstronomy to convert these into the angular scale radius and the deflection angle at the scale radius which constitute the remaining parameters of the NFW mass profile in lenstronomy. left and middle panels of Figure 13 show the distribution of individual (os) shear and convergence produced by the 𝑁 haloes on the LOS. We can see that the individual tidal effects are generally small, except for a handful of haloes producing convergences and shears of percent order.

Lastly, in each plane containing a halo, we add a uniform component of negative density in order to compensate for the added mass under the form of haloes. This procedure mostly aims to avoid producing a too strong net convergence along the LOS. In practice, we set the negative convergence in each of those planes ℎ so as to exactly compensate the positive convergence 𝜅 ℎ (0) produced by the halo on the optical axis. This choice is not strictly consistent from a physical point of view, but it has the virtue of simplicity.

Is the tidal regime a good approximation?

A first, qualitative, assessment of the tidal approximation can be made by examining the variations of the convergence and shear fields due to the haloes only, i.e. in the absence of the main lens. This is shown in Fig. 14, where we observe that the (os) convergence and shear due to the haloes display variations on the order of 100 % on the scale of the main lens's Einstein radius. The convergence map in particular (left panel of Fig. 14) reveals the presence of a handful of haloes sitting very close to the optical axis, hence strongly evading the tidal regime, while remaining subcritical.

This diagnostic can be further quantified by computing the relative variation Δ of the (os) shear 𝛾 as we move away from the optical axis by one unit of the main lens's Einstein radius, 𝜃 E = 1 ,

Δ 2 = 𝛾(𝜃 E x) -𝛾(0) 𝛾(0) 2 + 𝛾(𝜃 E ŷ) -𝛾(0) 𝛾(0) 2 , (22)
where ( x, ŷ) is an orthonormal basis of the image plane. Note that if the gradients of the shear can be considered constant across the image (flexion regime), then Δ can be expressed in terms of the two complex flexions defined in Equation ( 20) as

Δ 2 = 1 2 𝜃 2 E |𝛾| -2 |F | 2 + |G| 2 . ( 23
)
We first consider the value Δ ℎ of Δ for each halo ℎ taken individually. The distribution of Δ ℎ for the 𝑁 haloes is shown in the right panel of Fig. 13. We can see that a considerable number of haloes produce a shear with significant variations on the scale of 𝜃 E . Figure 12 shows in red the haloes for which Δ ℎ ≥ 0.1; they are unsurprisingly clustered around the optical axis. We find that about 14 % of the total sample (4.8 × 10 3 haloes) have these especially large values of Δ ℎ .

As expected from the above consideration, the combined effect of all the haloes does escape the tidal regime. We indeed find Δ tot = 0.44 for the total (os) shear, whose central value is 𝛾 os = -2.4 × 10 -2 + 3.8i × 10 -3 . This number drops to Δ tot = 0.045 if we remove all -2 the haloes ℎ with Δ ℎ ≥ 0.1, in which case the cumulative shear becomes 𝛾 os = -1.8 × 10 -2 -9.8i × 10 -4 . This indicates that the haloes escaping the tidal regime are not negligible in the total shear produced along the LOS.

0 2 θ x [arcsec] -3 -2 -1 0 1 2 3 θ y [arcsec] -2 0 2 θ x [arcsec] γ 1 -2 0 2 θ x [arcsec]
In order to quantify the impact of beyond-tidal effects on the measurement of the LOS shear, we generate two images: the first one is the image that would be observed from a system consisting of our main lens plus 𝑁 haloes and the second one the image that would be observed from the same system but without the haloes ℎ such that Δ ℎ ≥ 0.1. We do this using the multi-plane formalism in lenstronomy. The two images are shown in the top and bottom left-hand panels in Fig. 15.

We then fit these two images using the minimal model, plus an EPL mass profile for the main lens. For this part of the work, we changed our MCMC parameter inference method to use ensemble slice sampling (Karamanis & Beutler 2021) as implemented in the zeus package 6 (Karamanis et al. 2021), which we again modified lenstronomy to include. We made this change due to the increased difficulty of the inference which required a faster and more powerful sampling method than the affine-invariant ensemble sampler provided in emcee.

In Fig. 16, we show the one-and two-dimensional marginalised posterior distributions of the (od) and LOS shears that result from fitting the images in Fig. 15. The dashed lines indicate expected values of 𝛾 od (0) and 𝛾 LOS (0) = 𝛾 os (0) + 𝛾 od (0) -𝛾 ds (0) exactly along the optical axis, in the absence of the main lens. We firstly consider the 6 https://github.com/minaskar/zeus. left-hand panel of this figure, which shows the results from fitting the image produced in the presence of the full population of haloes. The expected LOS shear is recovered with good accuracy, but with a precision degraded by about two orders of magnitude compared to the results of Section 2, whose image was generated from a strictly tidal LOS perturbation model. Furthermore, as we can see from the top right panel of Fig. 15, there is an obvious strong residual between the observed image and the reconstructed image. The dipolar morphology of these residuals is reminiscent of the F -type flexion signal predicted by Fleury et al. (2021a) (see Fig. 7 therein). Further investigation is nevertheless required in order to assess the detectability of such a signal in a realistic set-up.

Taking now the contour plot in the right-hand panel of Fig. 16, which results from fitting the image with the strongly non-tidal haloes removed, we can see that once again the LOS shear is recovered with good accuracy, while the precision is not significantly improved. The residuals exhibited in the case with the full halo population are still present with a similar amplitude, as can be seen in the bottom right panel of Fig. 15.

These results indicate that, quite surprisingly, significant violations of the tidal regime do not entirely prevent one from measuring a notion of LOS shear from a strong lensing image. The tidal model is thus extremely resilient. Furthermore, the best-fit LOS shear appears to be unbiased with respect to the expected values on the optical axis. 7Beyond-tidal effects then significantly increase the uncertainty on the shear measurement, which may be attributed to the fact that the shears are simply not single-valued across the image. They also produce visible patterns in the residuals of the best fit. A natural question is whether those patterns can be attributed to a combination of flexion signals; we shall address this point in future work.

We stress that the quantitative conclusions of this section should be taken with a pinch of salt. Indeed, in the scenario we constructed, haloes were initially randomly distributed in space without any clustering. However, in reality, halo positions are strongly correlated, with many small haloes found in the vicinity of larger ones, and following the cosmic two-point correlation function on large scales. Such correlations imply the presence of large coherent structures far from the LOS that could produce a rather uniform shear. In other words, we expect the role of beyond-tidal effects to be overestimated in the present analysis. An important next step for future work will be to simulate a distribution of haloes that follows a cosmologically motivated two-point correlation function, possibly by obtaining a distribution of haloes and their physical characteristics such as mass and concentration from an 𝑁-body simulation.

CONCLUSIONS

In this paper, we have provided a clear demonstration that the minimal model, based on the dominant lens approximation and the tidal regime, to describe line-of-sight effects in strong gravitational lensing, evades the degeneracies inherent in other formalisms which try to address the same problem. Armed with this powerful model, we showed how the line-of-sight shear can be measured with percentlevel accuracy from complicated strong lensing images even when source perturbations are not modelled. Furthermore, we tested the validity of the tidal approximation by simulating the effect of around 3.4 × 10 4 dark matter haloes on a strong lensing image, finding that while the approximation may not be completely robust, LOS shear measurements are nevertheless possible.

Concretely, in Section 2, we laid out the theoretical background to this work, beginning with the treatment of LOS effects in the tidal regime. We argued that the degeneracies present between the three shear terms in the lens equation in this regime means that they cannot be measured independently of each other or of the model parameters describing the main lens, but that a so-called minimal model can be constructed in which these degeneracies are eliminated. Using lenstronomy -modified for this work to include this LOS formalism -we demonstrated the advantage that the minimal model has over the full model, by constructing a simple lens model with shear and fitting the resulting image with both the full and minimal model. As expected, the aforementioned degeneracies were present in the full model but not in the minimal model.

Equipped with this demonstration of the advantage of the minimal model, we applied it to the inference of the LOS shear from a set of 64 simulated strong lensing images, each produced by a rich composite lens model consisting of an elliptical dark matter halo offset from an elliptical baryonic core, along with LOS shear and perturbations to the source light, in Section 3. We fitted these images with a series of models, among which the optimal performance was obtained in the case where all the features of the simulated lenses are included save the foreground shear. This model is very similar to the composite lens models used to fit the H0LiCOW lens SDSS 1206+4332 by Birrer et al. (2019) and the STRIDES lens DES J0408-5354 by Shajib et al. (2020), both re-analysed by the TDCOSMO collaboration, again with a composite lens model (Millon et al. 2020), and is similar to the multi-baryonic component lens models used by [START_REF] Williams | Two-component mass models of the lensing galaxy in the quadruply imaged supernova[END_REF] and Nightingale et al. (2019), for example. In this case, we concluded that the LOS shear could be measured with an average absolute uncertainty of 1 %, without any systematic bias or degeneracy with the lens model parameters. We also showed how the inference of the LOS shear worsens with decreasing model complexity, as is to be expected.

Finally, in Section 4, we explored the validity of an important assumption: that LOS perturbers can be treated in the tidal regime. By simulating a scenario in which a volume around the LOS is populated with dark matter haloes, we analysed the images produced in the case where potentially beyond-tidal haloes were present and in the case where these haloes were removed. We observed a clear breakdown of the tidal regime, albeit with no effect on the accuracy of the measurement of the LOS shear. However, the precision was significantly worsened with respect to the best-case scenario of a very simple lens model in Section 2. As a result of the way we modelled the distribution of the dark matter haloes, it is likely that we have overestimated their beyond-tidal effects, making this a very conservative estimate of the achievable precision of the LOS shear.

Lastly, we would like to highlight the important debate in the literature regarding whether external shear is truly external or notin other words, if the distortions to a strong lensing image which have been traditionally attributed to the presence of perturbers external to the main lens are in fact due to another, perhaps poorly understood or controlled-for effect. One example of this was brought to light by Van de Vyvere et al. (2020), who showed how incorrect truncation at the edges of a mass profile for a main lens can induce a spurious external shear. While this effect is important for reconstructions of lensing potentials from mass maps, it does not play a role in our current work due to our use of exact mass profiles which essentially extend to infinity. In that work, the spurious shear was also found to decrease to negligible magnitudes provided a large enough region for the mass map was used. The existence of external shear was recently further called into question by Etherington et al. (2023), who showed how the best-fit value of an external shear parameter does not faithfully represent the actual shear if the lens mass is modelled by a simplistic elliptical power law. This finding is completely in agreement with the results of the current work.

In conclusion, the prospects for the use of the LOS shear as a new cosmological probe are excellent -provided that the main lens mass is modelled correctly. The LOS shear, which as we have shown is free from the usual degeneracies that plague typical strong lensing observables, will be able to be measured not only in current strong lensing data but from the many thousands of new strong lenses which are expected to be observed in the next five years by LSST and JWST, and further afield, by Euclid. The LOS shear thus has the potential to be employed in synergy with weak lensing surveys, magnifying the information that we receive from the inhomogeneous sky and boosting the constraints on cosmological parameters that are already in our hands. density in the lens plane in units of 𝑀 /arcsec 2 and the normalisation 𝐼 is the integrated convergence for the given lens when 𝜅 eff = 1 in units of arcsec 2 . Finally, the centre of each lens is fixed to zero, 𝑥 = 𝑦 = 0, i.e. exactly aligned with the optical axis.

B2 Lens light: Sérsic ellipse

We model the lens light with an identical Sérsic ellipse profile so that the lens light traces the lens mass exactly. We use the SERSIC_ELLIPSE light profile in lenstronomy, which takes the same parameters as the SERSIC_ELLIPSE_POTENTIAL mass profile, less the effective convergence and with the addition of an apparent magnitude.

We start by computing the absolute magnitude of each lens,

𝑀 = 𝑀 -2.5 log 10 M Υ , (B5)
where 𝑀 = 4.74 is the absolute magnitude of the Sun and Υ = 2 is the mass-to-light ratio for the lens galaxy. The apparent magnitude is then

𝑚 = 𝑀 + 5 log 10 𝐷 1 Mpc + 25, ( B6
)
where 𝐷 is the luminosity distance to the lens, computed from the redshift and angular diameter distance,

𝐷 = (1 + 𝑧) 2 𝐷 od . (B7)
The remainder of the lens light parameters for each lens are the same as those used for the baryonic core.

B3 Dark matter halo: NFW ellipse

We model the dark matter halo using the NFW_ELLIPSE mass profile in lenstronomy, which requires the scale radius, 𝑅 s , the deflection angle at the scale radius 𝛼, the ellipticity components 𝑒 1 , 𝑒 2 and the centre position, 𝑥, 𝑦.

We start by drawing the ratio of dark matter to baryons in the lens galaxy from a uniform distribution, 𝜆 ∼ U (0.3, 1.0), and thus compute the mass of each halo via

M halo = M/𝜆. (B8)
We draw the concentration of each halo randomly from a log-normal distribution, 𝑐 ∼ ln N (𝑐 mean , 0.2), where 𝑐 mean = 10 1.1 . We then use lenstronomy to compute the scale radius and deflection angle for each halo from the mass and concentration. We follow the same procedure as for the ellipticity of the baryonic component to compute the ellipticity of each NFW halo. Lastly, we displace the centre of each halo from the baryonic core by drawing an offset (in parsecs) from a normal distribution, 𝑜 ∼ N (0, 300). This distribution is motivated by the findings of Kuhlen et al. (2013), who found a 300-400 pc offset in a simulation of a Milky Way-like galaxy. An even larger offset of 1.62 kpc has been observed in a central galaxy in the Abell 3827 cluster (Massey et al. 2015).

B4 Source light: Sérsic ellipse

As with the lens light, we model the source component using the SERSIC_ELLIPSE light profile in lenstronomy, and our methodology for computing the model parameters for each source is the same as followed for the lens light, except for two differences: the half-light radius and apparent magnitude are computed using the angular diameter distance from observer to source, 𝐷 os , rather than observer to lens, 𝐷 od ; and the centre of the source, 𝑥, 𝑦, is allowed to be offset from the optical axis. We draw the source centre randomly from a uniform distribution on a disk, with the maximum allowed displacement from the optical axis being given by whichever is smallest: the half-light radius or the effective Einstein radius of the lens.

B5 Einstein radius

We estimate the Einstein radius of our composite lenses by neglecting the ellipticity and offset of their components, i.e. as if the lenses were axially symmetric; 𝜃 E is, thus, the solution of 𝜃 E -𝛼 Sérsic (𝜃 E ) -𝛼 NFW (𝜃 E ) = 0 , (B9)

where 𝛼 Sérsic and 𝛼 NFW are the displacement angles individually produced by the baryon and dark-matter component of the lens. After drawing randomly their parameters, we compute 𝜃 E by solving Equation (B9) numerically, and we only keep the lenses such that 𝜃 E > 0.5 , to ensure that the images we produce have arcs or rings which are reasonably well separated from the lens light.

B6 LOS shear

We add the LOS shear to the images using the LOS model that we implemented in lenstronomy. The two components of the od, os and ds shears, 𝛾 1 and 𝛾 2 , are drawn from uniform distributions on a disk. This is done by drawing the square of each shear component randomly from a uniform distribution, 𝛾 2 ∼ U (0, 𝛾 2 max ), where we set 𝛾 max = 0.025 to ensure that the LOS shear, 𝛾 LOS , is never larger than 5% -as we discussed in Section 2, large shears can induce rotation in the image with respect to the source.

The two components of each shear are then computed via 𝛾 1 = 𝛾 cos(2𝜙), (B10)

𝛾 2 = 𝛾 sin(2𝜙), ( B11
)
where 𝜙 is drawn randomly from a uniform distribution, 𝜙 ∼ U (0, 𝜋).

Figure 1 .

 1 Figure 1. When LOS perturbations to a strong lens can be treated in the tidal regime, they intervene in the lens equation (2) via three distortion matrices: A od , A os , A ds . For each couple ab, A ab represents the distortions to a small source placed at (b) as observed from (a) due to the perturbers only.

Figure 2 .

 2 Figure 2. The image produced by a simple lens model consisting of an elliptical power law mass profile and a elliptical Sérsic lens light profile, with an elliptical Sérsic source (model parameters listed in Table1).

  and Fig.4. The ellipticities, positions and shear parameters were all drawn at random. The lens light position and ellipticity are fixed to be the same as those of the EPL profile. The magnitudes quoted are apparent magnitudes. The expected values of the minimal LOS quantities 𝛾 LOS , 𝜔 LOS are derived from Equation (8) using the decomposition of Equation (3

Figure 3 .

 3 Figure 3. Inference with the full LOS model (2), from fitting the simulated image in Fig. 2. Left panel: one-and two-dimensional marginalised posterior distributions of the (os) and (ds) shears. Right panel: Same for the (od) shear and lens ellipticity. The dashed lines show the input value of each parameter.

Figure 4 .

 4 Figure 4. Inference with the minimal model (7), from fitting the simulated image in Fig. 2 with the minimal model. Left panel: one-and two-dimensional marginalised posterior distributions of the effective LOS shear and lens ellipticity. Right panel: Same for the (od) shear and lens ellipticity. The dashed lines show the input value of each parameter.

Figure 6 .

 6 Figure 6. Our catalogue of 64 simulated images produced according to the protocol presented in Section 3.1.1.

  , we show the input LOS shear components as computed from the three individual shears of the full model versus the output LOS shear components as recovered from the images by the MCMC parameter inference. The colour bars show the log of the quality criterion 𝑄. Darker colours thus represent larger values of 𝑄 and hence larger SNR values. Only the results from the productive images in each model are shown in this figure.

Figure 7 .

 7 Figure 7. Left panel: the 44 th image in our catalogue. Right panel: the one-and two-dimensional marginalised posterior distributions for 𝛾 LOS 1 and 𝛾 LOS 2 obtained by fitting this image with the power law model (blue), aligned halo model (green), no foreground shear model (magenta) and comprehensive model (orange). The dashed lines represent the expected values.

Figure 8 .

 8 Figure 8. Input 𝛾 LOS,in versus output 𝛾 LOS,out in the comprehensive model (orange), no foreground shear model (magenta), aligned halo model (green) and power law model (blue). The left-hand column shows the results for 𝛾 LOS 1 and the right-hand column shows 𝛾 LOS 2 . The dashed line in each panel is where 𝛾 LOS, in = 𝛾 LOS, out . The brightness of each point indicates the log of the quality criterion 𝑄 for the associated system. A larger 𝑄 (darker colour) corresponds to a larger SNR, or in other words, a better quality image.

Figure 9 .

 9 Figure 9. Left panel: histogram of the uncertainty 𝜎 on 𝛾 LOS , as defined in Equation (13) for the four models. Right panel: histogram of the difference between the input and best-fit output values of 𝛾 LOS in units of 𝜎.

Figure 10 .

 10 Figure 10. Left panel: quality criterion 𝑄 versus the estimated uncertainty 𝜎 on 𝛾 LOS [Equation (13)] in the comprehensive model. The grey dashed line represents a linear fit, log 10 𝜎 ≈ -0.76 log 10 𝑄 + 1.01, highlighting the correlation between these quantities. Right panel: quality criterion 𝑄 versus the residuals of 𝛾 LOS in units of 𝜎 in the no foreground shear model. There is no significant correlation between these quantities.

Figure 11 .

 11 Figure 11. The Tinker et al. (2008) halo mass function between 10 8 𝑀 and 10 15 𝑀 .

Figure 12 .

 12 Figure 12. The cylindrical region populated with our sample of 𝑁 haloes, shown in projection in blue (for clarity only a tenth of the population is plotted), of which around 14 % have a relative variation of shear Δ ℎ greater than 0.1 across 𝜃 E (of which we again plot a tenth), shown in red.

Figure 13 .

 13 Figure 13. Left panel: the histogram of the convergence 𝜅 ℎ produced by individual haloes ℎ on the optical axis. Middle panel: same for the modulus of the shear |𝛾 ℎ |. Right panel: same for the relative variation of shear Δ ℎ on the angular scale 𝜃 E . These are all the (os) quantities.

Figure 14 .

 14 Figure 14. Convergence field 𝜅 os (𝜽) (left panel) and shear field 𝛾 os (𝜽) (middle and right panels) produced by the placed haloes around the LOS. Dashed circles indicate the Einstein radius of the main lens.

Figure 15 .Figure 16 .

 1516 Figure 15. The image produced in the presence of the full population of haloes (top left) along with the reconstruction from the best-fit parameters found by the MCMC (top middle) and the normalised residuals of the difference between the two (top right). The bottom panels show the same but without the haloes ℎ escaping the tidal regime (Δ ℎ ≥ 0.1).

Table 2 .

 2 Summary of the properties and results of the four models used to fit the mock images of Fig. 6. Comprehensive model No foreground shear model Aligned model Power law model

	Features
	Foreground shear
	LOS shear
	Baryon ellipticity

Table A1 .

 A1 Example output of the SinglePlaneLOS and MultiPlane classes for comparison.

		SinglePlaneLOS	MultiPlane			Difference
	Displacement angle	(0.92460116, 0.77805976)	(0.92460116, 0.77805976)			4.44089210 × 10 -16 , 3.33066907 × 10 -16
	Hessian Time delay	      -0.01004559 0.15331121 0.11183298 -0.0151521       6.19293093 days	      -0.01004559 0.15331119 0.11183299 -0.0151521 6.19293093 days	     	      	6.50405486 × 10 -9 -1.94965391 × 10 -10 -1.18127730 × 10 -11 days -3.40569642 × 10 -10 -2.69573809 × 10 -8	      

MNRAS 000,1-19 (2023)

https://github.com/dfm/emcee.MNRAS 000,1-19 (2023)

https://github.com/Samreay/ChainConsumer.

These definitions follow the conventions of Fleury et al. (2021a), which differ by a factor two from those of Bacon et al. (2006).MNRAS 000,1-19 (2023)

https://bitbucket.org/bdiemer/colossus.

This conclusion regarding the impact of higher-order LOS perturbations can be compared with the recent finding that the boxiness or diskiness of strong lenses do not significantly bias 𝐻 0 measurements (Van de Vyvere et al. 2022).MNRAS 000,1-19 (2023)

ACKNOWLEDGEMENTS

We thank all the participants of the LOS workshop held in Montpellier in June 2022 for the lively and enriching discussions about this work; we are particularly grateful to Simon Birrer for detailed discussions about lenstronomy and to Daniel Johnson for the implementation of LOS effects in the Fermat potential which allowed for the time delay comparison in Appendix A. We also thank Daniel Gilman for reviewing our lenstronomy pull request, Anowar Shajib for providing information about TDCOSMO lens modelling, Sherry

Suyu for suggesting that we consider a notion of signal-to-noise ratio to quantify the quality of our images, and Dominique Sluse and Liliya Williams for their detailed comments on the paper. We are also grateful to the anonymous referee for their comments which led to a substantial improvement of our results. Some of the colourmaps used in this work were provided by the CMasher package 8 (Van der Velden 2020); other colour palettes were taken from colorbrewer2. Some of the numerical computations for this work were done on the Sciama High Performance Compute (HPC) cluster, which is supported by the ICG, SEPNet and the University of Portsmouth. NBH is supported by a postdoctoral position funded by the French Commissariat à l'énergie atomique et aux énergies alternatives (CEA). JL's work was supported by a starting grant from the University of Montpellier and the French Ministry of Higher Education, Research and Innovation. MM acknowledges funding by the Agenzia Spaziale Italiana ( ) under agreement no. 2018-23-HH.0 and support from INFN/Euclid Sezione di Roma.

DATA AVAILABILITY

The software lenstronomy was modified during this work, and the modifications are now part of the public version of that software. The full modification history of the code can be found in the pull request on GitHub: https://github.com/lenstronomy/ lenstronomy/pull/363. Our wrapper for lenstronomy built to produce and analyse the mock image catalogue presented here is publicly available at https://github.com/nataliehogg/analosis, along with the image data itself, while Jupyter notebooks for the rest of the analysis in this work can be found at https: //github.com/nataliehogg/los_proof_of_concept.

CREDIT AUTHORSHIP CONTRIBUTION STATEMENT

Natalie B. Hogg: Methodology, software, investigation, data curation, writing -original draft, visualisation. Pierre Fleury: Conceptualisation, methodology, software, validation, investigation, writing -original draft, visualisation. Julien Larena: Conceptualisation, methodology, writing -original draft. Matteo Martinelli: Resources, writing -review & editing.

APPENDIX A: MODIFYING lenstronomy

The lenstronomy software is a modular code designed primarily for modelling strong gravitational lenses, with numerous cosmological applications. The software contains a number of packages and subpackages which deal with the separate aspects of this task. For example, the LensModel package contains subpackages such as Profiles, which in turn contains a class for each lens mass profile which is available in lenstronomy.

The modifications to lenstronomy made for this work follow this modular structure. We introduced a new subpackage under the LensModel package called LineOfSight, which contains a new class, SinglePlaneLOS. This new class contains functions to compute the displacement angle, Hessian and Fermat potential in the LOS formalism. The class also inherits the ray shooting function (i.e. the computation of the source position 𝜷), the lensing potential function and the functions for computing the lens mass enclosed in a given sphere or projected radius and lens mass density from the native lenstronomy class SinglePlane, as these quantities are unchanged by LOS effects.

Furthermore, our new LineOfSight package contains a subpackage called LOSModels, inside which are the LOS and LOS_MINIMAL classes. The functions encoding the LOS effects are defined in the LOS class, and these functions are then called in the SinglePlaneLOS class. The additional parameters associated with the LOS effects are also defined in the LOS class, and the parameters associated with the minimal model are defined in LOS_MINIMAL. The additional parameters are: The final modifications we made to lenstronomy were to introduce additional unit tests which ensure the validity of the functions and the quantities which they return each time the code is built. We created two types of unit test. The first checks that when all the additional LOS parameters are fixed to zero, the SinglePlaneLOS functions return the same values as the SinglePlane functions for a given source and lens. The second exploits the correspondence that exists between the LOS formalism and multi-plane lensing with three shear planes. This allows us to construct unit tests which check that the LOS modifications return the same values as those returned by the equivalent multi-plane set-up using the MultiPlane class of lenstronomy.

For example, in Table A1, we show the output of the altered singleplane functions that we implemented in the SinglePlaneLOS class with the equivalent MultiPlane settings for a simple lens model consisting of an elliptical power law lens model with three shear planes. The lens and shear parameters, image positions, and the lens and source redshifts were chosen at random to produce these results.

The displacement angle computed by the SinglePlaneLOS class is identical to that computed in the multi-plane case to machine precision, Δ𝛼 ∼ 10 -16 arcsec. The elements of the Hessian matrix in the LOS case are the same as in the multi-plane case at the level of at least 10 -8 . Finally, the time delay resulting from the modified Fermat potential in the SinglePlaneLOS class matches that of the multi-plane case at the level of O (10 -11 ) days.

APPENDIX B: MOCK IMAGE CATALOGUE MODEL PARAMETERS

In this Appendix, we describe the full set of model parameters used to create the images in our mock catalogue. The numbers used to describe galaxy masses, half-light radii and Sérsic indices are freely inspired from Suess et al. (2019).

B1 Baryonic core: Sérsic ellipse

This component is modelled using the SERSIC_ELLIPSE_POTENTIAL in lenstronomy, which requires the angular half-light radius 𝑅 Sérsic , the Sérsic index 𝑛 Sérsic , the ellipticity components 𝑒 1 , 𝑒 2 , the effective convergence 𝑘 eff and the centre position 𝑥, 𝑦.

We begin by drawing a redshift for each lens from a uniform distribution, 𝑧 ∼ U (0.4, 0.6) and from this compute the angular diameter distance to each lens, 𝐷 od . Next, we draw the mass of the baryonic component randomly from a log-normal distribution, M ∼ ln N [ln(M mean ), ln(2)/2], where the mean mass M mean = 2 × 10 11 M . This distribution ensures that 95% of the lens galaxies we simulate have a mass which is at most a factor of two smaller or larger than the mean mass. From this we compute the half-light radius of the baryonic core in arcseconds,

We also draw the values for the Sérsic index from a log-normal distribution, 𝑛 Sérsic ∼ ln N [ln(𝑛 mean ), ln(1.5)/2], where 𝑛 mean = 4. We compute the ellipticity components using the following expressions,

)

where 𝑞 is the aspect ratio and 𝜙 the orientation angle of the ellipse. We draw the aspect ratio and orientations for each lens from uniform distributions, 𝑞 ∼ U (0.7, 1.0) and 𝜙 ∼ U (0, 2𝜋) respectively. We compute the effective convergence at the half-light radius,

where M is the lens mass in units of 𝑀 , Σ crit is the critical surface