On simple analytic models of microlensing amplification statistics

TL;DR

This paper reviews four simple analytic models for microlensing amplification statistics, introduces a corrected derivation for the multiplicative model, and evaluates their accuracy against simulations, highlighting limitations at high optical depths.

Contribution

It presents a new derivation of the multiplicative amplification distribution and compares four models' effectiveness against ray-shooting simulations.

Findings

All models perform well at low optical depth.

Models fail to capture features at high optical depth.

Lens-lens coupling is crucial at large optical depths.

Abstract

Gravitational microlensing is a key probe of the nature of dark matter and its distribution on the smallest scales. For many practical purposes, confronting theory to observation requires to model the probability that a light source is highly amplified by many-lens systems. This article reviews four simple analytic models of the amplification probability distribution, based on different approximations: (i) the strongest-lens model; (ii) the multiplicative model, where the total amplification is assumed to be the product of all the lenses' individual amplifications; (iii) a hybrid version of the previous two; and (iv) an empirical fitting function. In particular, a new derivation of the multiplicative amplification distribution is proposed, thereby correcting errors in the literature. Finally, the accuracy of these models is tested against ray-shooting simulations. They all produce excellent results as long as lenses are light and rare (low optical depth); however, for larger optical depths, none of them succeeds in capturing the relevant features of the amplification distribution. This conclusion emphasizes the crucial role of lens-lens coupling at large optical depths.

Figures

Tables

Table 1. Table 1: Parameters used for the simulations: optical depth  κ 𝜅 \kappa , number of lenses  N l subscript 𝑁 l N_{\text{l}} , number of rays shot  N 𝑁 N . We also indicate the mean amplification  ⟨ A ⟩ delimited-⟨⟩ 𝐴 \left\langle A\right\rangle , and the parameters  κ eff , b subscript 𝜅 eff 𝑏 \kappa_{\text{eff}},b of Rauch’s fitting formula ( 101 ).

$\kappa$	$N_{\text{l}}$	$N$	$⟨A⟩$	${\kappa }_{\text{eff}}$	$b$
0.01	5	${10}^8$	$1.019$	$0.0101$	$3340$
0.1	20	${10}^8$	$1.189$	$0.106$	$41.5$
0.5	50	${10}^8$	$1.859$	$0.719$	$2.53$
1	100	$2\times {10}^7$	$2.602$	$1.96$	$0.818$

Equations

$$\mathbf{\beta }=\mathbf{\theta }-\frac{{\epsilon }^2}{\mathbf{\theta }},$$

$${\epsilon }^2\equiv \frac{4Gm{D}_{\text{ds}}}{D_{\text{d}}{D}_{\text{s}}}$$

$${\theta }_{\pm }=\frac{1}{2}\left(\beta \pm \sqrt{{\beta }^2+4{\epsilon }^2}\right),$$

$$A_{\pm }^{-1}=\left|\det \frac{\partial \mathbf{\beta }}{\partial \mathbf{\theta }}\right|({\mathbf{\theta }}_{\pm })=\left|1-\frac{{\epsilon }^4}{{\theta }_{\pm }^4}\right|={\left|\frac{1}{2}\pm \frac{1}{2}\frac{u^2+2}{u\sqrt{u^2+4}}\right|}^{-1},$$

$$A(u)=A_{+}+A_{-}=\frac{u^2+2}{u\sqrt{u^2+4}}.$$

$$u^2=\frac{2A}{\sqrt{A^2-1}}-2.$$

$$u=\frac{b}{r_{\text{E}}},$$

$$r_{\text{E}}\equiv D_{\text{d}}\epsilon =\sqrt{4GM\mathcal{D}},\mathcal{D}\equiv \frac{D_{\text{d}}{D}_{\text{ds}}}{D_{\text{s}}}.$$

$$p_1(m,z,b)\mathrm{d}m\mathrm{d}z\mathrm{d}b=\frac{2b\mathrm{d}b}{R^2}\Theta (R-b)p(m,z)\mathrm{d}m\mathrm{d}z.$$

$$\frac{\mathrm{d}\Sigma }{\mathrm{d}z}\equiv \frac{{\mathrm{d}}^{3}N}{{\mathrm{d}}^{2}S\mathrm{d}z}=\frac{N}{\pi R^2}{\int }_0^{\infty }p(m,z)\mathrm{d}m,$$

$$\Sigma ={\int }_0^{z_{\text{s}}}\frac{\mathrm{d}\Sigma }{\mathrm{d}z}\mathrm{d}z=\frac{N}{\pi R^2}.$$

$$p_1(m,z,u)$$

$$p_1(u)=\frac{2u}{R^2}{\int }_0^{z_{\text{s}}}\mathrm{d}z{\int }_0^{M(u,z,R)}\mathrm{d}mp(m,z)r_{\text{E}}^2(m,z).$$

$$u_{\text{c}}=\frac{R}{\sqrt{4G{m}_{\text{max}}{\mathcal{D}}_{\text{max}}}},$$

$$\forall u\le u_{\text{c}}{p}_1(u)=\frac{2u⟨r_{\text{E}}^2⟩}{R^2},$$

$$⟨r_{\text{E}}^2⟩\equiv {\int }_0^{z_{\text{s}}}\mathrm{d}z{\int }_0^{m_{\text{max}}}\mathrm{d}mp(m,z)r_{\text{E}}^2(m,z).$$

$$p_{\text{s}}(u)\mathrm{d}u=\sum _{i=1}^N{p}_i(u)\mathrm{d}u\times \text{Prob}(u_{j\ne i}\ge u),$$

$$p_{\text{s}}(u)=N{p}_1(u){\left[1-{\int }_0^u\mathrm{d}{u}^{\prime }{p}_1(u^{\prime })\right]}^{N-1}.$$

$$N{p}_1(u)\to 2\kappa u,$$

$$\kappa \equiv \Sigma \pi ⟨r_{\text{E}}^2⟩={\int }_0^{z_{\text{s}}}\mathrm{d}z{\int }_0^{\infty }\mathrm{d}m\frac{{\mathrm{d}}^2\Sigma }{\mathrm{d}z\mathrm{d}m}\pi r_{\text{E}}^2(m,z).$$

$$p_{\text{s}}(u)=2\kappa u{\mathrm{e}}^{-\kappa u^2}=-\frac{\mathrm{d}}{\mathrm{d}u}\left({\mathrm{e}}^{-\kappa u^2}\right).$$

$$p_{\text{s}}(A)$$

$$=\frac{2\kappa }{(A^2-1{)}^{3/2}}\exp \left[-2\kappa \left(\frac{A}{\sqrt{A^2-1}}-1\right)\right],$$

$${⟨A⟩}_s\equiv {\int }_1^{\infty }\mathrm{d}AA{p}_{\text{s}}(A).$$

$${⟨A⟩}_{\text{s}}$$

$$=2\kappa {\mathrm{e}}^{2\kappa }{K}_1(2\kappa ),$$

$${⟨A⟩}_{\text{s}}=1+2\kappa +\mathcal{O}({\kappa }^2).$$

$$\kappa$$

$$=4\pi G{\int }_0^{z_{\text{s}}}\mathrm{d}z\frac{D_{\text{d}}(z)D_{\text{ds}}(z)}{D_{\text{s}}}\frac{{\rho }_{\text{l}}(z)}{(1+z)H(z)},$$

$$\kappa =\frac{3}{2}{H}_0^2{\int }_0^{{\chi }_{\text{s}}}\mathrm{d}\chi \frac{f_K(\chi )f_K({\chi }_{\text{s}}-\chi )}{f_K({\chi }_{\text{s}})}\frac{{\rho }_{\text{l}}(z)}{(1+z{)}^2}$$

$${⟨A⟩}_{\text{smooth}}=1+2\kappa +\mathcal{O}({\kappa }^2),$$

$${⟨A⟩}_{\text{s}}({\kappa }_{\text{s}})=2{\kappa }_{\text{s}}{\mathrm{e}}^{2{\kappa }_{\text{s}}}{K}_1(2{\kappa }_{\text{s}})={⟨A⟩}_{\text{true}}.$$

Full text

I Introduction

When, in 1936, Einstein was convinced by Mandl to publish the outcome of a short calculation about the lens-like action of a star by the deviation of light in the gravitational field Einstein (1936), he could not imagine the astronomical potential of his finding. That research note already proposed that the apparent brightness of a light source can be highly amplified by compact forms of matter on the line of sight. Somehow excavated 27 years later by Liebes Liebes (1964), this topic mostly concerned, until the mid-1980s, the effect of stars in galaxies acting as lenses for distant quasars Chang and Refsdal (1979); Gott (1981); Young (1981); Chang and Refsdal (1984). The term microlensing seems to have been introduced in 1986 by Paczyński Paczyński (1986a), who also proposed to use it as a probe of MAssive Compact Halo Objects (MACHO) in our galaxy Paczyński (1986b). This idea led to several surveys in the 1990s-2000s: the MACHO experiment Alcock et al. (1996); the Expérience pour la Recherche d’Objets Sombres (EROS) Afonso et al. (2003); and the still ongoing Optical Gravitational Lensing Experiment (OGLE) Wyrzykowski et al. (2015).

While focus was later displaced towards the detection of exoplanets Mao and Paczyński (1991); Bond et al. (2004); Ranc (2016), the concept of microlensing as a probe of the nature and distribution of dark matter never ceased to be enriched with new ideas, from multiply-imaged quasars Mediavilla et al. (2017) to supernova lensing Seljak and Holz (1999); Metcalf and Silk (1999, 2007). More recently, interest in that matter was revived by the observation of Icarus Kelly et al. (2018), a single star visible through cosmological distances thanks to a huge gravitational amplification, on the order of $10^{3}$. The very possibility of such an event was attributed to the disruption of a strong lens’ caustic by its own substructure Oguri et al. (2018); Diego et al. (2018), thereby opening a new branch in gravitational-lensing science Venumadhav et al. (2017); Diego (2019); Dai et al. (2018); Kaurov et al. (2019); Dai and Miralda-Escudé (2019).

In many concrete microlensing problems, a central observable is the amplification probability distribution function (PDF), $p(A)$. From the theoretical side, the difficulty consists in accurately relating this PDF to the matter distribution producing the amplifications. A significant research endeavor Katz et al. (1986); Deguchi and Watson (1987) was conducted in that direction during the 1980s-1990s, presumably with the hope to explain the variability of quasars with microlensing Schneider and Weiss (1987); Seitz and Schneider (1994); Seitz et al. (1994); Hawkins (2010); Goodman and Sun (2014). In that context, it was first understood in Ref. Peacock (1982) that $p(A)$ generally displays a long algebraic tail, of the form $A^{-3}$, when microlensing is at work. Several attempts to further characterize the full PDF Vietri and Ostriker (1983); Peacock (1986); Marchandon and Nottale (1991); Pei (1993) then led to a better understanding of the nonlinear coupling between lenses, when these are very numerous Nityananda and Ostriker (1984); Blandford and Narayan (1986); Schneider (1987a); Lee and Spergel (1990); Mao (1992); Kofman et al. (1997); Lee et al. (1997). Theoretical works were also guided by and tested with numerical simulations Kayser et al. (1986); Wambsganss (1992); Rauch (1991); Lewis and Irwin (1995). Finally, Refs. Petters et al. (2009a, b); Petters and Werner (2010) must be mentioned for a number of results on the mathematics of gravitational lensing.

On the shoulders of giants, the goal of the present article is rather modest. The aforementioned works on amplification statistics have led to several modeling techniques, the simplest ones being currently used to set or forecast constraints on the nature of dark matter Zumalacárregui and Seljak (2018); García-Bellido et al. (2018); Diego (2019). However, it seems that the performance of these simple models for $p(A)$ has never been properly assessed. We propose here to fill this gap. Such a comparative analysis will also be the occasion to clarify and correct some theoretical points of the existing literature.

Specifically, four models will be reviewed throughout this article. In Sec. II, we introduce microlensing fundamentals, and consider the strongest-lens model, where the amplification due to a set of lenses is strictly due to the strongest one. In Sec. III we consider the multiplicative approach, where the total amplification, due to many lenses, is assumed to be the product of all their individual amplifications. We use this opportunity to correct the corresponding derivation of $p(A)$ with respect to earlier works. We also propose a hybrid model between the strongest-lens and the multiplicative models. Finite-source corrections to the previous calculations are considered in Sec. IV. Finally, in Sec. V, the above two approaches are compared with an empirical expression for $p(A)$ and confronted to numerical ray-shooting simulations. We conclude in Sec. VI.

We adopt units in which the speed of light is unity, $c=1$. Bold symbols, such as $\boldsymbol{\theta}$, stand for two-dimensional vectors. A hatted vector $\hat{\boldsymbol{\theta}}$ denotes the unit counterpart to the nonhatted one, $\hat{\boldsymbol{\theta}}=\boldsymbol{\theta}/\theta$, where $\theta\equiv||\boldsymbol{\theta}||$ is the Euclidean norm of $\boldsymbol{\theta}$.

II Microlensing by individual lenses

II.1 Amplification by a point lens

Consider an infinitesimal light source lensed by a point mass. Throughout this article, we will work in the weak gravitational field and geometric optics regimes, assuming the small-angle and flat-sky approximations—see Ref. Fleury et al. (2019) for a comparative discussion of these assumptions. Let the lens be at the origin of a celestial coordinate system, and call $\boldsymbol{\beta}$ the unlensed position of the source (see Fig. 1). The image position $\boldsymbol{\theta}$ then satisfies the lens equation Schneider et al. (1992)

[TABLE]

where $1/\boldsymbol{\theta}\equiv\hat{\boldsymbol{\theta}}/\theta$, while

[TABLE]

denotes the Einstein radius of the lens. In Eq. (2), $m$ denotes the mass of the lens; $D_{\text{d}},D_{\text{s}}$ are respectively the angular-diameter distance of the lens (deflector), source, as seen from the observer; and $D_{\text{ds}}$ is the distance to the source as seen from the lens. The precise expression of these distances depend on a choice of background, i.e., a fiducial no-lensing situation.

The lens equation (1) has two solutions $\boldsymbol{\theta}_{\pm}=\theta_{\pm}\hat{\boldsymbol{\beta}}$

[TABLE]

which are the positions of the two images of the source at $\boldsymbol{\beta}$. The luminous intensity of each image reads $I_{\pm}=A_{\pm}I_{0}$, where $I_{0}$ is the unlensed intensity, i.e. the source’s apparent luminosity in the absence of lensing. Due to surface-brightness conservation, the individual inverse amplifications (or magnifications) read

[TABLE]

with $u\equiv\beta/\varepsilon$.

If the two images are not resolved by the telescope, i.e. if $|\theta_{+}-\theta_{-}|$ is smaller than the telescope’s resolution, they are called microimages. Lensing then only manifests through the apparent amplification of the luminosity of the macro-image, $I=I_{+}+I_{-}$. This is a microlensing event, and the total amplification reads

[TABLE]

This relation can also be inverted to get

[TABLE]

Finally, it is useful to rewrite $u$ as

[TABLE]

where $b=D_{\text{d}}\beta$ is the physical impact parameter of the unlensed light path, and

[TABLE]

Thus, $u$ must be understood as a reduced impact parameter, i.e. expressed in units of the lens’ cross-sectional radius. Indeed, it is customary to designate $\pi r_{\text{E}}^{2}$, or $\pi\varepsilon^{2}$, the cross section of the lens. This is because $\pi\varepsilon^{2}$ is the area of the sky where a light source gets amplified by a factor $A>A(1)\approx 1.34$.

II.2 From one to many: the strongest-lens approximation

The single-lens case is the only one which is analytically solvable. However, in many physically relevant situations, a given source may be affected by many lenses. These cases include microlensing by planetary systems (lenses are a star and its planets) Mao and Paczyński (1991); but also nearby supernovae or quasars observed through galaxies (lenses are stars and globular clusters) Vietri and Ostriker (1983); or distant sources observed through clusters of galaxies (lenses are galaxies). This question is especially relevant in a scenario where a significant fraction of dark matter would be made of compact objects, such as primordial black holes Carr and Hawking (1974); Garcia-Bellido et al. (1996); Bird et al. (2016); Carr et al. (2016); Clesse and García-Bellido (2018).

How to model the combined effect of many lenses? Specifically, we aim, here, to evaluate the probability density function (PDF) of the microlensing amplification, $p(A)$, due to all the lenses potentially located between the source and the observer. If the lenses are rare and not too massive, then the sum of their individual cross sections $\pi\varepsilon^{2}$ occupy a small fraction of the sky. Equivalently, $\varepsilon\ll\Delta\theta$, where $\Delta\theta$ is the typical angular separation between two lenses. In that context, referred to as the low-optical-depth regime, a typical source is mostly affected by a single lens: the one with smallest reduced impact parameter $u$.

To the best of our knowledge, the explicit expression of $p(A)$ in the strongest-lens approximation was first derived by Peacock in 1986 Peacock (1986). Another derivation, originally due to Nottale, can can also be found in Ref. Marchandon and Nottale (1991). In this section, we propose an alternative, step-by-step, calculation of $p(A)$, whose formalism will be useful in the remainder of the article.

II.2.1 Setup and notation

Consider, as depicted in Fig. 2, a tube with physical radius $R$ between the observer and the source, filled with a finite number $N$ of lenses. These lenses may have different masses, and their distribution may be inhomogeneous along the line of sight. Our only assumption is that, on each disk $z=\mathrm{cst}$, their spatial distribution is Poissonian.

For a single lens, the probability that its mass is in $[m,m+\mathrm{d}m]$, its redshift in $[z,z+\mathrm{d}z]$, and its physical impact parameter with respect to the line of sight is in $[b,b+\mathrm{d}b]$, reads

[TABLE]

In Eq. (9), the Heaviside function $\Theta$ ensures that the lens lies within the tube, while $p(m,z)$ is an arbitrary distribution of masses and redshifts, which may be correlated. The macroscopic number density of lenses, per unit area and redshift, is related to $p(m,z)$ as

[TABLE]

while the mean surface density, obtained by stacking all the planes $z=\mathrm{cst}$, simply reads

[TABLE]

When, at the very end of the calculation, we will take $N,R\rightarrow\infty$, the surface density $\Sigma$ will be kept constant, so that $N=\mathcal{O}(R^{2})$.

II.2.2 Amplification probability for one lens

Let us first consider the case of a single lens in the entire tube, and determine the corresponding amplification PDF, $p_{1}(A)$. Because the $A$ only depends on the reduced impact parameter $u$, it is equivalent to determine $p_{1}(u)$. For that purpose, the first step consists in translating Eq. (9) in terms of $u=b/r_{\text{E}}(m,z)$,

[TABLE]

The PDF of $u$ only is then obtained by marginalizing over $m,z$. In that operation, the subtlety is how to handle the Heaviside function, which depends on both $z,m$ via $r_{\text{E}}$. The simplest way consists in encoding it in an upper limit for the integration over $m$. We define $M(u,z,R)$ such that $R=u\,r_{\text{E}}(M,z)$, which represents the mass above which the impact parameter $b$ should be larger than $R$ to have the correct value of $u$. With this notation,

[TABLE]

Of course, in reality, the lens mass distribution does not extend to infinity; there exists some maximum lens mass $m_{\text{max}}$ such that $p(m>m_{\text{max}},z)=0$. Since the function $u\mapsto M(u,z,R)\propto u^{-2}$ is monotonically decreasing from $(0,\infty)$ to $(0,\infty)$, there exists a critical value $u_{\text{c}}(R)$ such that $\forall z\in[0,z_{\text{s}}]\quad M(u\leq u_{\text{c}},z,R)\geq m_{\text{max}}$. Its explicit expression is easily found to be

[TABLE]

where $\mathcal{D}_{\text{max}}$ is the maximum value of $\mathcal{D}$, which is typically reached when the lens lies midway between the observer and the source. Thus, for $u\leq u_{\text{c}}$, we can substitute $M$ with $m_{\text{max}}$ in Eq. (13) and get

[TABLE]

where we introduced

[TABLE]

It is not necessary to explicitly determine the case $u\geq u_{\text{c}}$ in order to proceed.

II.2.3 The strongest out of $N$ lenses

Let us now consider the case where $N$ lenses are in the tube, and let us determine the PDF of the strongest amplification. Again, since $A$ only depends on $u$, the strongest lens is the one with the smallest reduced impact parameter $u$. Let us call $p_{\text{s}}(u)$ the associated PDF (subscript “s” stands for “strongest”). Assuming that the $N$ lenses are independent, we have

[TABLE]

where $p_{i}$ is the unconstrained PDF of $u$ for the $i$th lens. Since all the lenses are characterized by the same PDF $p_{1}(u)$, Eq. (17) becomes

[TABLE]

II.2.4 Infinite number of lenses

We now take the limit of an infinite Universe containing an infinite number of lenses, i.e. $N,R\rightarrow\infty$, while keeping the surface density $\Sigma=N/(\pi R^{2})$ fixed. In that limit, $u_{\text{c}}\rightarrow\infty$, and hence, for any $u$,

[TABLE]

where we introduced the microlensing optical depth

[TABLE]

Besides, the limit of the square-bracket term in (18) yields an exponential, so that

[TABLE]

This result is easily translated in terms of amplifications using Eq. (6), and yields the PDF of the strongest amplification

[TABLE]

in agreement with Refs. Peacock (1986); Marchandon and Nottale (1991).

II.3 Discussion: optical depth, Weinberg, and Peacock

The most striking property of Eq. (22) is that it only depends on a single parameter: the optical depth $\kappa$. This parameter notably controls the amplitude of the high-amplification tail of the distribution, $p_{\text{s}}(A\gg 1)\sim 2\kappa/A^{3}$, which is a well-known behavior Vietri and Ostriker (1983). Alternatively, one could say that the entire PDF is controlled by its mean,

[TABLE]

With the change of variable $\cosh x=A/\sqrt{A^{2}-1}$, one finds

[TABLE]

where $K_{1}$ denotes the first modified Bessel function of the second kind.

Since $\kappa\sim(\varepsilon/\Delta\theta)^{2}$, we expect the strongest-lens approximation to hold only for $\kappa\ll 1$. In that regime, the mean amplification reads

[TABLE]

This result deserves to be noticed, because it relates to the long-standing question of flux conservation in gravitational lensing. As shown by Weinberg in 1976 Weinberg (1976), the mean amplification caused by a very sparse set of isolated, Poisson-distributed, point lenses is identical to the amplification which would be observed if the matter constituting these lenses was homogeneously distributed in space. Equation (27) generalizes Weinberg’s result to a nonhomogeneous distribution of lenses.

To understand that claim, it is useful to rewrite $\kappa$ by substituting the expression (8) of $r_{\text{E}}$ into Eq. (20); we obtain

[TABLE]

where $\rho_{\text{l}}(z)$ is the contribution of the point lenses to the mean mass density111We are not necessarily talking about a homogeneous Universe here; in particular, $\rho_{\text{l}}(z)\not\propto(1+z)^{-3}$ in general. at $z$, while a factor $(1+z)H(z)=\mathrm{d}z/\mathrm{d}d$ appears to convert redshifts into proper distances. Furthermore, assuming that angular-diameter distances $D_{\text{d}},D_{\text{ds}},D_{\text{s}}$ can be computed as in a homogeneous and isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) Universe, we find

[TABLE]

where $\chi$ denotes comoving distances, $K$ is the Universe’s spatial curvature parameter, and $f_{K}(\chi)\equiv\sin(\sqrt{K}\chi)/\sqrt{K}$. Equation (30) is the expression of the weak-lensing convergence which would be due to our point lenses, if their mass were smoothly distributed, instead of being concentrated into compact objects. In that case, we would have

[TABLE]

which, indeed, agrees with Eq. (27) in the limit $\kappa\ll 1$.

In Ref. Peacock (1986), Peacock proposed to heuristically extend the applicability of $p_{\text{s}}(A)$ to higher optical depths, $\kappa\sim 1$, by conjecturing flux conservation. The idea consists in treating the optical depth involved in $p_{\text{s}}(A)$ as a free parameter, denoted $\kappa_{\text{s}}$, and fix it by ensuring that $\left\langle A\right\rangle_{\text{s}}(\kappa_{\text{s}})$ produces the true mean amplification, i.e.222Peacock seems to have assumed $\left\langle A\right\rangle_{\text{true}}=(1-\kappa)^{-2}$, because he found, e.g., $\kappa_{\text{s}}=199$ for $\kappa=0.8$. This choice delivered results in good agreement with the numerical simulations of Paczyński Paczyński (1986a). However, the results of the present article do not support this conclusion; in Sec. V we rather find $\left\langle A\right\rangle_{\text{sim}}\approx(1+\sqrt{1+4\kappa})^{2}/4$, although the discrepancy shall be explained by the finite extension of our simulated map.

[TABLE]

As will be seen in Sec. V, this procedure would not succeed to capture the features of $p(A)$ due to collective effects for $\kappa\sim 1$.

III Multiplicative amplifications

The strongest-lens approximation is arguably crude, because it completely neglects the long-distance effects due to the other lenses. In order to allow for all the lenses together, one may conjecture that amplifications are multiplicative Vietri and Ostriker (1983); namely, the total amplification caused by two lenses reads $A=A_{1}A_{2}$, where $A_{1},A_{2}$ are the individual amplifications that each lens would produce in the absence of the other.

III.1 Can we really multiply amplifications?

An intuitive justification of the multiplicative model is the following. Suppose that the two lenses are well separated along the line of sight, lens 1 being the closest to the observer and lens 2 the closest to the source. Let $\boldsymbol{\theta}$ be an image. Each lens is individually endowed with an image-to-source mapping of the form (1), $\boldsymbol{\beta}_{1}(\boldsymbol{\theta}),\boldsymbol{\beta}_{2}(\boldsymbol{\theta})$. Considering their successive effect, and assuming that, for lens 2, $\boldsymbol{\beta}_{1}(\boldsymbol{\theta})$ plays the role of an intermediate image, the combined mapping would be $\boldsymbol{\beta}(\boldsymbol{\theta})=\boldsymbol{\beta}_{2}(\boldsymbol{\beta}_{1}(\boldsymbol{\theta}))$. The chain rule then yields the amplification as

[TABLE]

Leaving aside the important fact that $A_{2}$ must be evaluated at a non-trivial intermediate image position, the above reasoning is incorrect anyway. The reason is that gravitational lensing is a somewhat non-local phenomenon: the properties of a given lens depend on everything that happens to light before and after it. For example, the exact lens map corresponding to the combination of two point lenses is

[TABLE]

where $\varepsilon_{1},\varepsilon_{2}$ are the Einstein radii of the two lenses, $\boldsymbol{\lambda}_{1},\boldsymbol{\lambda}_{2}$ their positions, and $\varepsilon_{12}\equiv\sqrt{4Gm_{1}D_{\text{d}_{1}\text{d}_{2}}/D_{\text{o}\text{d}_{1}}D_{\text{o}\text{d}_{2}}}$. Equation (34) must be compared with

[TABLE]

Since $\varepsilon_{1}\not=\varepsilon_{12}$, the ansatz $\boldsymbol{\beta}_{2}(\boldsymbol{\beta}_{1}(\boldsymbol{\theta}))$ does not properly account for lens-lens coupling. This fundamentally prevents amplifications from being multiplicative. Another, equally important, obstacle, is that one must sum the individual amplifications of each image of a given source, which requires to determine these images in the first place.

Of course, when the multiplicative model was first proposed by Ref. Vietri and Ostriker (1983) in 1983, its authors were aware that it was only an approximation. A year later, a major aspect of lens-lens coupling was emphasized by Refs. Chang and Refsdal (1984); Nityananda and Ostriker (1984), namely shear. If a light source is (even weakly) sheared by a first lens, then the effect of a second lens turns out to be quite different from the no-shear case. In particular, the geometry of its caustics is significantly affected, as well as the amplification statistics. This problem has been thoroughly investigated in the 1980s-1990s, both analytically and numerically Kayser et al. (1986); Schneider (1987a); Mao (1992); Kofman et al. (1997); Lee et al. (1997). More recently, this coupling between the weak-lensing effect of the cosmic web with the properties of strong lenses Schneider (1997); Birrer et al. (2017) has been proposed as a possible measure of cosmic shear using Einstein rings Birrer et al. (2018).

III.2 Practical interest of multiplicativity

Albeit inexact, the multiplicative model is practically very convenient, as far as statistics are concerned. Let two lenses, or groups of lenses, generate random amplifications $A_{1},A_{2}$. Assuming multiplicativity, the total amplification reads $\ln A=\ln A_{1}+\ln A_{2}$. Hence, if $L\equiv\ln A$, then the PDF of $L$ reads

[TABLE]

where $P_{12}(L_{1},L_{2})$ is the joint PDF of the individual logarithmic amplifications. If, furthermore, $A_{1}$ and $A_{2}$ are independent, then $P_{12}(L_{1},L_{2})=P_{1}(L_{1})P_{2}(L_{2})$, so that Eq. (36) becomes a convolution product, $P=P_{1}*P_{2}$. This equality can be translated in terms of amplifications, using that

[TABLE]

which yields

[TABLE]

Equation (38) is a simple prescription to combine independent multiplicative amplifications, which will be useful for the remainder of this section.

III.3 Derivation of the amplification PDF

Assuming multiplicativity, an exact expression for the amplification PDF, $p_{\text{m}}(A)$ (subscript “m” stands for “multiplicative”) was derived in 1993 by Pei Pei (1993), hereafter P93. However, this derivation contains a few errors and inaccuracies—we refer the curious reader to Appendix A where Pei’s method is reproduced and commented. This encouraged us to propose an alternative proof in the present section. Fortunately, the final result is left unchanged.

III.3.1 From 1 to $N$ lenses in a finite Universe

Consider the same setup as described in Sec. II.2.1, with $N$ lenses in a finite tubular Universe. As suggested by the discussion of Sec. III.2, when amplifications are multiplicative, it is more convenient to work with their logarithm. Let $P_{N}(L)$ be the PDF of $L\equiv\ln A=\ln A_{1}+\ldots\ln A_{N}$, we then have

[TABLE]

where $P_{1}(\ln A)=Ap_{1}(A)$ logarithmic amplification PDF for a single lens. This quantity can be directly deduced from the expression (15) of $p_{1}(u)$, namely,

[TABLE]

as long as $A\geq A_{\text{c}}\equiv A(u_{\text{c}})$, where $u_{\text{c}}$ is given by Eq. (14). Again, it is not necessary to explicitly determine the case $A\leq A_{\text{c}}$ in order to proceed.

The convolution product (39) is more easily handled in Fourier space, where it becomes a regular product. We adopt the following convention for Fourier transforms:

[TABLE]

so that $\tilde{P}_{N}(K)=\tilde{P}_{1}^{N}(K)$. Therefore, by taking the inverse Fourier transform, and restoring the normal amplification variable $A=\exp L$, we obtain the formal expression

[TABLE]

III.3.2 An infinity of lenses

Let us finally determine the limit $p_{\text{m}}(A)$ of $p_{N}(A)$ as both the number $N$ of lenses and the size $R$ of the tube go to infinity. For that purpose, the first step consists in manipulating the logarithmic Fourier transform of $p_{1}$ as

[TABLE]

For $R\rightarrow\infty$, we have $u_{\text{c}}\rightarrow\infty$, and hence $A_{\text{c}}\rightarrow 1$. Specifically,

[TABLE]

Hence, the first integral of Eq. (45) can be expanded as

[TABLE]

where, in the second line, we used the normalization of $p_{1}(A)$. Substituting the expression of $p_{1}(A\geq A_{\text{c}})$, we find

[TABLE]

The lower limit $A_{\text{c}}$ of the integral of Eq. (49) can actually be replaced with $1$, because

[TABLE]

This leaves us with the analytically solvable integral

[TABLE]

which does not depend on $N,R$.

The last step consists in raising $\tilde{P}_{1}(K)$ to the $N$th power, which yields

[TABLE]

We recognize the optical depth $\kappa=N\langle r_{\text{E}}^{2}\rangle/R^{2}$, while the remainder $\mathcal{O}(N/R^{4})$ goes to zero as $N,R\rightarrow\infty$. Substituting that result into Eq. (43), we conclude that

[TABLE]

where $\Gamma$ denotes the usual Gamma function. The above Eq. (53) is equivalent to Eq. (29) in P93, although the present definition (20) of $\kappa$ is slightly more general, because it allows the lenses to have different masses, and to be inhomogeneously distributed along the line of sight.

III.4 Discussion: low-optical-depth behavior and moments

Albeit complicated due to its highly oscillatory integrand, the expression (53) of $p_{\text{m}}(A)$ has a simple limit in the low-optical depth limit $\kappa\ll 1$. Indeed, expanding the exponential, we have

[TABLE]

in agreement with P93. The latter result is obtained by using that, in Eq. (55), integration over $K$ yields an integrand proportional to $\delta(\ln A-\ln A^{\prime})-\delta(\ln A)$. Note that, as expected, this low-$\kappa$ behavior exactly coincides with the one of $p_{\text{s}}(A)$. It is not trivial, however, to which extent the $\mathcal{O}(\kappa^{2})$ terms can be neglected, especially for $A-1\ll 1$, because $(A^{2}-1)^{3/2}$ is nonintegrable on $(1,\infty)$.

The moments of $p_{\text{m}}(A)$ are more conveniently determined starting from Eq. (44), replacing $A^{-\mathrm{i}K}$ with $A^{n}$, and then following the same calculation as in Sec. III.3.2. The result is

[TABLE]

For $n\geq 2$, the integral diverges because of the upper limit $A\rightarrow\infty$. The mean amplification is, therefore, the only nonzero moment of $p_{\text{m}}(A)$, and reads

[TABLE]

again in agreement with P93.

For $\kappa\ll 1$, we find $\left\langle A\right\rangle_{\text{m}}=1+2\kappa+\mathcal{O}(\kappa^{2})$ just like in the strongest-lens model. It is instructive to further compare the behaviors of $\left\langle A\right\rangle_{\text{s}}$ and $\left\langle A\right\rangle_{\text{m}}$; in particular,

[TABLE]

anticipating on the results of Fig. 7. These inequalities are a hint that the strongest-lens model may systematically underestimate amplifications, while the multiplicative model may systematically overestimate them.

III.5 Improving the multiplicative model: a hybrid approach

One reason why the multiplicative model overestimates amplifications may be that it virtually accounts for two images per lens, which is not realistic. In order to simply allow for the combined effect of many lenses, without dramatically overestimating it, we propose the following hybrid approach between the strongest-lens and multiplicative model.

III.5.1 Hybrid model

Suppose that a light source is strongly affected by one lens, and weakly affected by all the others. Let us assume that each weak lens produces a single image—the principal image, with amplification $A_{+}$. Let us also neglect lens-lens coupling, so that the amplification $A_{\text{s}}$ of the strong lens is left unchanged. Then, the total amplification reads $A_{\text{s}}A_{\text{w}}^{+}$, where $A_{\text{w}}^{+}$ is the product of all the principal amplifications of the weak lenses.

Furthermore, let us adopt a mean-field approach333It is actually possible to proceed without this simplifying assumption. However, the final analytic result is too complicated to have any practical interest. We thus restrict to the mean-field case in this section., and replace the stochastic contribution $A_{\text{w}}^{+}$ by its statistical average $\bar{A}_{\text{w}}^{+}$. Our hybrid model is then defined by

[TABLE]

where $\bar{A}_{\text{w}}^{+}$ depends on $A_{\text{s}}$ because of the constraint that the lenses contributing to $A_{\text{w}}$ are weaker than the main strong lens producing $A_{\text{s}}$. In this model, the stochasticity of $A$ is entirely controlled by the strongest amplification $A_{\text{s}}$. In other words, the amplification PDF of the hybrid model reads

[TABLE]

The remaining task thus consists in computing $\bar{A}_{\text{w}}^{+}(A_{\text{s}})$.

III.5.2 Calculation of the mean weak amplification

Let us consider again the finite setup described in Sec. II.2.1. Suppose that the strongest lens has a reduced impact parameter $u_{\text{s}}$. Then all the other lenses must satisfy $u\geq u_{\text{s}}$, so that for any one of them

[TABLE]

which can be marginalized over $m,z$ to get $p_{1}(u|u_{\text{s}})$. The result can then be converted in terms of amplifications using Eq. (4), and we find

[TABLE]

where $A_{\text{c}}^{+}=A_{+}(u_{\text{c}})$, $u_{\text{c}}$ being given by Eq. (14) to ensure that $b\leq R$. The second Heaviside function corresponds to $u\leq u_{\text{s}}$; note that this is not equivalent to $A_{+}\leq A_{\text{s}}$ because $A_{+}(u)=1/2+A(u)/2\not=A(u)$. The mean weak amplification due to a single lens besides the strongest lens thus reads

[TABLE]

In the finite tubular Universe, the average weak amplification reads $\bar{A}_{\text{w}}^{+}=\left\langle A_{+}\right\rangle^{N}$. In the limit $N,R\rightarrow\infty$, following a similar computation as in Sec. III.3.2, we finally get

[TABLE]

Summarizing, in the hybrid model the amplification PDF reads

[TABLE]

In practice, one has to numerically invert the relation $A(A_{\text{s}})$ in order to compute $p_{\text{h}}(A)$.

III.5.3 Asymptotic behavior and moments

In the low-optical-depth regime, $A_{\text{s}}=A+\mathcal{O}(\kappa)$, and hence we find again

[TABLE]

just like $p_{\text{s}},p_{\text{m}}$. The amplitude of the high-amplification tail is also easily obtained whatever $\kappa$: for $A\gg 1$, $A\approx\mathrm{e}^{\kappa}A_{\text{s}}$, so that $p_{\text{h}}(A\gg 1)\approx 2\kappa\mathrm{e}^{2\kappa}/A^{3}=\mathrm{e}^{2\kappa}p_{\text{s}}(A\gg 1)$. Thus, in the optically-thick regime, the tail is significantly higher than for the strongest-lens approximation.

Just like in the other approximation schemes, the $A^{-3}$ algebraic tail for large amplifications implies that only the first moment of $p_{\text{h}}(A)$, i.e. the mean amplification, is finite. It can be computed analytically as

[TABLE]

where, from Eq. (72) to Eq. (73), we defined the variable $x=\sqrt{(A_{\text{s}}+1)/(A_{\text{s}}-1)}$, and $\Gamma(a,y)$ denotes an incomplete Gamma function,

[TABLE]

IV Extended sources

Hitherto, all the calculations were performed assuming that the light sources undergoing microlensing were point-like. However, it is quite instructive to investigate the case where this assumption is relaxed, especially if one is interested about the collective effect of many low-mass lenses. Seminal works on that specific topic include Refs. Nottale (1986); Schneider and Weiss (1987); Schneider (1987b); Marchandon and Nottale (1991).

In full generality, an extended source $\mathcal{S}$ can be modelled as a collection of infinitesimal patches. The apparent luminous intensity444We call luminous intensity the electromagnetic power per unit area received by the detector. of such an infinitesimal patch $\mathrm{d}^{2}\boldsymbol{\beta}$ of intrinsic intensity $\mathrm{d}^{2}I_{\text{s}}$ being $A(\boldsymbol{\beta})\mathrm{d}^{2}I_{\text{s}}$, the total apparent intensity of the extended source reads

[TABLE]

The quantity $\mathrm{d}^{2}I_{\text{s}}/\mathrm{d}^{2}\boldsymbol{\beta}$ is called the specific intensity of the source. From the above, we immediately conclude that the net amplification of the extended source reads

[TABLE]

In other words, accounting for the finite extension of sources is equivalent to smoothing amplification maps.

IV.1 Amplification of a disk source by a single one point lens

Let us assume, for simplicity, that $\mathcal{S}$ is a disk with apparent angular radius $\sigma$, and whose specific intensity is homogeneous. In that case, Eq. (77) simplifies into

[TABLE]

If the source is lensed by a single point lens, an analytical expression for $A_{\text{ext}}$ was found by Ref. Witt and Mao (1994) as below, and represented in Fig. 3

[TABLE]

where, as before, $u$ denotes the reduced impact parameter of the source; the new parameter $r\equiv\sigma/\varepsilon$ describes the size of the source in units of the lens’ Einstein radius. The functions $\mathrm{E},\mathrm{F},\Pi$ are the complete elliptic integrals of the first, second, third type, respectively, and

[TABLE]

Note that, in the above expressions, we used Wolfram’s convention555https://www.wolframalpha.com/input/?i=ellipticE for elliptic integrals, which differs from the Gradshteyn & Ryzhik convention used in Ref. Witt and Mao (1994).

The point-source case is shown to be recovered for $r\rightarrow 0$. Contrary to what happens in that case, when $r\not=0$, the amplification is bounded up. The maximal amplification happens when $u=0$, and reads

[TABLE]

As expected, $A_{\text{max}}$ diverges for $r\rightarrow 0$.

IV.2 Impact on amplification statistics

Let us now investigate how the finite extension of sources affect amplification statistics. We will restrict that analysis to the strongest-lens and multiplicative models.

IV.2.1 In the strongest-lens model

The computation of $p(A)$ with extended sources can be adapted from the point-source derivation of Sec. II.2. In doing so, the main difficulty is that $A$ is no longer a function of the reduced impact parameter $u$ only. We must now account for the dependence in $r=\sigma/\varepsilon$.

Let us start with the simplified case where all the lenses have the same Einstein radius $\varepsilon$. Then $r$ only takes a single value, and hence all the steps of Sec. II.2 can be repeated identically, except that $A(u)$ must be replaced by $A_{\text{ext}}(u,r)$. This yields

[TABLE]

where $A\mapsto u(A,r)$ is the inverse of $u\mapsto A_{\text{ext}}(u,r)$. Unfortunately, even in the special case of a homogeneous disk-source, Eq. (79) cannot be inverted analytically.

How is that result generalized to lenses with a distribution of Einstein radii? What prevents us from following the derivation of Sec. II.2 is that the strongest lens is no longer the lens with the smallest $u$; indeed, if the latter turns out to be characterized by a large ratio $r$, its amplification may be significantly damped, so that another lens with higher $u$ may be more efficient.

Nevertheless, this difficulty can be circumvented as follows. Let us organize the lenses in groups characterized by a given Einstein radius. Each group $i$ is thus equipped with an $\varepsilon_{i}$ (hence an $r_{i}=\sigma/\varepsilon_{i}$), and its importance is quantified by a partial optical depth $\kappa_{i}$. The strongest amplification is then the strongest of the strongest within each group, and hence the associated PDF reads

[TABLE]

where $\kappa$ is the sum of the partial optical depths of all the groups of lenses, and we introduced the optical-depth-weighed impact parameter

[TABLE]

in the continuous limit.

IV.2.2 In the multiplicative model

With the multiplicative approach, we have found in Sec. III.3.1 that $p_{\text{m}}(A)$ takes the form of an inverse Fourier transform

[TABLE]

with

[TABLE]

These expressions turn out to be independent of the actual expression of $A(u)$, and thus partly generalize from the point-source case to the extended-source case. Specifically, if all the lenses have the same Einstein radius $\varepsilon$, we simply have

[TABLE]

where, again, $u(A;r)$ denotes the inverse of $A_{\text{ext}}(u;r)$.

For a general population of lenses, we can proceed similarly to Sec. IV.2.1 and organize them in groups of identical Einstein radii. Due to the multiplicative assumption, the Fourier transform of the logarithmic total amplification reads

[TABLE]

where $u_{\text{eff}}^{2}$ is the same as defined in Eq. (86). Contrary to the point-source case, there is no analytic expression of the above integral. Transforming the integral over $A$ into an integral over $u^{2}$, we conclude that

[TABLE]

where $A_{\text{eff}}(u)$ is the inverse of $u_{\text{eff}}(A)$.

IV.3 Discussion

Just like in the infinitesimal-source case, both the strongest-lens model and the multiplicative model coincide in the limit $\kappa\ll 1$. Expanding the exponential of Eq. (93) at first order in $\kappa$, we find for $A>1$

[TABLE]

Beyond that regime, the strongest-lens and multiplicative models have qualitatively distinct behaviors, especially at large amplifications. Indeed, while $p_{\text{s}}^{\text{ext}}(A)$ experiences a sharp drop to zero as $A$ reaches its maximum value, there is no such drop for $p_{\text{m}}^{\text{ext}}(A)$. Instead, the high-$A$ tail of the latter decreases more rapidly than in the infinitesimal-source case, without being totally annihilated.

The above analysis of finite-size effects shows that they are controlled by the dimensionless parameter $r=\sigma/\varepsilon=R_{\text{s}}/r_{\text{E}}$, where $R_{\text{s}}$ is the physical size of the source. It is instructive to evaluate the order of magnitude of $r$ as expected in relevant setups. For a type-Ia supernova at $z=1$ lensed by an object of mass $M$, we find

[TABLE]

where $\mathcal{D}$ is the distance ratio defined in Eq. (8). In that situation, $r$ becomes comparable to unity for $M\lesssim 10^{-3}M_{\odot}$. For a quasar at $z=2$, we have

[TABLE]

so that finite-size effects already kick in for massive lenses.

V Comparison with numerical simulations

How well do the analytic models $p_{\text{s}}(A)$, $p_{\text{m}}(A)$, or $p_{\text{h}}(A)$ reproduce the actual amplification PDF for a given optical depth $\kappa$? This question can be partially addressed using numerical simulations. Section V.1 introduces the inverse ray-shooting technique used for that purpose; Secs. V.2, V.3 compare the three models’ performance in reproducing respectively the mean amplification and full PDF; and Sec. V.4 focuses on the impact of extended sources.

V.1 Numerical method: inverse ray shooting

V.1.1 Principle and setup

Inverse ray shooting Kayser et al. (1986); Mediavilla et al. (2006) is a conceptually simple method to generate amplification maps, PDFs, and lightcurves. Its principle relies on the fact that the lens map expresses the position of the source $\boldsymbol{\beta}$ of a given image $\boldsymbol{\theta}$ (whence the name inverse ray shooting). Consider $N$ fictitious images $\boldsymbol{\theta}$, arranged on a regular grid666Another option is to pick image positions randomly, but the resulting shot noise significantly reduces the accuracy of the code., with homogeneous surface density $n_{\text{im}}=\mathrm{d}^{2}N/\mathrm{d}^{2}\boldsymbol{\theta}$. Map these $N$ images to their respective sources, using $\boldsymbol{\beta}(\boldsymbol{\theta})$, so that the source plane is inhomogeneously filled with points. Call $n_{\text{s}}(\boldsymbol{\beta})=\mathrm{d}^{2}N/\mathrm{d}^{2}\boldsymbol{\beta}$ its surface density. Since the number of points is conserved during the lens mapping,

[TABLE]

Thus, the surface density of the fictitious sources directly tells us about the amplification in the source plane.

The simulations presented in this article were performed with identical lenses, randomly distributed on a plane (2D lensing). In that case, the lens map takes the form

[TABLE]

where $\varepsilon,\boldsymbol{\lambda}_{k}$ are respectively the Einstein radius and angular position of the $k$th lens, and $N_{\text{l}}$ is the total number of lenses. The 2D-lensing choice was mostly made for simplicity, although it would not be difficult to generalize it to a 3D distribution of lenses, e.g., using multi-plane lensing Blandford and Narayan (1986). Reference Lee et al. (1997) found that, for a given optical depth, 3D lensing tends to give rise to fewer caustic mergers than 2D lensing. This suggests that the impact of lens-lens coupling is weaker in 3D than in 2D. Therefore, the departure from the predictions of the strongest-lens model that will be described here shall be seen as an upper-bound estimate.

In practice, we consider a square map with edge $\beta_{\text{map}}$ in the source plane. Given a number $N_{\text{l}}$ of lenses, we then choose their Einstein radius $\varepsilon$ such that the desired optical depth $\kappa=N_{\text{l}}\pi\varepsilon^{2}/\beta_{\text{map}}^{2}$ is reached. The lens positions $\boldsymbol{\lambda}_{k}$ are then randomly picked in such a way that their individual Einstein disk $\pi\varepsilon^{2}$ is entirely comprised in the map. An example of inverse ray shooting, with $N_{\text{l}}=10$ and $\kappa=0.5$ is given in Fig. 4. Note the important contraction from the image plane to the source plane, which requires the area of the image plane to be larger than the considered area in the source plane.

V.1.2 Adaptive mesh refinement

In order to get a high-resolution amplification map while maintaining good accuracy, it is necessary to generate a very large number of points in the source plane. Then, several options are available to evaluate $n_{\text{s}}(\boldsymbol{\beta})$ in the source plane. Kernel density estimation algorithms can be very accurate, but they are computationally too expensive. We chose to simply count the number of sources which end in each pixel of the source plane. However, since amplification (and hence $n_{\text{s}}$) varies over several orders of magnitude throughout the map, it is necessary to adapt the pixel size depending on its position.

We addressed this issue with a simple adaptive mesh refinement (AMR) procedure. In cosmology, AMR is particularly useful to resolve high-density regions in $N$-body codes, such as Ramses Teyssier (2002). Its principle is depicted in Fig. 5. Starting from a regular coarse grid, count the number of sources in each pixel; if the number of sources exceeds a given threshold $N_{\text{max}}$, divide the pixel in four equal subpixels, and repeat this operation with the subpixels until all of them contain less than $N_{\text{max}}$ sources. The density of sources in a subpixel of area $\Omega_{\text{sub}}$ containing $N_{\text{sub}}$ sources is then estimated as $n_{\text{s}}=N_{\text{sub}}/\Omega_{\text{sub}}$. Since $\Omega_{\text{sub}}$ can be very small, AMR allows one to locally access very high values of the amplification without the need to divide the whole map into tiny pixels. An example of amplification map obtained from this procedure is given in Fig. 6.

V.1.3 Performance

The code has been written in Python. In its current, nonoptimized and nonparallelized form, for $N_{\text{l}}=100$, $N=2\times 10^{7}$, it takes roughly 6 hours to generate a full map on a laptop with Intel Core i5 CPU. Computing time is linear in $N_{\text{l}}N$. A file containing all source positions typically occupies hundreds of MB to a few GB, depending on $N$. This amount of disk space is divided by 20 once the sources are arranged on a refined map, and their exact positions are deleted.

V.2 Mean amplification

As a first application, we consider the mean amplification $\left\langle A\right\rangle$ over the simulated map. Note that $\left\langle\cdots\right\rangle$ represents an average over randomly distributed sources. This source-averaging procedure must be distinguished from directional averaging, which would give the same weight to random directions in the image plane Kibble and Lieu (2005); Kaiser and Peacock (2016); Bonvin et al. (2015); Fleury et al. (2017). From a simulated amplification map, such as the one depicted in Fig. 6, it is straightforward to estimate $\left\langle A\right\rangle$ as follows:

[TABLE]

where the sum runs over each (sub-)pixel $i$ of the map, $\Omega_{i}$ being the area of the (sub-)pixel, and $N_{i}$ the sumber of sources in it. The areas $\Omega_{i}$ cancel out, so that $\left\langle A\right\rangle$ is simply the ratio between the total number of sources $N_{\text{s}}$ in the map, and the number of images which were in the same region before lens mapping.

The average amplification can also be estimated as the relative enhancement of the map’s area from the source plane to the image plane, $\left\langle A\right\rangle=\Omega_{\text{im}}/\Omega_{\text{s}}$. Let us approximate the map as a disk with radius $\beta_{\text{map}}$ instead of a square, and deal with the set of lenses as a single point lens, with squared Einstein radius $\Theta_{\text{E}}^{2}=N_{\text{l}}\varepsilon^{2}=\kappa\beta_{\text{map}}^{2}$. Then the principal image of the map’s contour is a circle with radius \theta_{\text{map}}=\big{(}\beta_{\text{map}}+\sqrt{\beta_{\text{map}}^{2}+4\Theta_{\text{E}}^{2}}\big{)}/2. Thus,

[TABLE]

Note the importance of the finite extent of region containing the lenses in the above estimate. This is what allows us to treat the $N_{\text{l}}$ lenses as a single point lens producing an image of the map’s edge. If the lens distribution was infinite in extent, the expected result would be $\left\langle A\right\rangle=(1-\kappa)^{-2}$.

Figure 7 shows the simulated mean amplification $\left\langle A\right\rangle$ as a function of optical depth $\kappa$, and compares it with the three models (strongest lens, multiplicative, and hybrid) considered in this article. In the low-optical depth regime, $\left\langle A\right\rangle\approx 1+2\kappa$, and the three models are in excellent agreement with the simulation. For larger optical depths, all three models fail. As expected, the strongest-lens approximation underestimates $\left\langle A\right\rangle$ while the multiplicative approach overestimates it. The hybrid model lies in between, but it still overestimates $\left\langle A\right\rangle$. However, since the estimate (100) provides an excellent fit to the simulation, the failure of $\left\langle A\right\rangle_{\text{s}},\left\langle A\right\rangle_{\text{m}},\left\langle A\right\rangle_{\text{h}}$ to predict the correct $\left\langle A\right\rangle$ must be partly attributed to finite-size effects. Indeed, all three models assume an infinite number of lenses in an infinite Universe, whereas the numerical setup is finite.

V.3 Comparing PDF models

Let us now confront the three PDF models $p_{\text{s}}(A)$, $p_{\text{m}}(A)$, $p_{\text{h}}(A)$ to the simulations. For the sake of completeness, we add a fourth one, empirically proposed by Rauch in 1991 Rauch (1991). For low optical depths, the following expression was found to provide a good fit to Monte-Carlo simulations, aiming to determine microlensing amplification statistics in an expanding Universe filled with point lenses,

[TABLE]

In Eq. (101), $\kappa_{\text{eff}}$ and $b$ are two parameters fixed by the conditions that $p_{\text{R}}$ is normalized to unity, and that $\left\langle A\right\rangle_{\text{R}}$ gives the correct mean amplification (as given by the simulation). Rauch’s fitting formula, combined with flux conservation, has been used in Refs. Seljak and Holz (1999); Metcalf and Silk (1999, 2007), and even more recently in Ref. Zumalacárregui and Seljak (2018) to set constraints on the abundance of primordial black holes from supernova lensing (see also Ref. García-Bellido et al. (2018)).

We consider four values of the optical depth, $\kappa=0.01,0.1,0.5,1$. For each value, the parameters of the simulation, Rauch’s fitting function, and the mean amplification, are summarized in Table 1. The four numerical amplification PDFs are shown together in Fig. 8, and individually compared to the analytic models in Figs. 9, 10, 11, 12.

As expected, for low optical depths $(\kappa=0.01,0.1)$, all four models $p_{\text{s}},p_{\text{m}},p_{\text{h}},p_{\text{R}}$ essentially coincide, and are in very good agreement with numerical results. However, in the optically-thick regime ($\kappa=0.5,1$), all four models heavily fail to reproduce the actual behavior of $p(A)$. Specifically, the models tend to underestimate the probability of large amplifications ($A\sim 10$), and overestimate probability of very large amplifications ($A>100$).

The reason for this failure is nonlinear lens-lens coupling, which none of the four models really accounts for. The impact of lens-lens coupling is clearly visible on the amplification maps. As $\kappa$ increases, the map’s aspect changes from a set of small isolated regions, where the amplification can become extremely large, to an intricate and cuspy caustic network, where large areas are characterized by intermediate amplifications, but where it is difficult to access very high values of $A$. As already observed in the literature Nityananda and Ostriker (1984); Schneider (1987a); Mao (1992); Kofman et al. (1997), this explains the leaking from very large to large amplifications, compared to what would be naïvely expected from approaches where the lenses are independent.

Of the four analytic models, the multiplicative approach is probably the worst, while the strongest-lens approach may be considered the least bad. However, in order to accurately model the amplification PDF in the optically-thick regime, it is necessary to properly tackle the problem of lens-lens coupling. Such a program, already initiated by other authors Schneider (1987a); Kofman et al. (1997); Lee et al. (1997), is beyond the scope of the present article, but shall be addressed in a future work.

V.4 Extended sources

Let us finally consider the impact of the finite size of sources. In Sec. IV, this effect was discussed for homogeneous disk-sources with angular radius $\sigma$. In the present numerical setup, however, square sources are more easily implemented. In order to connect the former theoretical results to the latter numerical ones, we assume that a disk source with angular area $\pi\sigma^{2}$ is mostly equivalent to a square source with the same area.

As mentioned in Sec. IV, the microlensing amplification map for extended sources is a smoothed version of the the map obtained with point sources. This property is clearly illustrated in Fig. 13 for $\kappa=0.5$. Four different source sizes are represented, $r=0,10^{-2},10^{-1},1$, where $r=\sigma/\varepsilon$ is the reduced size of the source. We are using square sources with edge $\ell=\sqrt{\pi}\sigma$, so that $\ell^{2}=\pi\sigma^{2}=\pi r^{2}\varepsilon^{2}$. As the source’s radius increases, the map gets smoother; in particular, the maximum amplification decreases. The impact of that smoothing on the amplification PDF is depicted in Fig. 14.

What is the performance of the strongest-lens and multiplicative models to describe $p(A)$? For low optical depths, we know from Eq. (94) that both models coincide up to terms of order $\kappa^{2}$. Their common prediction is in excellent agreement with the numerical results, as we can see in Fig. 15 for $\kappa=10^{-2}$. However, for higher optical depths, both models fail just like in the point-source case, see Fig. 16. This is especially true for $r=1$, where the strongest-lens model underestimates the maximum amplification by almost a factor $5$. This is due to the fact that, in the strongest-lens model, the light beam can only enclose one lens, while for high optical depths the beam would typically enclose multiple lenses. The multiplicative model does not predict any maximum amplification, thereby missing an important property of the actual $p(A)$.

VI Conclusion

This article reviewed simple analytic models for the probability distribution of microlensing amplifications due to a set of microlenses. All the models are parameterized by a single quantity, namely the optical depth $\kappa=\Sigma\pi\langle r_{\text{E}}^{2}\rangle$, which depends on the surface density of lenses $\Sigma$ and their mean Einstein radius $r_{\text{E}}$. The optical depth quantifies the expected amount of lensing caused by the system. It also coincides with the weak-lensing convergence which would be observed if the lenses’ mass was smoothly distributed.

In the strongest-lens model, the amplification $A=A_{\text{s}}$ is only due to the lens whose dimensionless impact parameter with respect to the line of sight is the smallest. In order to allow for the combined effect of many lenses, we considered a multiplicative model. In that approach, the net amplification is assumed to be the product of every individual lens amplification, $A=A_{1}A_{2}\ldots A_{N}$. Since the existing derivation of the amplification PDF in that context was incorrect, we proposed an alternative one, and found the same final result. Finally, we considered a hybrid model where the net amplification is the product of the strongest individual amplification with the mean principal amplification of the other lenses, $A=A_{\text{s}}\bar{A}^{+}_{\text{w}}$. While calculations were mostly conducted assuming point sources, we also derived the finite-source corrections to the strongest-lens and multiplicative models.

These three simple models were confronted to two-dimensional inverse-ray-shooting simulations. For low optical depths, $\kappa<0.1$, they essentially coincide and are in excellent agreement with numerical results. They also predict the correct behavior of the mean amplification, $\left\langle A\right\rangle=1+2\kappa+\mathcal{O}(\kappa^{2})$. However, for large optical depths, all three models fail to reproduce both the mean amplification and the main features of its PDF. In particular, the probability of large amplifications ($A\sim 10$) is underestimated, while the probability of very large amplifications ($A>100$) is overestimated. While the misestimate of the mean amplification can easily be attributed to the finite size of the map, the mismatch of the full PDF comes from a deeper modeling limitation. Namely, none of the models properly account for lens-lens coupling, which cannot be neglected in the optically-thick regime.

An accurate description of the microlensing amplification PDF is essential for using this observable as a reliable test of the nature and small-scale distribution of dark matter. This is particularly true for large optical depths, where most of the signal is expected. In future work, we aim to efficiently model the impact of lens-lens coupling, in order to produce realistic amplification statistics in a clumpy Universe.

Acknowledgements.

PF thanks Jonathan Blazek and David Harvey for discussions and encouragements. We also thank the anonymous referee of Physics of the Dark Universe, for their constructive comments that significantly enriched the content of this article. PF acknowledges the financial support of the Swiss National Science Foundation. PF received the support of a fellowship from “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/PI19/11690018. JGB thanks the CERN TH-Division for hospitality during his sabbatical, when this project was initiated, and acknowledges support from the Research Project FPA2015-68048-03-3P [MINECO-FEDER] and the Centro de Excelencia Severo Ochoa Program SEV-2012-0597. He also acknowledges support from the Salvador de Madariaga Program, Ref. PRX17/00056.

Appendix A Pei’s erroneous derivation of $p_{\rm m}(A)$

This appendix reproduces Pei’s derivation of the multiplicative amplification PDF $p_{\text{m}}(A)$, where some mistakes are corrected, and the remaining weaknesses are pointed out.

A.1 Setup and notation

Consider a source located at redshift $z_{\text{s}}$, and randomly distributed lenses with redshifts $z\in[0,z_{\text{s}}]$, mean number density $n(z)$, and possibly different masses. We introduce the following notation:

• •

$p(A,z)\,\mathrm{d}A$ is the probability that the lenses within $[0,z]$ yield a total amplification between $A$ and $A+\mathrm{d}A$. Thus, $p(A)=p(A,z_{\text{s}})$.

• •

$q(A,\mathrm{d}z)\,\mathrm{d}A$ is the probability that the lenses within $[z,z+\mathrm{d}z]$ yield a total amplification between $A$ and $A+\mathrm{d}A$.

• •

$f(A,z)\,\mathrm{d}A\,\mathrm{d}z$ is the number of lenses777This quantity corresponds to $\rho(A,z|z_{\text{s}})$ in the notations of P93. between $z$ and $z+\mathrm{d}z$ generating an individual amplification between $A$ and $A+\mathrm{d}A$.

It is tempting to identify $q$ with $\partial p/\partial z$. However, this cannot be true, since $q$ is a probability density, hence normalized to 1, while the integral of $\partial p/\partial z$ over $A$ vanishes.

A.2 Integrodifferential equation for $p(A,z)$

The actual relationship between $p$ and $q$ comes from the prescription (38) for the multiplicative combination between two groups of lenses. Namely, splitting the lenses of $[0,z+\mathrm{d}z]$ between $[0,z]$ and $[z,z+\mathrm{d}z]$, we find

[TABLE]

The tricky part then consists in relating $q$ to $f$. One can arguably decompose $q$ as follows,

[TABLE]

where $q_{k}(A,\mathrm{d}z)$ is the probability that exactly $k$ lenses within $[z,z+\mathrm{d}z]$ significantly contribute to the amplification $A$. As $\mathrm{d}z\rightarrow 0$, it is rather intuitive that $q_{k}=\mathcal{O}(\mathrm{d}z^{k})$; thus, at first order in $\mathrm{d}z$, we must have

[TABLE]

where $C$ is a constant to be determined. Indeed, if no lens contributes ($k=0$), then $A=0$, so that $q_{0}(A,\mathrm{d}z)\propto\delta(A-1)$. The identification between $q_{1}$ and $f\mathrm{d}z$ then comes from the definition888P93 erroneously considers $q_{1}(A,\mathrm{d}z)=\rho(A,z|z+\mathrm{d}z)$ of $f$. We determine $C$ using the normalization of $q$,

[TABLE]

Note that, since $\bar{f}$ represents the total number of lenses per unit redshift, this quantity is generally infinite! However, its combination with $\delta(A-1)$ and $f(A,z)$ should remain finite. Substituting the expression of $q$ in Eq. (102), we get the following integrodifferential equation:999Equation (107) is corrected with respect to Eq. (4) of P93, which includes confusions between $z$ and $z_{\text{s}}$, as well as a mathematical mistake about derivatives of dependent integrals.

[TABLE]

Alternatively, one can use logarithmic probabilities, $P(L)=Ap(A)$, $F(L)=Af(A)$, with $L=\ln A$, to get

[TABLE]

A.3 Solving the equation

The main difficulty of Eq. (108) is the convolution product on its right-hand side. This is greatly simplified in Fourier space,

[TABLE]

where we used that

[TABLE]

Equation (109) is easily integrated, and the result translated in terms of $p$, yielding

[TABLE]

which is Eq. (8) of P93.

A.4 Determining $\tilde{F}(K,z)$

The last step of the derivation consists in explicitly computing $\tilde{F}(K,z)$. For that purpose, we first determine $f(A,z)$, which by definition reads

[TABLE]

For point lenses, we have seen in Sec. II.1 that the amplification only depends on the reduced impact parameter $u=b/r_{\text{E}}$. Hence, for one lens, the region of the plane $z=\mathrm{cst}$ such that its amplification lies in $[A,A+\mathrm{d}A]$, has an area

[TABLE]

The quantity $\mathrm{d}\sigma_{A}/\mathrm{d}A$ must be understood as a differential cross section. Besides, from Eq. (6), one finds

[TABLE]

Suppose that an infinity of such lenses are randomly distributed within $[z,z+\mathrm{d}z]$, with density $\mathrm{d}\Sigma/\mathrm{d}z=\mathrm{d}^{3}N/\mathrm{d}^{2}S\mathrm{d}z$. The average number of lenses within $\mathrm{d}z$ causing an amplification within $\mathrm{d}A$ is $\mathrm{d}^{2}N=\mathrm{d}\Sigma\,\mathrm{d}\sigma_{A}$, whence

[TABLE]

Here, we have implicitly assumed that all the lenses have the same Einstein radius $r_{\text{E}}$, i.e. the same mass. However, it is straightforward to generalize the above rationale for lenses with a spectrum of masses. The idea consists in replacing $\mathrm{d}\Sigma/\mathrm{d}z$ with $\mathrm{d}^{2}\Sigma/\mathrm{d}z\mathrm{d}m$, the density of lenses per unit mass, and then integrating over $m$.

The Fourier transform

[TABLE]

does not exist, because the integral does not converge at the limit $A=1$. However, the combination $\tilde{F}(K,z)-\tilde{F}(0,z)$, as it appears in Eq. (111) does exist. Precisely,

[TABLE]

Therefore, the final result is

[TABLE]

with

[TABLE]

Equation (119) is the main result of P93. It is remarkable that it perfectly matches the outcome of Sec. III.3, despite the several issues of the original derivation.