Neutrino and Photon Lensing by Black Holes: Radiative Lens Equations and Post-Newtonian Contributions

We extend a previous phenomenological analysis of photon lensing in an external gravitational background to the case of a massless neutrino, and propose a method to incorporate radiative effects in the classical lens equations of neutrinos and photons. The study is performed for a Schwarzschild metric, generated by a point-like source, and expanded in the Newtonian potential at first order. We use a semiclassical approach, where the perturbative corrections to neutrino scattering, evaluated at one-loop in the Standard Model, are compared with the Einstein formula for the deflection using an impact parameter formulation. For this purpose, we use the renormalized expression of the graviton/fermion/fermion vertex presented in previous studies. We show the agreement between the classical and the semiclassical formulations, for values of the impact parameter $b_h$ of the neutrinos of the order of $b_h\sim 20$, measured in units of the Schwarzschild radius. The analysis is then extended with the inclusion of the post Newtonian corrections in the external gravity field, showing that this extension finds application in the case of the scattering of a neutrino/photon off a primordial black hole. The energy dependence of the deflection, generated by the quantum corrections, is then combined with the standard formulation of the classical lens equations. We illustrate our approach by detailed numerical studies, using as a reference both the thin lens and the nonlinear Virbhadra-Ellis lens.


Introduction
According to classical general relativity (GR) massless particles follow null spacetime geodesics which bend significantly in the presence of very massive sources. The gravitational lensing enforced on their spatial trajectories provides important information on the underlying distributions of matter and, possibly, of dark matter, which act as sources of the gravitational field. Several newly planned weak lensing experiments such as the Dark Energy Survey (DES) [1], the Large Synoptic Survey Telescope (LSST)[2], both ground based, or from space with the Wide-Field Infrared Survey Telescope (WFIRST) [3] and Euclid [4], are expected to push forward, in the near future, the boundaries of our knowledge in cosmology. In the analysis of the deflection by a single compact and spherically symmetric source, one significant variable, beside the mass of the source, is the impact parameter of the incoming particle beam, measured respect to the center of the source, which determines the size of the deflection. It is very convenient to measure the impact parameter (b), which is typical of a given collision, in units of the Schwarzschild radius r s ≡ 2GM , denoted as b h ≡ b/r s . In the Newtonian approximation for the external background, this allows to scale out the entire mass dependence of the lensing event.
For an impact parameter of the beam of the order of 10 5 − 10 6 , the corresponding deflection is rather weak, of the order of 1-2 arcseconds, as in the case of a photon skimming the sun. Stronger lensing effects are predicted as the particle beam nears a black hole, with deflections which may reach 30 arcseconds or more. These are obtained for impact parameters b h of the order of 2×10 4 . Even larger deflections, of 1 to 2 degrees or a significant fraction of them, are generated in scatterings which proceed closer to the event horizon [5]. In fact, as we are going to show, for closer encounters, with the beam located between 20 and 100 b h , such angular deflections are around 10 −2 radians in size, as predicted by classical GR. A high energy cosmic ray of 10-100 GeV will then interact with the field of the source by exchanging momenta far above the MeV region, and will necessarily be sensitive to radiative effects, such as those due to the electroweak corrections. Interactions with such momentum exchanges cannot be handled by an effective Newtonian potential, as derived, for instance, from the (loop corrected) scattering amplitude. We recall that, in general, in the derivation of such a potential, one has to take into account only non-analytic terms in the momentum transfer q. These are obtained from a given amplitude and/or gravitational form factor of the incoming particle after an expansion at small momentum. The analytic terms in the expansion correspond to contact interactions which are omitted from the final form of the potential, being them proportional to Dirac delta functions. As one can easily check by a direct analysis, non-analytic contributions originate from massless exchanges in the loops, which approximate the full momentum dependence of the radiative corrections only for momentum transfers far below the MeV region. Therefore, the validity of the method requires that the typical impact parameter of the beam, for particle of the energy of few GeV, be of the order of 10 6 Schwarzschild radii and not less. For such a reason, if we intend to study a lensing event characterized by a close encounter between a cosmic ray and a black hole, we need to resort to an alternative approach, which does not suffer from these limitations. Finally, with the photon sphere located at b h ∼ 2.5 for a Schwarzschild metric, one expects that very strong deflections are experienced by a beam for scattering events running close to such a value of the impact pa-rameter. This is also the radial distance from the black hole center at which the scattering angle diverges. A Standard Model and V the electromagnetic current. In the fermion case (f), the corresponding vertex is the T f f , with f denoting a neutrino. The comparison between the classical and the semiclassical formula for the deflection derived by this method can then be performed at numerical level, as shown in [5] for the photons. The energy dependence of the bending angle, for a given impact parameter of the photon beam, though small, is found to become more pronounced at higher energies, due to the logarithmic growth of the electroweak corrections with the energy.
The goal of our present work is to propose a procedure which allows to include these effects in the ordinary lens equations, illustrating in some detail how this approach can be implemented in a complete numerical study. We mention that our semiclassical analysis is quite general, and applies both to macroscopic and to microscopic black holes. In the case of macroscopic black holes the procedure has to stop at Newtonian level in the external field. In fact, post-Newtonian corrections, though calculable, render the perturbative expansion in the external (classical) gravitational potential divergent, due to the macroscopic value of the Schwarzschild radius. On the other hand, in the case of primordial black holes, the very same corrections corrections play a significant role in the deflection of a cosmic ray, and bring to a substantial modification of the classical formulas.

Organization of this work
In the first part of our work we will extend a previous analysis of photon lensing [5], developed along similar lines, to the neutrino case, presenting a numerical study of the complete one-loop corrections derived from the electroweak theory. The formalism uses a retarded graviton propagator with the effects of back reaction of the scattered beam on the source not included, as in a typical scattering problem by a static external potential. In this case, however, because of the presence of a horizon, we search for a lower bound on the size of the impact parameter of the collision where the classical GR prediction and the quantum one overlap. Indeed, above the bound the two descriptions are in complete agreement. As already mentioned above, both in the fermion as in the photon case [5], this bound can be reasonably taken to lay around 20 b h , which is quite close to the horizon of the classical source. For smaller values of b h (4 < b h < 20), the two approaches are in disagreement, since the logarithmic singularity in the angle of deflection, once the beam gets close to the photon sphere, starts playing a significant role. This is expected, given the assumption of weak field for the gravitational coupling, which corresponds to the Newtonian approximation in the metric. The second part of our work deals with the implementation of the semiclassical deflection within the formalism of the classical lens equations. We use the energy dependence of the angular deflection to derive new lens equations, which are investigated numerically. We quantify the impact of these effects both in the thin lens approximation, where the trigonometric relations in the lens geometry are expanded to first order, and for a nonlinear lens. As an example, in this second case, we have chosen the Virbhadra-Ellis [12] lens equation. The observables that we discuss are limited to solutions of these equations and to their magnifications, although time delays, shears and the light curves of a typical microlensing event can be easily included in this framework. We anticipate that the effects that we quantify are small and cover the milliarcsecond region, remaining quite challenging to detect at experimental level. We hope though, that the framework that we propose can draw further interest on this topic in the future, both at theoretical and at phenomenological level.
In the third part of our study we discuss the post Newtonian formulation of the impact parameter formalism, and apply it to the case of a compact source with a microscopic Schwarzschild radius. This is the only case in in which the gravitational corrections to the Newtonian cross section can be consistently included in our approach in a meaningful way. We then summarize our analysis and discuss in the conclusions some possible future directions of possible extensions of our work.

Gravitational interaction of neutrinos
We start our analysis with a brief discussion of the structure of the gravitational interaction of neutrinos, building on the results of [13,14], to which we refer for additional details, and that we are going to specialize to the case of a massless neutrino. An analysis of gravity with the fermion sector is contained in [15]. We simply recall that the dynamics of the Standard Model in external gravity is described by the Lagrangian This includes the Einstein term S G , the S SM action and a term S I involving the Higgs doublet H [16], called the term of improvement. S SM , instead, is obtained by extending the ordinary Lagrangian of the Standard Model to a curved metric background. The term χ is a parameter which, at this stage, is arbitrary and that at a special value (χ ≡ χ c = 1/6) guarantees the renormalizability of the model at leading order in the expansion in κ.
Deviations from the flat metric η µν = (+, −, −, −) will be parametrized in terms of the gravitational coupling κ, with κ 2 = 16πG N and with G N being the gravitational Newton's constant. At this order the metric is given as g µν = η µν + κh µν , with h µν describing its fluctuations. We will consider two spherically symmetric and static cases, corresponding to the Schwarzschild and Reissner-Nordstrom metrics. The first, in the weak field limit and in the isotropic form is given by In this case the fluctuation tensor takes the form The inclusion of higher order terms in the weak field expansion will be discussed in the following sections. The coupling of the gravitational fluctuations to the fields of the Standard Model involves the EMT, which is defined as with a tree-level coupling summarized by the action where T µν is symmetric and covariantly conserved. The complete expression of the EMT of the Standard Model, including ghost and gauge-fixing contributions can be found in [17].
The Higgs field is parameterized in the form in terms of h, φ and φ ± , which denote the physical Higgs and the Goldstone bosons of the Z and W s respectively. v is the Higgs vacuum expectation value. The terms of the Lagrangian S I , generate an extra contribution to the EMT which is given by the term of improvement, which can be multiplied by an arbitrary constant (χ). As mentioned above, it is mandatory to choose the value χ = 1/6 for any insertion of the EMT on the correlators of the Standard Model.
These are found to be ultraviolet finite only if T µν I is included [16,17,18]. We will be dealing with the T ff vertex, where T denotes the EMT and f ≡ ν f a neutrino of flavour f , and work in the limit of zero mass of the neutrinos. The vertex, to lowest order, is obtained from the EMT of the neutrino. For instance, the explicit expression of the EMT for the electronic neutrino is given by with In momentum space, in the case of a Dirac fermion, the vertex takes the form We refer to appendix G for a list of the relevant Feynman rules necessary for the computation. We will denote withT the corresponding invariant amplitude, a notation that we will use also at one-loop level in the electroweak expansion. We introduce the two linear combinations of momenta p = p 1 + p 2 and q = p 1 − p 2 to express our results. It has been shown that the general T ff vertex, for any fermion f of the Standard Model, decomposes into six different contributions [13], but in the case of a massless neutrino only three amplitudes at one-loop level are left, denoted asT Figure 1: The one-loop Feynman diagrams of the neutrino vertex in a gravitational background. The dashed lines can be Z and W .
In the expression above, the subscripts indicate the contributions mediated by virtual Z and W gauge bosons, while CT indicates the contribution from the counterterm. We show in Fig. 1 some of the typical topologies appearing in their perturbative expansion. Two of them are characterized by a typical triangle topology, while the others denote terms where the insertion of the EMT and of the fermion field occur on the same point. The computation of these diagrams is rather involved and has been performed in dimensional regularization using the on-shell renormalization scheme. Neutrinos interactions, in the limit of massless neutrinos, involve only few of the structures of the T ff tensor decomposition presented in [13]. In this case we are left with only one tensor structure and hence only one form factor for each sectorT where we have defined the vertex The counterterms needed for the renormalization of the vertex can be obtained by promoting the counterterm Lagrangian of the Standard Model from a flat spacetime to the curved background, and then extracting the corresponding Feynman rules, as for the bare one. We obtain where we have denoted with Σ L the neutrino self-energy which is a combination of the self-energy contributions with expressed in terms of the scalar form factor B 0 , given in appendix E together with all the other relevant scalar integrals. We have denoted with m Z and m W the masses of the Z and W gauge bosons; with q 2 the virtuality of the incoming momentum of the EMT and m f is the mass of the fermion of flavor f running in the loops. The explicit expressions of the form factors appearing in (13) is given by with C 0 denoting the scalar 3-point function, and with the form factor f W 1 related to the exchange of the W ' s given by Being the computations rather involved, the correctness of the results above has been secured by appropriate Ward identities, whose general structure has been discussed in [17]. As an example, by requiring the invariance of the generating functional of the theory under a diffeomorphic change of the spacetime metric, one derives the following Ward identity where Γf f (p) is the fermion two-point function, diagonal in flavor space [13]. From this equation one obtains which, as one can check, are identically satisfied by the explicit expressions of f Z and f W given above.
In the case of MeV neutrinos, the expressions of the two form factors simplify considerably, since the typical momentum transfers q 2 = −4E 2 sin 2 (θ/2) may be small. These expansions, in fact, are useful in the case of scattering and lensings of neutrinos far from the region of the event horizon, of the order of 10 3 − 10 6 horizon units. As we are going to see, an expansion in q 2 provides approximate analytical expressions of the b h (α) relation, connecting the impact parameter to the angle of deflection α, valid at momentum transfers which are smaller compared to the electroweak scale, i.e. q 2 /m 2 W 1. We will come back to illustrate this point more closely in the following sections. In these cases the expression of the renormalized f Z form factor takes the form while the W form factor is slightly lengthier 3 Neutrino scattering at leading order in the external field Using the results of the previous section, we proceed by giving the expression of the neutrino gravitational cross section in the electroweak theory. The approach followed has been described in previous works [13,14]. We just recall that the scattering matrix element is written as where V is the integration volume where the scattering occurs, which gives Denoting with i and f the initial and final neutrino, we have introduced plane waves normalized as and similarly for ψ f , while V denotes a finite volume. The E 1 (E 2 ) are the energy of the incoming (outgoing) particle respectively. In momentum space the matrix element is given by (29) in terms of the gravitational fluctuations in momentum space h µν (q). For a static external field the energies of the incoming/outgoing fermions are conserved (E 1 = E 2 ≡ E).
The Fourier transform of h µν in momentum space is given by which for a static field can be expressed as in terms of a single form factor h 0 ( q) The leading order cross section of a neutrino scattering off a Schwarzschild background, in the weak field limit, is given by which is energy independent. Its expression at 1-loop level, instead, takes the form whose explicit expression has been given in appendix C. In the massless approximation for the neutrino masses, loop corrections do not induce flavor transition vertices, such as those computed in [14].
In the case of neutrinos of an energy E in the MeV range, the expression above simplifies considerably and takes the form We show in Fig. 3 three plots of the tree level and one-loop cross sections for an energy of the incoming neutrino beam of 1 MeV, for 2 different angular regions (plots (a) and (b)), together with a global plot of the entire cross section (plot (c)) for the rescaled differential cross section dσ/dΩ ≡ 1/r 2 s dσ/dΩ. Notice that the tree-level and one-loop results are superimposed. We can resolve the differences between the two by zooming-in in some specific angular regions of the two results, varying the energy of the incoming beam. The result of this analysis is shown in Fig. 4, where in plots (a) and (b) we show the rescaled cross section dσ/dΩ as a function of the scattering angle θ, for three values of the incoming neutrino beam equal to 1 GeV, 1 TeV and 1 PeV. PeV neutrinos events are rare, due to the almost structureless cosmic ray spectrum, which falls dramatically with energy. They could be produced, though, as secondaries from the decays of primary protons of energy around the GZK [19,20] cutoff, and as such they are part of our analysis, which we try to keep as general as possible. It is clear from these two plots that the tree level and the one-loop result are superimposed at low energies, with a difference which becomes slightly more remarked at higher energies. A similar behaviour is noticed in the cross section for scatterings at larger angles. Also in this case the radiative corrections tend to raise as the energy of the incoming beam increases. This behaviour is expected to affect the size of the angle of deflection α as we approach the singular region of a black hole. In fact, α is obtained by integrating the semiclassical equation (36), introduced below, and large deviations are expected as the impact parameter b h reaches the photon sphere. As we are going to illustrate in the next sections, the b h (E, α) relation is significantly affected by the behaviour of the cross section at large θ as b h → 3/2r s . This is the closest radial distance allowed to a particle approaching the black hole from infinite distance without being trapped. Therefore, these differences inσ for large θ are going to render b h sensitive on the changes in energy of the neutrino beam for such close encounters of the neutrinos with a black hole.
4 Impact parameter formulation of the semiclassical scattering As pointed out in previous studies [5,7,14,21], the computation of the angle of deflection for a fermion or a photon involves a simple semiclassical analysis, in which one introduces the impact parameter representation of the specific classical cross section and equates it to the quantum one. The classical/semiclassical scattering process is illustrated in Fig. 2, with α denoting the angle of deflection. By assuming that the incoming particle is moving along the z direction, with the source localized at the origin, and denoting with θ the azimuthal scattering angle present in the quantum cross section, we have the relation between the impact parameter b and θ, as measured from the z-direction. This semiclassical equation [7,21] allows to relate the quantum and the classical features of the interaction between the particle beam and the gravitational source. The explicit expression of b(α), at least for small deflection angles, which correspond to large values of the impact parameter, can be found either analytically, such as at Born level and, for small momentum transfers also at one-loop, but it has to be obtained numerically otherwise. The solution of (36) takes the general form with b 2 h (θ) denoting the constant of integration. The semiclassical scattering angle α is obtained from (37) as a boundary value of the integral in θ of the quantum cross section. As discussed in [5], the integration constant derived from (37) has to be set to zero (forθ = π) in order for the solution of (36) to match the classical GR result for a very large b h .
In the case of a point-like gravitational source and of neutrino defection, one obtains from (37) the differential equation Notice that the variation of b with the scattering angle θ is negative, since the impact parameter decreases as θ grows, as we approach the center of the massive source. A comparison of this expression with the analogous relation in the photon case (γ) shows that the two equations differ by a simple prefactor The solution of (38) takes the form and in the small α (i.e. large b) limit takes the asymptotic form which allows us to identify the deflection angle as in agreement with Einstein's prediction for the angular deflection. This is the result expected from the classical (GR) analysis. The inversion of the asymptotic expansion (41) generates the asymptotic behaviour which corresponds to the general functional form As shown in Fig. 5, the analytic inversion of (41), given by (43), is very stable under an increase of the order of the asymptotic expansion over a pretty large interval of b h , from low to very high values. Solutions (43) and (44) can be obtained by an iterative (fixed point) procedure, which generates a sequence of approximations implemented after a Laurent expansion of (41) and the use of the initial condition α 0 = 2/b h . The approach can be implemented also at one-loop and with the inclusion of the post-Newtonian corrections, if necessary. The logarithmic corrections present in (44) are a genuine result of the quantum approach and, as we a re going to discuss below, are not present in the classical formula for the deflection. Radiative and post-Newtonian effects, not included in (43), give an expression for α(b h ) which coincides with the form (44), with specific coefficients (a n , d n ) which are energy dependent. This is at the origin of the phenomenon of light dispersion (gravitational rainbow) induced by the quantum corrections, which is absent at classical level. Eq. (44) will play a key role in our proposal for the inclusion of the radiative corrections in the classical lens equation. Such equation will relate the angular position of the source in the absence of lensing, β, to α(b).
Moving to the one-loop expression given in (34), we can derive an analytic solution of the corresponding semiclassical equation (36) for b = b(E, α), in the limit of small momentum transfers. For this reason we perform an expansion of (34) in with the coefficients C, D, F and G are functions of the energy and of the masses of the weak gauge bosons. Their explicit expressions can be found in appendix D. The impact parameter b h (α), as shown in the same appendix, has a dependence on the angular deflection α which can be summarized by an expression of the form that we can invert in order to get α(E, b h ). This is given by We show in Fig. 6 some plots of the impact parameter b h as a function of the deflection angle in a range closer to the horizon of a black hole, computed using the Newtonian approximation derived from the metric (2). The region involved covers the interval between 20 and 100 horizons. The numerical results refer to the GR solution and to the full one-loop prediction respectively. The classical expression and the quantum one start differing as we approach the value of b h ∼ 20, and are characterized by a certain dependence on the energy of the incoming beam. Shown are the plots corresponding to neutrinos of energies in the TeV and the PeV range respectively. In these regions the lensing is very strong, corresponding to 10 3 arcseconds and larger. As the neutrino (or the photon) beam gets closer to the photon sphere (x 0 = 3/2r s ), which is the point of maximum approach, the angular deflection diverges. This is the impact parameter region where one expects the formation of relativistic images. The divergence can be parameterized by an integer n, with α n = 2πn, and n tending to infinity. The integer is the winding number of the beam path around the photon sphere. In the external neighborhood of the point of closest approach the beam still escapes to infinity, forming an infinite set of images which are parameterized by the same integer n [6].

1/b n contributions to the deflection
It is interesting to compare the classical GR prediction for the deflection with the result of (44), by resorting to a similar expansion for the deflection integral. This has been studied quite carefully in the literature, especially in the limit of strong lensing [22,23]. The 1/b n h expansion has been shown to appear quite naturally in the post-Newtonian approach applied to the Einstein integral for light deflection. We recall that Einstein's expression in GR is given by the integral and can be re-expressed in the form with the variable s ≡ r s /(2r 0 ) being related to the ratio between the Schwarzschild radius and the distance of closest approach between the particle and the source, r 0 . Th exact computation of this integral is discussed in appendix F, and involves elliptic functions. Additional information on α(r 0 ) is obtained via an expansion of the integrand in powers of s and a subsequent integration. This method shows that the result can be cast in the form with The coefficients a i differ from those given in [23] (up to a 7 ) just by a normalization. They are obtained by re-expressing s = s(r 0 ) in terms of the impact parameter b h using the relation between the impact parameter and the radial distance of closest approach, having redefined x 0 ≡ r 0 /(2 GM ). This can also be brought into the form An expression equivalent to (53) can be found in [5]. Eq. (53) can be given in a 1/b h expansion which will turn useful below. We can invert (50) obtaining the relation which differs from (44) by the absence of logarithmic terms in the impact parameter b h and by the energy independence of the coefficients. The inclusion of the extra contributions mentioned above, in the classical GR expression, becomes relevant in the case of strong lensing. We show in appendix B how the inclusion of the additional 1/b n h terms in the expansion of the angular deflection can be extended to the case of a continuous distribution of sources/deflectors. This provides a simple generalization of the standard approach to classical lensing for such distributions.

Lens equations and 1/b n corrections
The standard approach to gravitational lensing in GR is based on an equation, derived from a geometrical construction, which relates the angular position of the image (θ I ) to that of the source (β), with an intermediate angular deflection (α) generated on the lens plane. In this section we are going to briefly review this construction, which is based on the asymptotic expression for the angular deflection (α ∼ 2/b h ), and discuss its extension when one takes into account more general expansions of α(b h ) of the form given by Eq. (50). The extension that we consider covers the case of a thin lens and concerns only the extra 1/b n h terms derived from classical GR. The discussion is preliminary to the analysis of the next section, where we will consider the inclusion of the radiative effects, parameterized by (44), into the classical lens equation.

The lens geometry
We show in Fig. 7 the lens geometry in the case of a continuous distributions of sources and deflectors. A simplified picture of the geometry, with pointlike source and deflector is shown in Fig. 8. We indicate with β the oriented angle between the optical axis (OP ) (taken as the z axis) and the unlensed direction of the source (OS). θ I denotes the angle formed by the visual line of the image (OI) with the optical axis. We also denote with D OL the distance between the observer and the lens plane; with D LS the distance between the lens plane and the source plane and with D OS the distance of the source plane from the observer.α is the (oriented) angle of deflection, measured clockwise as all the other angles appearing in the geometrical construction. We also introduce the relations, valid for D LS , D OL much larger than the size of the lens, typical of a linear lens, The thin lens equation follows from the approximate geometrical relation Denoting with ξ a 2-D vector in the lens plane, it is convenient to introduce two scales η 0 and ξ 0 defined as Using the lens equation in the geometric relation we find the relation which defines the thin lens equation. It is possible to give a simpler expression to the equation above if we go back to (57) and perform simple manipulations on the angular dependence. On the lens plane ( Fig. 8) the equation takes the scalar form which can be extended to the case of stronger lensing by the inclusion of the contributions of the 1/b n corrections in α(b). Use of the Einstein relation α = 4GM/b and of the relation b ∼ θ I D OL brings (61) into the typical form which defines the thin lens approximation, with θ E being the Einstein radius. For a source S aligned on the optical axis together with the deflector and the observer O (see Fig. 9) -which is defined by the segment  connecting the observer, the lens and the plane of the source (with β = 0) -the images will form radially at an opening θ I = θ E and appear as a circle perpendicular to the lens plane. For a generic β, instead, the primary and secondary image solutions are given by the well-known expressions away. This generates a difference in the solid angles under which the source is viewed by the observer in the unlensed and in the lensed cases. In the simple case of an axi-symmetric lens the ratio between the two solid angles can be defined in the scalar form In the case of a thin lens (60), the analogous expression is given by For this lens the analysis simplifies quite drastically. Using the expression of the two images θ I± given in (63) one obtains the simple expression for the primary and secondary images where the Einstein angle is defined as usual It is convenient to measure the angular variables in terms of the Einstein angle θ E , asβ ≡ β/θ E ,θ ≡ θ I /θ E , withθ then the total magnification takes a rather simple form This equation is commonly used to calculate the light curve in the microlensing case. We refer to [24] for a short review on this point.

Nonlinear effects in strong deflections
In conditions of strong lensing, the linear approximations in the trigonometric expressions are not accurate enough and one has to turn to a nonlinear description of the geometry, expressed in terms of the angular variables which are involved. We illustrate this point by taking as an example a typical lens equation, which in our case is given by the Virbhadra-Ellis construction (VE) [25]. Following Fig. 8, we recall that the VE lens equation is based on the geometrical relation [25,26] P S = P I − SI, which gives  under the assumption that the point R in Fig. 8 lies on the vertical plane of the lens. θ I is the angle at which the image is viewed by the observer and β is the unlensed angular position of the source. Within this approximation we can use the geometric relation b = D OL sin θ I , which allows to relate the image position θ I to the angular deflection of the beam α. Notice that this approximate relation is justified by the fact that the distances D OL and D OS are very large compared to the radius of closest approach r 0 . In this limit the two segments V H and V T are treated as equal.
We remind that (73) is not the unique lens equation that one can write down, but, differently from Eq. (60), it can be used in the case of strong lensing. It takes into account the nonlinear contributions to the angular deflection by the introduction of the tan(β) and tan(θ I ) terms, which in (60) are not included. We refer to [26] for a review of possible lens equations.
7 Radiative effects and the geometry of lensing Turning to our case study, radiative effects in the lens equations can be introduced by replacing the expression of the angular deflection generated by the source on the source plane, which is a function of the impact parameter b (α = α(b)) with the new, energy dependent relation α(b, E) whose general form is given by (44).
For simplicity we consider a pointlike source, and a pointlike deflector, as shown in Fig. 8. We recall that for a massless particle the geodesic motion is determined in terms of the energy E and of the angular momentum L at the starting point of the trajectory. The gravitational deflection, however, can be written only as a function of the impact parameter b of the source, with b = E/L, which is an important result of the classical approach. For a further clarification of this aspect, which differs from the semiclassical analysis we are interested in, we briefly overview the classical case, using the lens geometry as a reference point for our discussion. For a source located on the source plane at an angular opening β (in the absence of the deflector), the initial conditions can be expressed in terms of the two components of the initial momentum p = (p r , p φ ) on the plane of the geodesic, or, equivalently, by the pairs (p r , E) or (p ϑ , E), with E the initial energy of the beam. We recall that for a Schwarzschild metric these are defined as We have denoted withẋ ≡ dx/ds the derivative respect to the affine parameter. p t and p φ related to the energy and to the angular momentum as p t = E and p φ = −L, and with the motion taking place on the plane ϑ = π/2 (p ϑ = 0). They are constrained by the mass-shell condition with (p r =ṙ, p φ =φ, p t =ṫ). The lens equation, usually written as L(β, θ I ) = 0, can also be written, equivalently, in the form of a constraint between β and b using (74). We can use any of the independent variables mentioned above. For a given initial momentum of the beam, emitted from the plane of the source, the lens equation will then determine the position of the source in such a way that the geodesic motion will reach the observer at its location on the optical axis. In particular, an interesting description emerges if we choose as initial conditions the angular position of the source (β) and the value of the impact parameter b. These two conditions fix the direction of the trajectory of the beam at its origin on the source plane. In these last variables, the lens equation will then determine one of the two in terms of the other in such a way that outgoing geodesic will reach the observer. The inclusion of an energy dependence in the angle of deflection α renders this picture slightly more complex. For instance, the lens equation will now depend on 3 parameters, which can be chosen to to be (β, θ I , E) or (β, p r , p φ ) or any other equivalent combination, with one of the three fixed in terms of the other two by the equation itself. For a monochromatic and spherical source of energy E, fixed at a position β, emitting a beam with a given impact parameter b respect to the deflector, the lens equation may not have a real solution, since the deflector may disperse the beam in such a way that it will never reach the observer. For a fixed spherical source which emits photons or neutrinos of any energy, one can look for solution in the reduced variables b, E. Being b related to the primary and secondary images θ I± , the beam that reaches the observer will be characterized by a unique energy E, assuming that the images are detected at angular positions θ I± . The argument above can be repeated by using any triple combination of independent kinematic variables among those mentioned above.
Having clarified this point, we now move to a description of the actual implementation of the lens equation is this extended framework. The angular location of the image θ I and the impact parameter are related in the geometry of the lens by Eq. (74), and this allows to search for solutions of the lens equation (73) where the ellipsis refer to the extra logarithmic contributions present in Eq.(44). The expression above is known analytically if we manage to solve explicitly the semiclassical equation (36), otherwise it has to be found by a numerical fit. However, it is clear that the ansatz for the fit has, in any case, to coincide with Eqs.
which is an obvious generalization of (64), the latter being valid only in the classical GR case. As we are going to illustrate below, (79) can be studied numerically for several geometrical configurations, which are obtained by varying the lensing parameters D LS and D OL . A similar approach can be followed for the VE or for any other classical lens equation. The insertion of α(θ I , E) given by (44) into (73) generates the radiative lens equation which takes into account also the quantum corrections and is now, on the contrary of (73), energy dependent. At this point it is clear that all the lens observables, such as magnifications, shears, light curves of microlensing etc. descend rather directly by this general prescription. For instance, we can determine for the Virbadhra-Ellis lens the expression for the magnification using the where α ≡ ∂α/∂θ I . As clear from Eqs. (80) and (81), both equations are very involved, although they can be investigated very accurately at numerical level. It is also possible to discuss the analytical form of the solutions within the formalism of the 1/b n expansion. In fact, we are entitled to expand all the observables of the nonlinear lens in the angular deflection α, and work at a certain level of accuracy in the angular parameters. In this work, however, we prefer to proceed with a direct numerical analysis of the full equations, both for the thin and for the VE lens, leaving the discussion of the explicit solutions to a future work.

Numerical analysis for neutrino and photon lensing
The analysis that we present in this section, especially in the neutrino case, is of exploratory nature. It has the aim to test the consistency of the theoretical approach presented in the previous sections, rather than being an explicit proposal for the detection of such effects. We start investigating the behaviour of the solutions for the VE lens equation, which are shown in Figs. 10, by plotting the angular position of the source (β) as a function of the location of the image θ I , for neutrino beams. In panel (a) we have chosen distances between lens/source and observer of galactic size. The branch of the solution with θ I > 0 describes a primary image, while for θ I < 0 the points (β, θ I ) on the curve describe configurations corresponding to secondary images. One can test the energy dependence of the solution by varying the energy of the original beam in the lens equation (80), which, for this specific geometry, is clearly unnoticeable. In general, in fact, specific geometries (D OL , D LS ) select impact parameters b h in the beam path which are quite large. In this case the impact parameter turns out to be pretty large (b h ∼ 10 5 ), corresponding to very weak deflections, and causes a superposition between the two curves, the one describing the classical GR solution of (73), and the semiclassical one, obtained by solving (80). The curves intersect the θ I (image) axis in two opposite points (θ I = ±θ E ), giving rise to the Einstein ring, which are obtained for β = 0, i.e. for  We show in Fig. 11 the solutions for the Sagittarius configuration. In panels (a) and (b) the solutions that we identify correspond to secondary images obtained for very small values of θ I , of the order of a milliarcsecond.  These are the only configurations which guarantee close encounters between the cosmic ray beam and the black hole, with 20 < b h < 100. As the distance D LS gets reduced, the equation has solutions with values of θ I which define primary images, being β > 0 (panel c). At the same time, as shown in panel (c), the angular position of the source moves from β < 0 to β > 0. A primary image is shown to form when D LS is 1 milliparsec (mpc). Simple considerations show that such lensing configurations are not unreasonable. For instance, for the supermassive black hole that we consider (with M=4.31 × 10 6 M ), this distance is of the order of 5 × 10 3 r s , with r s denoting its Schwarzschild radius (∼ 6 × 10 6 km). Being the center of our galaxy rather densely populated by massive compact sources, one could envisage a distribution of these covering a large array of possible distances from the center of the black hole. For instance, it has been found that stars may orbit the supermassive galactic black hole with orbital periods even of the order of T ∼ 11.4 ys, corresponding to orbital distances (R) from its center as close as few astronomical units (AU), R = T 2/3 ∼ 5 AU, i.e. R ∼ 125 r s . While such distances may correspond to realistic lensing configurations, the lensing resolutions of these specific events, which is of the order of a milliarcsecond, remains a challenging aspect of these studies. This is due to the strong quasi-alignment required between source, lens and observer along the optical axis, which might be difficult to measure. The energy dependence of the lensing configuration, extracted from Eq. (80), is illustrated by plotting the solutions for several values of the initial energy of the neutrino beam, corresponding to 1 GeV, 1 TeV and 1 PeV respectively. A comparison with the classical GR solution is included. Differences among the various predictions for the position of the source can be as large as 15%.
The analysis for the magnification of this lensing configuration is presented in Fig. 12. Panels (a) and (b) show the strong suppression of the secondary images identified in Fig. 11 (a) and (b). The magnification of the secondary and primary images of Fig. 11 (c) is shown in Fig. 12 (c). The µ < 0 and µ > 0 regions, in this figure, are separated by asymptotes for θ I = θ E , the Einstein angle of the lens. These regions in µ correspond to the secondary (β < 0) and (β > 0) branches of the panel (c) in Fig. 12 and therefore refer to secondary and primary pictures respectively. The dependence on the energy of the incoming neutrino beam appears in the form of 3 displaced curves and respective asymptotes. The three vertical asymptotes therefore characterize the dependence of the Einstein radius on the energy. The analysis is repeated in Figs. 13 and 14 for lensing events detected on Earth from the galactic center of the Andromeda galaxy. In this case as before, we vary the distance between the lens and the source, with D LS = 1 kpc, 1 pc and 1 milliparsec in panels (a), (b) and (c) respectively. The first two plots correspond to secondary images while panel (c) describes a primary image solution of the lens equation. The corresponding plots for the  Fig. 15 (c), corresponding to one primary and one secondary image, is associated with the magnification given in Fig. 16 (c). In these two sets of figures the energy dependence of the result is quite small. Finally, the numerical result for the photon lensing case, generated by the supermassive black hole in the Andromeda galactic center is discussed in Figs. 17 and 18. While the two secondary images found in 17 (a) and (b) give suppressed magnifications, the solution in panel (c) corresponds to a primary image. The corresponding magnification, shown in Fig. 18 has, obviously, a single branch, with an Einstein radius located at approximately 7 × 10 −5 arcsec.

Post-Newtonian corrections: the case of primordial black holes
We have seen in the previous sections that the b h (α) expression for the deflection does not suffer from any apparent divergence (from the gravity or external field side) due to well-defined structure of the Newtonian cross section. The expression given in (33), in fact, is similar to the ordinary Rutherford scattering encountered in electrodynamics.
The dependence of the resulting cross section on the scale GM/c 2 , the Schwarzschild radius, manifests as an overall dimensionful constant. Therefore, the inclusion of the electroweak corrections -and the logarithmic dependence on the energy of the terms in the expansion that follows -do not appear in combination with the macroscopic scale r s . This allows, in principle, an extension of the perturbative computation up to any order in the electroweak coupling constant α w . It is also clear that this result is expected to be valid for any renormalizable field theoretical model, when combined with an external static gravitational field of Coulomb type, as in the case of the Newtonian limit of GR. From now on, we will be using the notation nPN to indicate the (post-Newtonian) order in the potential at which we expand the Schwarzschild metric. For instance, contributions of a certain nPN order involve corrections in the external field proportional to Φ n+1 , with 0PN denoting the ordinary (lowest order) Newtonian (i.e. zeroth post- Newtonian) contributions proportional to Φ, as given in Eq. (32). The inclusion of the higher order corrections in the external potential modifies this simple picture due 1) to the need of introducing a cutoff regulator in the computation of the Fourier transform of the higher powers of the Newtonian potential and 2) to the presence of the Schwarzschild radius r s in the actual expansion. These features emerge already at the first post-Newtonian order (1PN) for an uncharged black hole and at order 0PN for the Reissner-Nordstrom (RN) metric (charged black hole). Both points 1) and 2) are, in a way, expected, since the microscopic expression for the transition matrix element given by (26), in fact, cannot be extrapolated to the case of a macroscopic source, with the presence of a macroscopic scale such as the black hole horizon. This seems to indicate that the success of the Newtonian approximation is essentially due to the rescaling of r s found in the expression of the cross section, which is a feature of this specific order, and is therefore limited to a 1/r potential. It is then natural to ask if there is any other realistic case in which the post-Newtonian corrections can be included in an analysis of this type.
Obviously, the answer is affirmative, as far as we require that r s is microscopic and that the energy of the beam, which is an independent variable of a scattering event, is at most of the order of 1/r s . Under these conditions, we are then allowed to extend our analysis through higher orders in Φ, with scatterings in which the dimensionless parameter r s q with q the impact parameter, is at most of O(1). This specific situation is encountered in the case of primordial black holes, where r s can be microscopic. We are going to illustrate this point in some detail, since it becomes relevant in the case of primordial black holes.

Post Newtonian contributions in classical GR
To illustrate this point we extend the expansion of the Schwarzschild metric at order 0PN given in (2). A similar expansion will be performed on the RN metric. For this purpose, it is convenient to perform a change of coordinates on the Schwarzschild metric in such a way that this takes an isotropic form. The radial change of coordinates is given by which allows to rewrite (82) as with Post-Newtonian (weak field) corrections can be obtained by an expansion of A and B taking M/ρ 1. Up to third order in Φ this is given by In the RN spacetime for a charged black hole the analysis runs similar. The interest in this metric is due to the fact that the lowest order potential, in this case, involves charge-dependent 1/r 2 contributions which, for an uncharged black hole, appear at first post-Newtonian order (1PN). The metric, in this case, is given by the expression with Q denoting the overall charge of the black hole. It has two concentric horizons which become degenerate in the maximally charged case. The two horizons are the solution of the equation The r i are the roots of the fourth order polynomial ordered so that r 1 > r 2 > r 3 > r 4 . The comparison between Schwarzschild and RN deflection angle is shown in Figure 19. The plots describe the behaviour of the angular deflection as a function of the impact parameter b h for a RN and Schwarzschild metric in the region with 10 < b h < 50 (top left) and 4 < b h < 10 (top right) for the maximally charged case. The differences tend to be very pronounced as we approach the horizon of the Schwarzschild metric.
As pointed out in [22] in the Schwarzschild case, the 1/b expansion for the deflection angle does not reproduce the photon sphere singularity of the Schwarzschild metric, which is achieved using the exact GR expression in terms of elliptic function given in (95), but it represents nevertheless an improvement respect to the 0P N order.
Expanding the RN metric in M/ρ 1 up to the third order, the 2P N approximation gives Inserting this expansion into the deflection integral, we can account in a systematic way of the 1/b corrections in the angle of deflection α The deflection (103) in the maximally charged case is given by the expression In the next subsection we are going to illustrate how the inclusion of these expansions at nPN order affects the computation of the quantum corrections to the angular deflection. The corrections are embodied in a geometric form factor whose expression is entirely controlled by the 1/b expansion.

Quantum effects at 2nd PN order
The inclusion of the PN corrections to the external background requires a recalculation of the cross section, with the inclusion of the additional terms in the fluctuation of the metric in momentum space. As usual we consider a static source, so that the metric is written as At leading order in the external field Φ both the timelike and the spacelike components are equal ( h 00 ≡ h ii ), while at higher orders they are expressed in terms of two form factors h 0 and h 1 which at higher order in the weak external field are given by where we have explicitly reinstated the dependence on the speed of light. Below we will conform to our previous notations in natural units, with c = 1.

• Neutrinos
The computation, at this stage, follows rather closely the approach of the previous sections, giving for the averaged squared matrix element in the neutrino case with T being the time of the transition, and the differential cross section We have denoted with dn f the density of final states in the transition amplitude, and with j i the incoming flux density. After integration over the final states, and using | p 1 | = | p 2 |, we obtain the expression dσ dΩ where we have introduced the gravitational form factor of the external source Notice the complete analogy between the corrections coming from a distributed source charge, for a potential scattering in quantum mechanics, and the gravity case. In the evaluation of F g in momentum space we are forced to introduce a cutoff Λ, being the Fourier transforms of the cubic contributions in Φ divergent. The singularity is generated by the integration around the region of r ∼ 0 in the Fourier transform of the potential. The relevant integrals in this case are given by with and with I 3 requiring a regularization with an ultraviolet cutoff in space (Λ) The choice of Λ is dictated by simple physical considerations. Given the fact that consistency of the expansion requires that r s q O(1), it is clear the appropriate choice in the regulator is given by the condition that this coincides with the Scwarzschild radius, i.e. Λ ∼ r s . Expressed in terms of the cutoff, we obtain for the geometric form factors the expressions where we have indicated with Ci the cosine integral function From the previous equations we obtain the cross section dσ dΩ which is valid at Born level and includes the weak field corrections up to the third order in Φ. In the expression of the cross sections, we use the subscript nPN, with n = 0, 1, 2 to indicate a n-th order expansion of the metric in the gravitational potential, while the superscripts ((0), (1) and so on) label the perturbative order in α w . The leading order cross section at order 2PN, for instance, takes the form dσ dΩ with where we have factorized the tree level result dσ/dΩ| (0) 0PN given in (33). The post-Newtonian form factor PN 2 (E, θ) induces an energy dependence of the cross section which is unrelated to the electroweak corrections. The analysis, in fact, can be extended at one loop in the electroweak theory. In this case, a lengthy computation gives the 2PN result dσ dΩ where we have inserted the one loop expression given in (34).
We can obtain an explicit solution of the corresponding semiclassical equation at order 1PN. Using the expression of the PN function at this order on the right hand side of (119) in order to generate the 1PN cross section at Born level, and solving the corresponding semiclassical equation (36) we obtain b 2 (0) 1PN = 4 (GM ) 2 − 1 + csc 2 α 2 + 2 ln sin At this point, we can invert Eq. (122) for α(b) obtaining for the tree level post Newtonian one. For the Reissner-Nordstrom geometry the situation is similar. The post-Newtonian form factor is then given by

• Photons
We can extend the analysis presented above for neutrinos to the photon case. Here the cross section takes the form dσ dΩ and, as in the neutrino case, we have dσ dΩ where we inserted the tree level cross section for the photon In the 0PN Newtonian limit, this cross section has been computed in [5], and takes the form where are the relevant electroweak form factors entering in the computation. In the previous equations the sum f k is over all Standard Model fermions, with m f k and Q f k their masses and charges. N c f k is 1 for leptons and 3 for quarks. Proceeding similarly to the neutrino case, the one loop cross section in the 2PN approximation takes the form dσ dΩ (1) with PN 2 given by (119), which can be inserted again in (36) and investigated numerically. Solving at order 1PN the analogous of (133), the solution of (36) gives b 2 (0) γ,1PN = 2 (GM ) 2 − 1 + 2 csc 2 α 2 + cos α + 8 ln sin E 2 (GM ) 4 π 2 11 + 12 cos α + cos 2α + 32 ln sin α 2 .
In the photon case the inversion formulae at orders 0PN and 1PN are given by and respectively.

Range of applicability
The structure of the one-loop 2PN result for neutrinos and photons shows the complete factorization between the quantum corrections and the background-dependent contributions. While the former are process dependent, the latter are general. Obviously, this result is not unexpected, and follows rather closely other typical similar cases in potential scattering in quantum mechanics. An example is the case of an electron scattering off a finite charge distribution characterized by a geometrical size R, where the finite size corrections are all contained in a geometric form factor. We recall that for a Coloumb interaction of the form V (r) = e 2 /r, the cross section is given in terms of the pointlike (p) amplitude with q = k− k and | q| = 2| k| sin θ/2 being the momentum transfer of the initial (final) momentum of the electron k ( k ) and charge e. The scattering angle is measured respect to the z-direction of the incoming electron. The charge of the static source has also been normalized to e. The corresponding cross section is given by and the modification induced by the size of the charge distribution (ρ(x)) is contained in  with For a uniform charge density, for instance, the geometrical form factor F ( q), which is the transform of the charge distribution, introduces a dimensionless variable qR in the cross section which is absent in the point-like (Coulomb) case, of the form The validity of the expression above is for qR 1, and the presence of the geometrical form factor is responsible for the fluctuations measured in the cross section as a result of the finite extension of the charged region.
In the analysis of the nPN corrections in gravity, the situation is clearly analogous, with the size of the horizon taking the role of the classical charge radius R. For ordinary (macroscopic) horizons (e.g. of a km size) r s ∼ GM invalidates the perturbative expansion due to the appearance of the GM E parameter in the expression of the post Newtonian factor PN (E, θ), which is small only if E ∼ 1/GM , a choice which is not relevant for our analysis, since it applies to particle beams whose energy is in the very far infrared. By imposing that the cutoff Λ coincides with the Schwarzschild radius (Λ ≡ r s ), one can immediately realize that the post-Newtonian expansion gets organized only in terms of this parameter (GM E). In the regions of strong deflections, which are those that concern our analysis, we can reasonably assume that y ≡ sin θ/2 ∼ O(1), if we use the GR prediction to estimate the bending angle. This allows to discuss the convergence of the PN expansion only in terms of the energy E of the incoming beam and of the size of the horizon. The analogous of the charge oscillations given by (141), in the gravitational case, are then uniquely related to the post-Newtonian function PN , and hence to the size of the parameter ΛE ∼ r s E which defines its expansion in powers of the gravitational potential. Assuming a small value of x ≡ GM E, we can indeed rewrite (119) via a small-x expansion, obtaining This expression can be used to investigate the range of applicability of these corrections in terms of the two factors appearing in x, the energy of the incoming beam and the size of the horizon of the gravitational source.
The requirement that such a parameter be small defines a unique range of applicability of such corrections in the quantum case. One possible application of the formalism which renders the PN corrections to the gravitational scattering quite sizable is in the context of primordial black holes [28], which have found a renewed interest in the current literature [29,30]. We just mention that primordial black holes (PBHs) have been considered a candidate component of dark matter since the 70's, and conjectured to have formed in the early universe by the gravitational collapse of large density fluctuations, with their abundances and sizes tightly constrained by various theoretical arguments. These range from Hawking radiation, which causes their decay to occur at a faster rate compared to a macroscopic black hole (of solar mass); bounds from their expected microlensing events; their influence on the CMB, just to mention a few [31]. For instance, the mechanism of thermal emission by Hawking radiation sets a significant lower bound on their mass (∼ 5 × 10 14 g), in order for them to survive up to the present age of the universe. This bound satisfies also other constraints, such as those coming from the possible interference of their decay with the formation of light elements at the nucleosynthesis time. With the launch of the FERMI gamma ray space telescope [11], the interest in this kind of component has found new widespread interest. The unprecedented sensitivity of its detector in the measurement of interferometric patterns generated by high energy cosmic rays (femtolensing events), such as Gamma Ray Bursts [10], has allowed to consider new bounds on their abundances [32]. The hypothesis of having PBHs as a dominant component of the dark matter of the universe provides remarkable constrains on their allowed mass values, except for a mass range 10 18 kg < M PBH < 10 23 kg, where it has been argued that they can still account for the majority of it. In other mass ranges several analyses indicate that the PBH fraction of dark matter cannot exceed 1% of the total [31]. PN corrections turn out to be significant for PBH in this mass range, due to the large variation induced on the PN function by the 1PN and 2PN terms. These may play a considerable role in a PBH mediated lensing event. We illustrate this behaviour by showing plots of the post Newtonian behaviour of the relevant expressions for lensing. In Fig. 20 (a) we plot the angular deflection as a function of the impact parameter for the Newtonian 0PN, and relative post Newtonian corrections. We have considered a primordial black hole with a mass of 10 −16 M , which carries a microscopic Schwarzschild radius (300 fm) and chosen E = 1/(GM ) = 0.6 MeV for the incoming photon beam. The impact of the corrections on the gravitational cross section are quite large, as one can easily figure out from panel (b), where we plot the factor PN as a function of the Schwarzschild radius for these compact massive objects, for b h ∼ 1 fm. For a more massive primordial black hole, with 200 < b h < 1000, the pattern is quite similar, as shown in panel (c). In both cases the post-Newtonian corrections appear to be significant, of the order of 15-20 % and could be included in a more accurate analysis of lensing for these types of dark matter candidate solutions.

Conclusions and perspectives
We have presented a discussion of neutrino lensing at 1-loop in the electroweak theory. In our approach the gravitational field is a static background, and the propagating matter fields are obtained by embedding the Standard Model Lagrangian on a curved spacetime, as discussed in previous works [14,17]. As in a previous study [5], also in our current case the field theoretical corrections to the gravitational deflection are in close agreement with the predictions of general relativity. The agreement holds both asymptotically, for very large distances from the center of the black hole, of the order of 10 6 horizon sizes (b h ), but also quite close to the photon sphere (∼ 20 b h ). In this respect, the similarity of the results for photons and neutrinos indicates the consistency of the semiclassical approach that we have implemented. Various types of lens equations have been formulated in the past using classical GR, and we have illustrated the modifications induced on their expressions by the inclusion of the suppressed 1/b n corrections in the impact parameter to the angular deflections. We have then developed a formalism which allows to include the semiclassical results, due to the radiative effects in the propagation of a photon or a neutrino, in a typical lensing event. We have considered both the case of a linear lens, where the dependence on the angle of the lensing geometry is linear, and the nonlinear case, taking as an example the Virbhadra-Ellis lens equation. Radiative and post-Newtonian effects induce a dependence of the angle of deflection with the appearance of extra 1/b n suppressed contributions and of extra logarithms of the impact parameter, that we have studied numerically for some realistic geometric configurations. In general, radiative effects are significant only for configurations of the source/lens/observer which involve small impact parameters in the deflection (b h ∼ 20), and require angular resolutions in the region of few milliarcseconds. Our results are valid for a Schwarzschild metric, considered both in the Newtonian and in the post-Newtonian approximation, but they can be extended to other metrics as well.
We have also discussed the consistency of the post-Newtonian approach. We have shown that such corrections can be consistently taken into account in the case of microscopic horizon sizes, such as primordial black holes. These corrections have been shown to factorize and be accounted for by a post-Newtonian function. Our analysis can be extended in several directions, from the case of Kerr-Newman metrics to the study of microlensing and Shapiro delays, and to dynamical gravity. We hope to return on some of these topics in a future work.
A Fluctuations to first order in the Newtonian potential The solution for the static massive source is obtained from the linearized equation where h ≡ h µν η µν , and can be rewritten as The external field h µν is obtained by convoluting the static source with the retarded propagator normalized as G R (x, y) = δ 4 (x − y).
The solution of Eq.(144) takes the form h ext µν (x) = κ d 4 y G R (x, y) S µν (y) with the EMT of the external localized source, defining S µν , given by For a compact source of mass M at rest at the origin, with P µ = (M, 0), we have which gives and where the field generated by a local (point-like, L) mass distribution has the typical 1/r (r ≡ | x|) behaviour. The fluctuations are normalized in such a way that h µν has mass dimension 1, as an ordinary bosonic field, with κ of mass dimension −1.
B 1/b n corrections to lensing for discrete and continuous mass distributions in GR The method of extracting the 1/b n corrections to the angular deflection can be extended to the case of a countinuos mass distribution in the lens plane, as shown in Fig. 7. For this purpose we can consider the case of a distributed lens with a surface density Σ( ξ) = ρ( ξ, z) dz.
In the Newtonian approximation we have the usual relation or, after the rescaling ξ = ξ 0 x,ˆ The deflection angle is and is the gradient of the lensing potential where we have introduced the critical surface density

F Deflection integral
The deflection integral for the GR solution can be recast in the form with x 0 being the point of closest approach between the deflector and the beam. It is an elliptic integral of first kind. In the indefinite form, the general expression of these types of integrals is given by with the polynomial at the denominator R(z) = c 0 z n + c 1 z n−1 + · · · + c n , being of degree n. For n = 3, which is the GR case, the Weierstrass form of (167) is obtained by introducing the roots (e 1 , e 2 , e 3 ) and given by By the transformation t = e 3 − e 1 − e 3 z 2 it can be brought to the Legendre form where k is its modulus. In the case of a finite integration, the form that is needed is also re-expressed as by a simple change of the integration variable. For (166) the corresponding roots are and the transformation (170) is given by which takes to the Legendre form with Keeping into account the finite integration region, (166) becomes with The boundaries of integration in (178) can be expressed in the Jacobi form by the substitution z → ζ = 1/z, obtaining with τ (x 0 ) = (∆(x0)) −1 e λ 2 (x 0 ) = (k 2 (x 0 )) −1 . After reabsorbing a factor in front of the integral in the definition of Σ(x 0 ), we obtain α(x 0 ) = −4Σ(x 0 )F (φ(x 0 ); λ(x 0 )) − π, with Σ(x 0 ) = φ(x 0 ) = arcsin(τ (x 0 )) (182b) τ (x 0 ) = 2

G Feynman Rules
We collect here all the Feynman rules involving an external gravitational field that have been used in this work. All the momenta are incoming • graviton -gauge boson -gauge boson vertex where V stands for the vector gauge bosons g, γ, Z and W .