{'fraction_non_alphanumeric': 0.04343212163309408, 'fraction_numerical': 0.019313884386817935, 'mean_word_length': 4.302741968487825, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New Approximation of Magnification Statistics for Random Microlensing of Magnified Sources

Gravitationally lensed extragalactic sources are often subject to statistical microlensing by stars in the galaxy or cluster lens. Accurate models of the flux statistics are required for inferring source and lens properties from flux observations. We derive an accurate semi-analytic approximation for calculating the mean and variance of the magnification factor, which are applicable to Gaussian source profiles and arbitrary non-uniform macro lens models, and hence can save the need to perform expensive numerical simulations. The results are given as single and double lens-plane integrals with simple, non-oscillatory integrands, and hence are fast computable using common Monte Carlo integrators. Employing numerical ray-shooting experiments, we examine the case of a highly magnified source near a macro fold caustic, and demonstrate the excellent accuracy of this semi-analytic approximation in the regime of multiple micro images. Additionally, we point out how the maximum persistent magnification achievable near a macro caustic is fundamentally limited by the masses and number density of the foreground microlenses, in addition to the source's physical size.


INTRODUCTION
Strong gravitational lensing of extragalactic sources by intervening galaxies or galaxy clusters are exquisite astrophysical probes. The observed fluxes of these sources are often strongly affected by random microlensing, collectively caused by many compact masses embedded within the foreground lens (Kayser et al. 1986;Wambsganss 1992). Observations of multiply imaged quasars first engendered theoretical interests in statistical microlensing (see Wambsganss (2006) for a review), with the major aim to resolve the size and structure of the quasar accretion disk (Blackburne et al. 2011) and to probe compact lenses of planetary to stellar masses.
To quantify the flux variability under random microlensing, a statistical approach is appropriate. For observations carried out at random epochs, the key observables include the mean and variance of the magnification factor. The scenario of microlenses embedded in constant background convergence and shear has been extensively studied, mostly in the context of multiply-imaged quasars. In this case, the magnification factor averaged over random microlens realizations equals the macro one regardless of the source size, while the variance is difficult to compute analytically and was the focus of many past works. Lenses solely comprised of compact objects were studied in the pioneering works of Deguchi & Watson (1987) and Refsdal & Stabell (1991), and the effects of a shear and a diffuse surface mass component were investigated in follow-up studies Refsdal & Stabell 1997). The statistical formalism was further developed in more recent works (Neindorf 2003;Tuntsov et al. 2004;Goodman & Sun 2014), and other methods have been explored (Fleury & García-Bellido 2020). When analytical results are unavailable, unreliable, or cumbersome to apply, direct numerical simulations come to rescue. Important numerical techniques that greatly enhance efficiency and accuracy include inverse ray-shooting (Kayser et al. 1986), the hierarchical tree algorithm (Wambsganss 1999), and the image tracking method (Lewis et al. 1993;Witt 1993).
Recently, theoretical interests in statistical microlensing have been reinvigorated by the detections of high-z stars (Miralda-Escudé 1991) magnified by a spectacular ∼ 10 2 -10 3 folds near critical curves of cluster lenses Rodney et al. 2018;Chen et al. 2019;Kaurov et al. 2019). Owing to the extreme macro magnification factors realized in these cases, even a low convergence of intracluster stars around the cluster Einstein radius κ ∼ 10 −3 -10 −2 (compared to κ ∼ 0.1-1 around galaxy Einstein radii) leads to frequent microlensing brightening episodes (Venumadhav et al. 2017;Diego et al. 2018;Oguri et al. 2018;Diego 2019), making these outstanding probes of compact constituents of the lens mass. Such stochastic microlensing may act also on a cluster of stars if the cluster as a whole is highly magnified Dai 2021). With new observational data, other highly magnified candidate sources: lensed quasar images suspected of large magnification (Fujimoto et al. 2020) or flux anomaly (Glikman et al. 2018) have been reported; a magnified "knot" in a Cosmic Noon starburst showing perplexing flux anomalies (Vanzella et al. 2020).
While semi-analytic scaling laws have offered much insight into the behavior of statistical microlensing in the extreme magnification regime, more accurate modeling of flux statistics have so far relied on large yet artful numerical simulations (Venumadhav et al. 2017;Diego et al. 2018;Diego 2019;Dai 2021). Sometimes, it may even be computationally prohibitive to simulate for realistic parameter values! In this work, with a focus on the macro caustic vicinity, we aim to develop a general and practical analytic model, which will both deepen our understanding of microlensing in the high optical depth regime and facilitate scans of large parameter space.
Microlensing has a significant impact on the magnifications achievable near a macro caustic (Venumadhav et al. 2017;Diego et al. 2018;Diego 2019). Microlenses of a characteristic Einstein radius θ induce a stochastic deflection component κ 1/2 θ in the ray equation (Katz et al. 1986), up to a multiplicative "Coulomb" logarithm reflecting the long-ranged nature of point lens deflections. This sets an effective smoothing scale on the source plane, independently of the source's angular extent σ W . If the magnification caused by the macro lens is uniform across this scale, stochastic microlensing conserves the mean magnification, but only induces fluctuations around it. However, when the macro magnification varies significantly, a situation that inevitably arises near a caustic, even the mean magnification is modified. Indeed, Venumadhav et al. (2017) showed that a sharp caustic induced by a macro lens of a smooth mass profile is "smeared out" across a width of ∼ κ 1/2 θ by microlenses.
In the absence of microlenses, the maximum magnification is realized when a finite source grazes a sharp marco caustic, µ max ∼ (σ W d) −1/2 , where d is typically the inverse of the characteristic angular scale of the macro lens that produces the caustic. Sub-galactic substructure lenses produce secondary caustics (Dai et al. 2018) that have larger values of d compared to what smooth galaxy-scale or clusterscale lenses can produce, and hence decreased values of µ max . In the presence of microlenses, the sharp macro caustic is disrupted, and a corrugated network of micro caustics form instead. In this situation, the random deflection scale κ 1/2 θ comes into play. Now the mean magnification can reach a maximum value ∼ (σ eff d) −1/2 where σ eff σ 2 W + κ θ 2 , i.e. the grazing scale can get as small as σ W or κ 1/2 θ , whichever is larger.
In real astrophysical contexts, a source may be sufficiently large (e.g. quasar, SN) to overlap multiple micro caustics, i.e. σ W θ κ −1/2 /(σ eff d) −1/2 , when the density of micro caustics is the highest in the proximity of a macro caustic. As a result, its flux fluctuations around the mean do not exceed the mean by any large factor, and hence (σ eff d) −1/2 is often a fair estimate for the maximal possible magnification after accounting for fluctuations. For smaller sources, it may still be true that σ eff θ κ −1/2 /(σ eff d) −1/2 , so that flux fluctuations remain mild as many disconnected micro images form and contribute uncorrelated magnification fluctuations. Even in this regime, (σ eff d) −1/2 is not an underestimate of the maximal achievable magnification. Flux fluctuations can greatly exceed the mean value only for very small sources (e.g. individual stellar photospheres), and for a sufficiently low surface number density of microlenses. During these short events of micro caustic crossing, the peak magnification is dominated by just a pair of micro images. A conceptual obstacle to fully analytic calculations of the microlensing flux statistics has to do with the "ultraviolet" (UV) and "infrared" (IR) divergences, which are terms borrowed from field theory. These logarithmic divergences arise because the deflection due to any single point microlens scales inversely with the impact parameter, analogous to the Coulomb divergence associated with the inverse-square force law in three-dimensonal space. The UV divergence is traced back to arbitrarily large ray deflections at arbitrarily small impact parameters to any single point microlens; consequently, the probability distribution functions (PDFs) for the random deflections have divergent higher-order moments, and the corresponding characteristic functions (CFs) contain non-analytic logarithmic terms. We seek a prescription in which the UV divergence is regulated by the finite source size. On the other hand, the IR divergence implies that the statistics can be sensitive to far-away microlenses distributed over the largest scales on the lens plane. We require a prescription in which the dependence on the IR cutoff on statistical averages is manifest and unambiguous.
The main result of this work is an accurate semi-analytic approximation for the mean and variance of the magnifications (Eq. (29) and Eq. (37) respectively), for a Gaussian source profile, and in the regime that multiple micro images form. The gist of the approximation is about capturing the Gaussian bulk of the deflection distribution, while judiciously neglecting the non-Gaussian tail of large deflections which are only important for either a small portion of the source or a minority of the micro images. The approximation determines the UV and IR logarithms in a physical way, and respect the expected translational and rotational symmetries of the deflection distributions. The results only require evaluating two-and four-dimensional integrals with simple, well-behaved integrands, and are applicable to any macro lens model. In practice, answers can be obtained in less than a second by employing standard Monte Carlo integrators. The new results can be directly used to efficiently and accurately quantify the magnification statistics, for microlenses embeded in a variety of macro lens models and as a function of microlens abundance and source size, thus saving the need to perform expensive yet tricky numerical simulations. Conversely, the semi-analytic answers enable calibration to numerical codes.
The remainder of the paper is organized as follows. In Section 2, we introduce the general theoretical framework to compute magnification statistics, in particular, the mean and the variance. In Section 3, we examine a toy model in which microlensing deflections are Gaussian random variables. We will develop useful intuition into the problem as this toy model is exactly solvable. In Section 4, we turn to the real problem of discrete microlenses, and derive our key results. We then demonstrate that our semi-analytic approximation agrees well with numerical ray-shooting experiments, for various parameter choices ranging from Gaussian to non-Gaussian flux variability behaviors. We will discuss our results in Section 5, before we give concluding remarks in Section 6. Additional technical details are presented in Appendices for reference. Results of Section 3 and Section 4 are presented in dimensionless angular units, and can be easily scaled to the appropriate physical units in any specified astrophysical context.

THEORETICAL FRAMEWORK
A dozen of previous studies have treated microlensing using a statistical theory. Those include the earlier works of Deguchi & Watson (1987), Katz et al. (1986),  and , as well as a more general formulation in Neindorf (2003). For clarity, we reintroduce this theoretical framework, which is based on the concept of multivariate probability distribution functions and the corresponding characteristic functions.
We decompose the deflection field into a background component α B (x) and a fluctuating component α(x). The former varies smoothly as a function of the image plane position x and is given. The latter, being stochastic in nature with specific spatial correlations, will be given a statistical treatment. We adopt the general assumption that the α(x)'s have correlation functions that are invariant under spatial translations and rotations, as is the case if α(x) is generated by point-like microlenses that are uniformly distributed on the lens plane.
For a point source at the source-plane position y, the lens equation is y = x − α B (x) − α(x). Let W (y) be the normalized surface brightness profile of a finite-sized source centered at the source-plane origin, satisfying d 2 y W (y) = 1. If unresolved, the total magnification summed over all geometric images can be written as Throughout, we use the · · · notation to indicate averaging over random realizations of α(x). Inserting the Fourier decomposition of W (y), we obtain the mean magnification factor (Venumadhav et al. 2017) where is the Fourier wave vector conjugate to the real-space angular variable, and W ( ) = d 2 y W (y) e i ·y is the Fourier transform of W (y). This expression depends on the characteristic function (CF) involving the stochastic deflection at one image-plane position x. By statistical homogeneity, e i ·α(x) is independent of x. The result can be recast into the form where P W 1 (α) is the one-point probability distribution function (PDF) for the random deflection α(x) at any x, P 1 (α) (see e.g. Katz et al. (1986)), convoluted with the source profile By transforming the integration variable from y to α = x − y − α B (x), Eq.
(3) can also be written as (see e.g. Venumadhav et al. (2017)) where µ B (y) is the magnification of a point source at y due to only the background deflection. Eq. (5) can be interpreted as smoothing the point-source background magnification pattern with a "point spread function" P W 1 (α), whose characteristic width is set by that of the source profile W (y) or that of the random deflections P 1 (α), whichever is larger. µ W (y) equals the background magnification µ B (y) if the latter is approximately uniform over the smoothing scale set by P W 1 (α). A particularly interesting situation arises near a lensing caustic where µ B (y) often varies rapidly; as a result µ W (y) can differ substantially from µ B (y).
The simplest way to quantify the scatter in the value of µ W (y) due to random microlens realizations is the second moment, which is given by (Neindorf 2003) The CF involving the stochastic deflections at two different image-plane positions x 1 and x 2 needs to be computed. By statistical homogeneity, e i 1 ·α(x 1 ) e i 2 ·α(x 2 ) only depends on x 2 − x 1 , but not on x 1 and x 2 individually. The magnification factor has a standard deviation Std[µ W (y)] = µ W (y) 2 − µ W (y) 2 (Neindorf 2003).

TOY MODEL: GAUSSIAN RANDOM DEFLECTIONS
Before we consider discrete, point microlenses as realistic random deflectors, we would like to first solve a toy model in which α(x) behaves strictly as a Gaussian random vector field on the image plane. The toy model allows for the analytic calculation of many results, and hence will offer us much insight into how the source-convoluted magnification factor behaves in the presence of stochastic deflections. We note that this Gaussian deflection model precisely describes the collective lensing effects of axion minihalos on lensed extragalatic stars crossing micro caustics (Dai & Miralda-Escudé 2020).
For full analytic tractability, we choose to consider a Gaussian source profile with half width σ W , W (y) = [(2π) σ 2 W ] −1 exp(−(1/2) |y| 2 /σ 2 W ) throughout. While many commonly studied sources such as individual stars are more appropriately modeled as uniform disks, the assumption of a Gaussian source profile does not impose a fundamental limitation of our derivation, as we will see later. The Fourier transform of this Gaussian profile is given by We can generally assume that α(x) is the gradient of a scalar potential (valid for deflection by a single lens plane), α(x) = ∇ ψ(x), and that in the Fourier domain the scalar potential has an isotropic power spectrum ψ( ) ψ * ( ) = (2π) 2 δ D ( − ) P ψ ( ). The two-point correlation for α(x) can be decomposed into a longitudinal component C (r) and a transverse component C ⊥ (r): where r 12 = x 2 − x 1 and r 12 = |r 12 |. These are related to the potential power spectrum through: For a well-behaved P ψ ( ), we can define, at zero separation r = 0, For Gaussian random deflections, the mean magnification factor is given by an integral over a single image plane, The second moment is given by an integral over double image planes, To condense the expression, we have constructed a four-component vector We also introduce the four-by-four identity matrix I, and define a four-by-four covariance matrix where the shorthand notations c := cos ϕ and s := sin ϕ are used, under the polar coordinate parametrization r = r [cos ϕ, sin ϕ]. Given the matrix determinant det[C(r)] = (C(0) + C (r)) (C(0) − C (r)) (C(0) + C ⊥ (r)) (C(0) − C ⊥ (r)), C(r) is a positive-definite matrix if and only if C(0) ± C (r) > 0 and C(0) ± C ⊥ (r) > 0. In our toy model, these are strictly guaranteed by Eq. (9) and Eq. (10) provided that P ψ ( ) > 0. If the off-diagonal elements in Eq. (14) were all vanishing, we would have concluded that µ W (y) 2 = µ W (y) 2 . Hence, the off-diagonal matrix elements of C(r) is the reason for the nonzero variance for µ W (y). For giving numerical examples, we specify a simple analytic form for the potential power spectrum: The resultant root-mean-square (RMS) deflection is |α(x)| 2 = √ 2 α 0 . The parameter σ ψ is introduced to set the coherent scale of the deflections on the image plane. The deflection correlation functions can be explicitly computed, To study the situation of a rapid varying µ B (y), let us consider the background deflection of a fold caustic, which is the most commonly encountered caustic. We use the following parameterization for the fold caustic (Schneider et al. 1992), The parameters are the local background convergence κ 0 , and a gradient vector d = [d 1 , d 2 ] which is related to the third-order derivative of the background lensing potential. We introduce d = d 2 1 + d 2 2 and α = − tan −1 (d 1 /d 2 ) following the notation of Venumadhav et al. (2017). The coordinate system is conveniently chosen so that the caustic aligns with the y 2 axis on the source plane, and the corresponding critical curve on the image plane intersects the x 1 axis at an angle α. As an example, we consider the case d 1 > 0 and d 2 = 0, so that d = d 1 and α = π/2. Several important scales can be identified on the source plane. First, the characteristic source size is σ W . Second, in the vicinity of the ideally smooth background caustic, the random deflections cannot be neglected. This is relevant within a source-plane width ∼ α 0 from the background caustic, as set by the root mean square (RMS) deflection. Within this proximity of the background caustic, the background magnification reaches µ f ∼ (α 0 d) −1/2 , which we assume to be large. When µ f is greater than the inverse of the typical convergence fluctuation ∼ (α 0 /σ ψ ), micro caustics join to form a network within this narrow band. The random deflections have a coherent scale ∼ σ ψ on the image plane. When ray-traced onto the source plane following the background ray equation, this on average maps to a compressed scale ∼ σ ψ /µ f = σ ψ (α 0 d) 1/2 for the corrugated micro caustic pattern, along the direction perpendicular to the background caustic. When the micro caustic network does arise, parametrically we have the hierarchy α 0 > σ ψ (α 0 d) 1/2 . Figure 1 shows a numerical example of how a smooth macro caustic is replaced by a corrugated micro caustic pattern due to the Gaussian random deflections.
We evaluate the multi-dimensional integrations Eq. (11) and Eq. (12) using the widely used Monte Carlo algorithm vegas (Lepage 1978(Lepage , 2020. In Figure 2, we present numerical examples that verify the intuitions we have developed based on our analytic scrutiny of the Gaussian random deflection model. We see that the mean magnification µ W (y) as a function of source center y indeed has a smoothed behavior at the background caustic. The maximal mean magnification and the width of the smoothed curve is set by the size of random deflections α 0 if that exceeds the source size α 0 > σ W . From one realization to another, µ W (y) fluctuates around the mean µ W (y) . However, such fluctuations are suppressed when the source size σ W is larger than the characteristic separation of micro caustics.

RANDOM DEFLECTIONS FROM COMPACT MICROLENSES
In most astrophysical situations, random deflections are due to compact lenses such as individual stars. Figure 3 shows an example of the corrugated micro caustic network cast by random point lenses, which looks qualitatively different from that formed due to Gaussian random deflections.
An important difference between the case of point microlenses and that of Gaussian random deflections is that very large deflections are generated with a low but non-negligible probability. This happens when the ray encounters a microlens at a small impact parameter. This not only greatly enhances the flux fluctuations, but also renders the distribution of α(x) heavy-tailed (Katz et al. 1986). As a consequence, α(x) has divergent second-order and higher-order moments, which hinders exact analytic treatment.
For a simple discussion, we assume that all microlenses have the same mass, and hence the same angular Einstein scale θ . We will comment on the generalization to the case of a continuous microlens mass distribution toward the end of the this Section. Other angular scales of the problem can all be expressed in units of θ . We assume that the microlenses are uniformly distributed on the lens plane and contribute a mean convergence κ > 0. . Total magnification factor µ W for a finite source with a Gaussian profile as a function of distance to the background caustic (vertical dotted line at y 1 = 0). We set parameter values κ 0 = 0.7, d = 10 −5 , α 0 = 0.2 and σ ψ = 10, as defined in the text, and consider increasing source sizes σ W = 0.005, 0.05, 0.2, 0.5 (from top to bottom). We show both the mean magnification µ W (red solid) and its standard deviation µ 2 W − µ W 2 due to random deflections (cyan band; ±1 σ). The magnification factor for an idealized point source without random deflections is also shown for comparison (blue dashed). For σ W < α 0 , the height and width of the peak µ W are set by α 0 , while the magnification fluctuations decreases as σ W increases. For σ W > α 0 , the height and width of the peak µ W are instead set by σ W , and the magnification fluctuations are highly smoothed out.
The microlensing deflection can be written as In the second term, we sum over contributions from all microlenses I = 1, 2, · · · , N , which are located at z I 's, respectively. In the first term, we subtract the mean deflection due to a uniform mass sheet of convergence κ . Including the first term enforces that α ml (x) = 0. Unlike in the case of Gaussian random deflections, analytic calculations of the CFs for α ml (x) are difficult for the    21). We compare the exact numerical evaluation (black) and the leading-order approximation (red), −(1/2 t 2 ) (1 − γ E + ln 2t), in the (1/t)-expansion. The leading-order approximation works remarkably well for t 1.
case of point microlenses. Still, the CFs can be reduced to those corresponding to a single microlens, provided that microlenses have independent positions on the lens plane. This forms the basis of the following analytic treatment.

One-point deflection statistics
Let P 1 [α ml (x)] be the PDF for α ml (x). Its Fourier transform gives the one-point CF: Here ϕ 1 [ ; x] is the one-point CF due to a single microlens. We specify that each microlens is uniformly distributed within a disk of radius R, and compute for α ml (0). We perform the following integral (Figure 9; left) where = | |. The integral I(t) can be expressed as a series expansion: where J 0 (x) is the Bessel function of the first kind, and γ E ≈ 0.577216 is the Euler-Mascheroni constant. For sufficiently large R, t 1, and we observe that the series rapidly converges. Keeping only the O(1/t 2 ) term gives a decent approximation ( Figure 4): The presence of the -dependent logarithm ln(2 R/θ 2 ) reflects that the deflection due to one microlens follows a heavy-tailed distribution; a Taylor expansion in around = 0 is not possible. Using κ = N θ 2 /R 2 , and setting x = 0, we derive Note that in the large N limit, the higher order terms in Eq. (21), i.e. ∼ t −2 k for k = 2, 3, · · · , are suppressed by ∼ 1/N k−1 , respectively. The variance of the deflection angle |α ml (0)| 2 ≈ κ θ 2 (1 − γ E + ln(2 R/θ 2 )) is always of order σ 2 eff = θ 2 κ . Additionally, there is a multiplicative factor that depends logarithmically on the IR cutoff scale R as well as the correlation scale of interest ∼ 1/ . By fixing R = R * and = * in the logarithm for some appropriate values of R * and * , we essentially approximate P 1 [α ml ] as a two-dimensional Gaussian distribution: where we introduce The key question is what the appropriate values for R * and * would be under this approximation. The R dependence originates from the collective influence of faraway microlenses. One straightforward choice for R * would be to use the full angular scale on the image plane over which the microlens convergence has nearly a uniform value κ . This scale can be orders of magnitude larger than the Einstein scale θ . For galaxy or intracluster microlensing, it may be taken to be the characteristic extent of the galactic stellar halo or the intracluster light halo, respectively. Katz et al. (1986) (hereafter KBP86) instead propose that R * should be the scale over which the microlenses produce incoherent deflections across the distribution of micro images, arguing that coherent deflections, for any fixed microlens realization, lead to an overall uniform shift of the micro images without affecting the magnification. Define an effective source size which accounts for intrinsic source size and effective "broadening" due to random microlensing deflections. Let µ B be the macro magnification near a background fold caustic as introduced in Section 3. The maximal value of µ B is limited by the effective source size σ eff , reaching µ f ∼ (σ eff d) −1/2 . Either macro image thus has an extent ∼ µ B σ eff . Following the choice of KBP86, we may set R * to µ B σ eff µ f σ eff = (σ eff /d) 1/2 . However, the choice of KBP86 is not entirely justifiable if the background lens model is nonuniform. The issue can be particularly non-trivial near a macro caustic where the (small) gradients of background convergence and shear set the strength of the caustic -large-scale Poisson fluctuations in the microlens number not only generate a coherent deflection, but also contribute small gradients of convergence and shear which can modify the effective caustic strength from realization to realization. In this case, we set R * to be the largest scale over which the microlenses are uniformly distributed.
The typical value of * to be set in the logarithm can be related to the inverse of the width of P W 1 [α], which is the amount needed to compensate for x − y − α B (x) = 0. This means that we can set * = 1/σ eff .
For the choice R * = µ f σ eff (KBP86) and * = 1/σ eff , the truncation of expansion Eq. (21) is valid if If σ eff is larger than the characteristic scale of the micro-caustic network ∼ θ κ −1/2 µ −1 B , then the first factor is smaller than unit. This corresponds to the situation where either the physical source extent σ W overlaps multiple micro caustics, or the typical microlensing broadening θ κ 1/2 overlaps multiple micro caustics (the regime of many micro images). The second factor is no larger than unity. Hence, θ 2 * /R * < 1 is satisfied in this situation. As we verify numerically in Figure 5, assuming the random microlensing deflections to follow approximately a two-dimensional Gaussian distribution can fairly accurately reproduce the exact source-profile convoluted distribution P W 1 (α ml ), provided that θ 2 * < R * . Eq. (25) can be recast into the following form where N * = κ R 2 * /θ 2 is the number of microlenses within an image-plane disk of radius R * . The first logarithmic term has been previously derived by KBP86 as ln(3.05 N 1/2 * ). Our result includes a second term ln(θ κ 1/2 * ) not found in KBP86, as we have argued that it is appropriate to set  . Source-profile convoluted distribution P W 1 (α ml ) = P W 1 (|α ml |). We consider the case in which the highest possible mean magnification µ W (y) is realized when the source center y is within a distance ∼ σ eff (see Eq. (26)) of the background caustic, and set R * = (σ eff /d) 1/2 due to Katz et al. (1986). Panels correspond to different values for the characteristic source size σ W . Good agreement is found between direct numerical evaluation of Eq. (4) (black solid curves), and our analytic model (Eq. (24); red dashed curves) which approximates the PDF for the microlensing random deflections, P 1 (|α ml |), as an isotropic two-dimensional Gaussian distribution. * = 1/σ eff (see Eq. (26)). For small source sizes, σ W < θ κ 1/2 and hence 1/ * = σ eff ≈ θ κ 1/2 , rendering this second logarithm numerically negligible. On the other hand, this second logarithm increases the effective Gaussian width for large sources, i.e. σ eff ≈ σ W > θ κ 1/2 .
From the above analysis, we obtain the following approximation formula for the mean magnification factor, where the value of σ 2 ml (R * , * ) must be appropriately set using Eq. (25).

Two-point deflection statistics
Apart from the mean magnification µ W (y) , we would like to know the statistical fluctuations in the magnification factor between random microlens realizations. As we have shown in Eq. (6), calculating the second moment µ W (y) 2 requires the knowledge of the two-point CF for microlensing deflection: Here ϕ 2 [ 1 , 2 ; x 1 , x 2 ] is the two-point CF due to a single microlens, and is given by the following integral part that has quadratic dependence on 1 and 2 (with additional logarithms as a result of nonanalyticity). Keeping only this quadratic part amounts to approximating the two-point joint PDF for the random deflections as a two-dimensional normal distribution, which will lead to analytic simplification of the expression for µ W (y) 2 . We argue that this approximation is reasonable if the effective source size overlaps multiple micro-caustics θ /(κ 1/2 µ B σ eff ) 1, which, as we have argued before with Eq. (27), necessarily implies that θ 2 * /R * 1. In this regime, 1 and 2 are typically on the order of * , and the separation between the two image-plane positions is typically comparable to the macro image extent r 12 ≡ |x 1 − x 2 |. With the condition θ 2 1,2 /|x 1 − x 2 | < 1, an analytic approximation for the integral Eq. (31) can be obtained. This is given by Eq. (A11). See Appendix A for the derivation. The result is the following approximation for the two-point CF: Terms of | 1 | 2 and | 2 | 2 encode auto correlation at a single image-plane point and are consistent with the one-point CF in Eq. (23), while the cross terms encode cross correlation between a pair of imageplane points. Since the cross terms only depend on r 12 = x 2 − x 1 , but not on x 1 and x 2 separately, Eq. (33) preserves the expected symmetries of the two-point microlensing deflection statistics under spatial translations and rotations on the image plane. Through a direct calculation of α ml,i (x 1 ) α ml,j (x 2 ) , we verify that the second line of Eq. (33) can be rewritten as − C ml (r 12 ; R * ) ( 1 · r 12 ) ( 2 · r 12 ) r 2 12 + C ml ⊥ (r 12 ; R * ) 1 · 2 − ( 1 · r 12 ) ( 2 · r 12 ) r 2 12 .
where C ml (r; R * ) and C ml ⊥ (r; R * ) are respectively the two-point correlation functions for random microlensing deflections parallel and perpendicular to the separation vector (as defined via decomposition Eq. (8)): Both functions are sensitive to the "infrared" cutoff scale R * and are logarithmically divergent in the limit r → 0, so the correlation at zero separation is ill-defined. This is unlike our toy model of Gaussian random deflections with a regular lensing potential power spectrum, as we have studied in Section 3, for which C (0) = C ⊥ (0) is finite. To our knowledge, Eq. (35) and Eq. (36) have not been presented before in the literature.
Based on Eq. (33), we derive an approximate formula for µ W (y) 2 in a form identical to Eq. (12): . (37) where we introduce the covariance matrix appropriate for random point microlenses: Here on the diagonal we use which depends on the choice for R * and * . Note that C ml (0) is not given by the r → 0 limit of C ml (r) or C ml ⊥ (r). The covariance matrix C(r) being positive definite requires C ml (0) ± C ml (r) > 0 and C ml (0) ± C ml ⊥ (r) > 0. Unlike in the Gaussian random deflection model, this is not strictly guaranteed here, which is a shortcoming of our analytic approximation. However, violation occurs in two regimes, r R * or r θ 2 * , both of which are not expected to have important contributions to µ W (y) 2 . In practice, we must regularize the logarithmic divergences in order for the integral Eq. (6) to be well behaved. For an example, we present in Appendix B one regularization scheme, which, as we numerically test out, renders the result for µ W (y) 2 insensitive to regularization parameter choices.
Our results are readily applicable to the special case of uniform background convergence and shear. This situation has already been intensively studied, mostly in the context of quasar microlensing. In Appendix C, we derive additional analytic results for this special case, and remark on comparisons to the literature.
While Eq. (37) assumes a Gaussian source, Eq. (33), Eq. (35) and Eq. (36) are generally valid independent of the specific source profile. For stellar photospheres, the uniform disk would be a more appropriate source model than the Gaussian one. However, the 1 -and 2 -integrals cannot be analytically carried out, unlike for a Gaussian source, which is a shortcoming of our results.

Numerical experiments
We validate the semi-analytic results we have derived for µ W (y) and µ W (y) 2 using numerical experiments. We set parameters θ = 1 and d = 10 −5 as adopted in Figure 3. Random microlenses with identical θ are generated within a circular disk that centers on the macro critical curve and has a radius R * = 1500. We efficiently compute the summed deflection from a large number of microlenses using the hierarchical tree algorithm (Wambsganss 1999). For κ = 0.004, 0.02, 0.1 that we simulate, we need to include for each realization N ≈ 9000, 45000, 225000 point lenses, respectively. We sample the image-plane vicinity of the macro critical curve with a large number of rays. These rays are inversely traced onto the source plane (Kayser et al. 1986;Wambsganss 1992), which can then be used to calculate µ W (y) for any source profile and central position.
We numerically derive µ W (y) and µ W (y) 2 by averaging over many random realizations for the microlenses. To be consistent with the averaging procedure we adopt in the numerical experiment, we always set R * = 1500 when evaluating C ml (0), C ml (r) and C ml ⊥ (r); this is different from setting R * to be the clustering size of the micro images as proposed in Katz et al. (1986). As we show in Figure 6, for a range of parameters our semi-analytic calculations agree with numerical results to high accuracy, provided that the source typically overlaps with multiple micro caustics. The magnification µ W (y) in fact can have a rather skewed, non-Gaussian distribution when the relative fluctuation is large, while the mean and variance are still accurately predicted by our semi-analytic formulae.
An interesting observation can be made from comparing the last three panels of Figure 6: µ W (y) and µ W (y) 2 in fact become insensitive to the source size σ W if θ κ 1/2 σ W , even though that the light curves are qualitatively distinct. The numerical results show that as σ W decreases, the light curve becomes increasingly non-Gaussian while preserves µ W (y) and µ W (y) 2 ! It is therefore reasonable to hypothesize that as long as the source size is much smaller than θ κ 1/2 , the results for µ W (y) and µ W (y) 2 are insensitive to the source profile either; if the Gaussian source is replaced with a uniform-disk source, our results for µ W (y) and µ W (y) 2 should remain correct. The cases we examine in Figure 6 all correspond to a sufficiently large source size σ W that overlaps multiple or at least order unity micro caustics, regardless of the size of the microlensing broadening θ κ 1/2 . As we can see from the "light curves", the fluctuations of the magnification factor are Gaussian or weakly non-Gaussian. To further test the range of validity of our approximation, in Figure 7 we examine cases where the number density of micro caustics are reduced and the source size σ W is made smaller. The fluctuations of the magnification factor become highly non-Gaussian and very dramatic, approaching the familiar behavior of small sources exhibiting intermittent "flares" at micro caustic crossings. In these cases, the physical source extent hardly overlaps multiple micro caustics, while our semi-analytic approximation remains successful. We note that in these cases multiple micro images still arise (albeit the number of micro images is small), which may explain the success of the approximation. Hence, we find robust numerical evidences that the semi-analytic approximation developed in this work is applicable to computing the mean and variance of the magnification factor over a wide range of parameters.

microlens mass distribution
So far, results have been derived assuming identical microlens masses. Generalization to an arbitrary distribution of Einstein radii is straightforward if microlenses of different masses thoroughly mix in space. Introduce a differential contribution of the microlens convergence dκ /d ln θ 2 . Formally, Eq. (25) must be modified to where θ 2 is the squared Einstein radius for the mean microlens mass. If we still set * = 1/σ eff , then in the defining equation for σ eff , Eq. (26), κ θ 2 also needs to be modified similarly to account for a distribution of Einstein radii. Following a similar logic, in using Eq. (33) we must replace the first line with −(1/2) (| 1 | 2 + | 2 | 2 ) multiplying Eq. (40), and replace κ θ 2 in the second line with the appropriate averaged quantity d ln θ 2 (dκ /d ln θ 2 ) θ 2 .  Figure 6. Statistics of the magnification factor µ W derived from numerical ray-shooting as a function of the distance y 1 to the macro fold caustic. We set parameter values θ = 1, d = 10 −5 and R * = 1500. Results are shown in separate panels for several choices of the microlens surface abundance κ and the source size σ W , all in the regime that the source overlaps multiple micro caustics. Theoretical calculations for the mean magnification µ W (solid blue) and its standard deviation Std[µ W ] (dash-dotted blue), all predicted by the semi-analytic model of this work, agree well with the numerical statistics (solid red curve for µ W and light red band for Std[µ W ]) derived from 200 independent microlensing realizations. Additionally, 10 random realizations are shown as the grey curves. In all panels, the macro caustic is located at y 1 = 0 (vertical dotted line), and the magnification factor for a point source is shown as the dashed black curve.
If microlenses do not differ in mass by orders of magnitude, the second integral is expected to be suppressed by the logarithmic factor, while the first integral is proportional to the average squared mass θ 4 (i.e. weighted toward the more massive microlenses). However, the second integral may not be small at all when there is a hierarchy in θ 2 . Interestingly, for galactic or intracluster stars a large mass hierarchy does exist between the sub-solar main-sequence (MS) dwarfs and the remnant black holes (BHs). For an old stellar population of which all stars with initial masses > 1 M have become stellar remnants, and assuming an IMF dφ ( Figure 6, but for a relatively low microlens surface density κ = 0.0004 and reduced source sizes. Such small sources do not overlap multiple micro caustics, except in the very vicinity of the macro caustic where the macro magnification factor is sufficiently high µ B 2000. In spite of the highly non-Gaussian nature of the magnification fluctuations, our semi-analytic predictions for the mean magnification and its variance still agree well with numerical simulations.
8 M for the BHs, the BHs can make a comparable contribution to θ 4 , if not more, than the MSs. This implies the importance of BH microlenses in broadening the "point spread function" of random deflections despite their low number fraction. A detailed investigation into a realistic mass function will be included in a future work.

DISCUSSION
As we have explained, stochastic microlensing has a profound effect on the magnification factor of a lensed source. Now we discuss this more in the context of strong lensing produced by galaxy and galaxy cluster lenses.
The characteristic scale of random microlensing deflection for a microlens mass M corresponds to a source-plane scale where D L , D S and D LS are the angular diameter distances to the lens plane, to the source plane, and from the lens plane to the source plane, respectively. We have multiplied by a factor of 3 to account for the Coulomb logarithm (see Eq. (25)). In galactic lenses, the surface density of stellar microlenses is high κ 0.1-1; hence R is smaller than the typical size of a star cluster, comparable to or larger than the sizes of optical quasars ∼ 10 3 AU (Blackburne et al. 2011), while certainly larger than individual stellar photospheres. In galaxy cluster lenses, the intracluster stars have a substantially lower surface density κ ∼ 0.001-0.01, and R can be reduced by up to a factor of ten.
As long as the source's physical size is smaller than R , the highest persistent magnification factor µ W (y) reached at a macro caustic is on the order of For galaxy lenses, the choices d −1 ∼ 1 and κ ∼ 0.1 are reasonable, and hence only sufficiently compact sources can possibly acquire a temporary magnification much higher than 1000, at micro caustic crossings. This still requires that µ W can fluctuate to a value much higher than µ W (y) . A conservative constraint is that the source size σ W is smaller than ∼ θ κ −1/2 / µ W max , the typical separation of micro caustics on the source plane. This limits the source size to for which only individual source stars meet the requirement. However, the fluctuation in µ W can still be significantly suppressed for small sources if many micro images form (σ eff > θ κ −1/2 / µ W max ), because only one pair of micro images are enhanced at each micro caustic crossing.
The best opportunities to have very high temporary magnifications for individual stars are to be found in galaxy cluster lensing with a small κ ∼ 0.001-0.01, for which the maximum values (at the tail of the distribution) can range from a few thousands to 10 4 (Diego 2019). The maximal mean magnification Eq. (42) can now reach µ W max ∼ 3000 (1−κ 0 ) −1 for κ = 0.01 and d −1 = 10 , and the constraint on the source size Eq. (43) is relaxed, which is also helped by the fact that the typical value of d is reduced in cluster lenses. For galaxy lenses with κ 0.1, we do not expect the magnification to strongly fluctuate and reach significantly higher than 1000, as the corrugated micro caustic network is too dense. This conclusion is also reached from the argument that the peak magnification at a micro caustic crossing scales as κ −3/4 (Venumadhav et al. 2017). From the perspective of this work, if we choose κ 0 = 0.7, d −1 = 1 , θ = 1 µas and κ = 0.3, our semi-analytic approximation for µ W and µ 2 W is applicable because σ eff > θ κ −1/2 / µ W max , and a small standard deviation ∼ 10-20% for the fractional magnification fluctuation around µ W ∼ 1000-2000 is predicted, insensitive to the source size. If we instead set κ = 0.005, for the same macro caustic, we find a larger standard deviation 30% for the fractional magnification fluctuation, around a much higher mean µ W ∼ 2000-6000. This example is shown in Figure 8.
It is worth to note that sub-galactic DM subhalos as substructure lenses tend to strongly perturb a galactic or cluster caustic and create secondary caustics under suitable conditions (Dai et al. 2018. An interesting consequence of these subhalos is then to increase the typical value of d (i.e. weaken the caustic strength) and hence further reduce the allowed maximum mean magnification, in both galaxy and cluster lenses.
Applying the same analysis to a variety of large sources tens of AUs, which include quasars, SNe (Kelly et al. 2015;Goobar et al. 2017), or bloated stellar photospheres due to outburst or mass ejection, µ W is expected to fluctuate only mildly around the mean value in Eq. (42), so magnifications significantly higher than ∼ 1000 are prohibited by microlensing, even for galaxy cluster lenses. While multiply-imaged quasars commonly have magnification factors on the order O(10), quasars magnified by a hundred to a thousand fold are rarely reported. Fujimoto et al. (2020) suggested a candidate lensed quasar with a total magnification ∼ 450. However, analysis of the proximity zone does not seem to support this idea (Davies et al. 2020). In another quadruply-imaged quasar, a ten-fold magnification anomaly was detected for one of the images, requiring a magnification factor as large as ∼ 100 (Glikman et al. 2018). These large magnifications are likely to be consistent with the maximal values permitted by microlensing effects, if the caustic strength is not dramatically reduced by subhalos, i.e. d −1 0.1 . We expect that microlensing effects dominate the truncation in the = 0.005 = 1 as d 1 = 1 Figure 8. The mean and variance of the magnification factor µ W in the vicinity of a macro caustics with κ 0 = 0.7 and d −1 = 1 computed using the semi-analytic approximation developed in this work. We contrast between a high microlens surface density κ = 0.3 (left), typical of galaxy lensing, and a low microlens surface density κ = 0.005 (right), typical of cluster lensing. Both the mean magnification and the fluctuation around it are significantly suppressed in the former case by the excessively high number of micro caustics. While µ W and µ 2 W are insensitive to the source size σ W as long as σ W θ κ 1/2 (true for stellar photospheres), the magnification distribution becomes increasingly non-Gaussian for smaller sources.
high magnification tail of the lensed quasars (especially for low-mass quasars and in the case of small Einstein radii), which may have implications for the impact of magnification bias on the luminosity function (Pacucci & Loeb 2019.

CONCLUSION
Gravitationally lensed sources exhibit stochastic fluxes as a result of random microlensing if compact masses contribute a fraction of the lens surface mass. Through a first-principle statistical treatment of microlensing deflections, we have in this work derived a semi-analytic approximation for the mean and variance of the magnification factor, for a finite Gaussian source and for arbitrary macro lens models. A theoretically appealing feature of the new result is that the UV and IR logarithms are physically determined.
These general results are in the form of single and double image-plane integrals with simple and well-behaved integrands, and hence are practically useful as these can be efficiently evaluated using Monte Carlo integrators. Our analytic derivations suggest that the results are good approximations if the source of an effective size σ eff overlaps multiple micro caustics, where the effective size is either the source's physical size σ W or the scale of random microlensing deflections θ κ 1/2 , whichever is larger. Using numerical ray-shooting with random microlens realizations, we have demonstrated the accuracy of the approximation, even in cases where the microlensing-induced light curves are highly non-Gaussian. While we have specifically examined highly magnified sources near a macro fold caustic, for which a small convergence from the microlenses can induce dramatic flux variance, our results are readily applicable to other macro caustics, such as a cusp caustic, or higher-order catastrophes (Feldbrugge et al. 2019). We have pointed out that the maximal magnification that can be realized at a macro fold caustic is not only limited by the source size σ W , but also by the characteristic scale of microlensing deflections θ κ 1/2 , especially for a source that overlaps multiple micro caustics or has many micro images. Future work may adopt the formalism here to study the correlation of the magnification factor between two different source-plane positions, i.e. µ w (y 1 ) µ w (y 2 ) . For a moving source, this translates to the temporal correlation of microlensing lightcurves (Wyithe & Turner 2002;Lewis & Irwin 1996;Neindorf 2003), and will be useful for interpreting cadence observations. Another interesting question regards the third-order moments of the magnification factor, as well as higher-order moments, which characterize the departure from Gaussian statistics. Our formalism may be applicable to the computation of these moments, which can help with the analyses of highly non-Gaussian light curves.
We only concern terms that have quadratic dependence on the Fourier wave vectors (the second line). First, let us examine the contribution: We are supposed to perform the integration over |z| < R excluding |z − x 1 | < r and |z − x 2 | < r. However, in the regime R r and |r 12 | r, we introduce an error suppressed by (r/R) 2 if instead we evaluate − θ 4 2 1 π R 2 r<|z|<R d 2 z ( 1 · z) 2 |z| 4 = − θ 4 2 1 2 R 2 ln R r . (A8) Here, we essentially neglect the |z − x 2 | < r disk, and move the |z − x 1 | < r disk to the center of the |z| < R disk. Similarly, we have another contribution −θ 4 2 2 /(2 R 2 ) ln(R/r) from the 1 ↔ 2 and x 1 ↔ x 2 crossed term.
We seek a regularization scheme such that: (1) for r θ 2 * , C ml (r) and C ml ⊥ (r) are always smaller than but approach C ml (0); (2) for r R * , C ml (r) and C ml ⊥ (r) are asymptotically zero. Specifically, we define the following regularized functions: C ml X,reg (r) := C ml (0) 1 + C ml (0) |C ml X (r)| e −(r/(ν R * )) 2 n −1/n sgn[C ml X (r)], where X is either or ⊥. We find that a suitable choice for the two parameters are n = 10 and ν = 1. An numerical example is plotted in Figure 10. As far as the computation of µ W (y) 2 is concerned, what is important is the regime θ 2 * < r < R * , where both C ml (r) and C ml ⊥ (r) are linear functions of ln r.

C. UNIFORM BACKGROUND CONVERGENCE AND SHEAR
In this Appendix, we apply the analytic approximation derived in Section 4 to a special case: the background lens has a constant convergence κ B and shear (γ B,1 , γ B,2 ). We can write where we have defined the background deformation matrix A B . Thus, the background lens model has translational invariance. Without loss of generality, we may set γ B,1 = γ B and γ B,2 = 0. The moments of µ W (y) are independent of y, so we are free to set y = 0.