††thanks: corresponding author

# Propagation and lensing of gravitational waves in Palatini $f(\hat{R})$
gravity

Sreekanth Harikumar<EMAIL_ADDRESS>National Centre for
Nuclear Research, Ludwika Pasteura 7, Warsaw 02-093, Poland  Laur Järv
<EMAIL_ADDRESS>Laboratory of Theoretical Physics, Institute of Physics,
University of Tartu, W. Ostwaldi 1, 50411 Tartu, Estonia  Margus Saal
<EMAIL_ADDRESS>Laboratory of Theoretical Physics, Institute of Physics,
University of Tartu, W. Ostwaldi 1, 50411 Tartu, Estonia  Aneta Wojnar
<EMAIL_ADDRESS>Department of Theoretical Physics & IPARCOS, Complutense
University of Madrid, E-28040, Madrid, Spain  Marek Biesiada
<EMAIL_ADDRESS>National Centre for Nuclear Research, Ludwika
Pasteura 7, Warsaw 02-093, Poland

###### Abstract

Accelerated expansion of the Universe prompted searches of modified gravity
theory beyond general relativity, instead of adding a mysterious dark energy
component with exotic physical properties. One such alternative gravity
approach is metric-affine Palatini $f(\hat{R})$ theory. By now routine
gravitational wave detections have opened a promising avenue of searching for
modified gravity effects. Future expected cases of strong lensing of
gravitational waves will enhance this opportunity further. In this paper, we
present a systematic study of the propagation and gravitational lensing of
gravitational waves in Palatini $f(\hat{R})$ gravity and compare it with
general relativity. Using the WKB approximation we explore the geometric-
optical limit of lensing and derive the corrections to the measured luminosity
distance of the gravitational source. In addition, we study the lensing by the
Singular Isothermal Sphere lens model and show that Palatini $f(\hat{R})$
modifies the lensing potential and hence the deflection angle. Then we show
that the lens model and chosen theory of gravity influences the rotation of
the gravitational wave polarization plane through the deflection angle. To be
more specific we discuss the $f(\hat{R})=\hat{R}+\alpha\hat{R}^{2}$ gravity
theory and find that the modifications comparing to general relativity are
negligible if the upper bound of $\alpha\sim 10^{9}\,$m2 suggested in the
literature is adopted. However, this bound is not firmly established and can
be updated in the future. Therefore, the results we obtained could be valuable
for further metric-affine gravity vs. general relativity tests involving
lensing of gravitational waves and comparison of luminosity distances measured
from electromagnetic and gravitational wave sources.

##  I Introduction

First detections LIGOScientific:2016aoc  followed by ongoing routine
observations LIGOScientific:2021djp  of gravitational waves (GW) opened a new
range of possibilities not only to investigate astrophysical phenomena
inherently hidden to the electromagnetic wave domain but also to test the
theory of gravity. In addition we have our first direction detection of GW
signal from the Pulsar Timing Array Antoniadis:2023rey ; NANOGrav:2023gor ;
Xu:2023wog ; Reardon:2023gzh . So far the results obtained only strengthened
our confidence in general relativity (GR) and severely constrained some of
alternative gravity theories scalartensor ; dead_ends ; Tessa_test . One of
the predictions of the GR is light bending by massive objects. This phenomenon
underlies the theory of gravitational lensing Gravitational_lenses1992 .
Hence, in the era of GW astronomy, the detection of GW signals lensed by
massive sources along the line of sight is highly anticipated. So far,
however, there is no observational evidence for the lensed GW signals in LVK
data Hannuksela:2019kle ; LIGOScientific:2021izm ; LIGOScientific:2023bwz .

The detection of lensed GW signals will open up unique opportunities for
precision cosmology Sereno:2011ty ; Cao:2019kgn ; Liao:2017ioi ; Li:2019rns ;
Hannuksela:2019kle ; Biesiada:2021pzo ; Grespan:2023cpa , detection of
Intermediate Mass Black Holes (IMBH)Lai:2018rto ; meena2023gravitational , as
well as detection of low mass halos and Primordial Black Holes Diego:2019rzc ;
Oguri:2020ldf . It can also lead to the tests of GR Ezquiaga:2020dao ;
Goyal:2020bkm ; Fan:2016swi ; Collett:2016dey ; Sharma:2023vme ;
Mishra:2023vzo  and other fundamental interactions Baker:2016reh . Next
generation of GW detectors like the Einstein Telescope (ET) Maggiore:2019uih ,
Cosmic Explorer (CE) Evans:2021gyd  and the space based detectors such as
DECIGO Kawamura:2020pcg  and LISA lisa_2017  are designed to observe distant
sources up to redshifts $z\sim 20$ and higher. Some of these sources located
at such large cosmological distances should be lensed. Currently the
literature discussing the predictions for gravitational wave lensing rate is
rich. For example, current generation detectors are expected to detect 1
lensed signal per year for A+ Li:2018prc . On the other hand, ET has a much
higher estimate of expected 50-100 strongly lensed events per year
Biesiada:2014kwa ; Yang_2019 ; Lilan2021  which is not surprising as the
accessible volume will be three orders of magnitude larger than in the current
generation detectors. Predictions for the DECIGO are a bit less optimistic due
to contamination by unresolved sources, yet 50 strongly lensed BH-BH systems
should be detected each year of DECIGO operation Piorkowska2021 .

In the much more familiar context of electromagentic lensing the approximation
of geometric optics (GO) is used, which is justified by the fact that the
wavelength is much smaller in comparison to the typical size of lens objects
(stars or galaxies). Regarding GWs, the frequency range currently probed by
ground-based detectors covers 10 Hz $<f<10$ kHz. Future space-borne detectors,
will probe 0.1 mHz $<f<100$ mHz – LISA and 1 mHz $<f<100$ Hz – DECIGO. This
corresponds to the GW wavelengths of $10^{4}$ m $<\lambda<10^{7}$ m in ground-
based and $10^{6}$ m $<\lambda<10^{12}$ m in space-borne detectors. In the
case of lenses whose Schwarzschild radii are comparable to $\lambda$ wave
optics (WO) approach is crucial Nakamura:1999uwi ; Takahashi:2003ix . Hence,
for ground-based detectors, WO should be used for the lenses less massive than
$10^{4}\;M_{\odot}$, while in the case of space-borne detectors, this upper
mass limit reaches $10^{9}\;M_{\odot}$.

Even though propagation and detection of GWs involves the weak field regime,
yet the structure of the waveform is created in violent processes in strong
gravity regime. Hence, the analysis of the waveforms offers a unique
opportunity to test competing theories of gravity. GR predicts only two
polarization modes known as the ‘plus’ and ‘cross’. However, other viable
metric theories of gravity typically predict more than these two polarizations
for a generic gravitational wave. In fact, the most general weak gravitational
wave is composed of six modes of polarization, expressible in terms of the six
“electric” components of the Riemann tensor Will_book ; Eardley73 ; Thorne73 .
Therefore the detection of extra polarizations has been considered as a
smoking gun for modified gravity. Yet, it has been shown that even in the GR
strong gravitational lensing may distort the GW waveform in in such a way that
extra polarizations appear Hou:2019wdg ; Dalang:2021qhu ; Cusin:2019rmt . This
is another motivation for studying GW lensing in alternative theories of
gravity. In what follows, we are going to examine the effect of modified
gravity on the lensed gravitational signal. To our knowledge, this is the
first work dealing with such a problem.

Here, we focus upon Palatini $f(\hat{R})$ gravity Olmo:2011uz , which is one
of the simplest extensions of GR. The theory does not introduce additional
degrees of freedom, and there are no instabilities of the kind found in the
metric $f(R)$ theories Sotiriou:2006sf . Palatini formulation has implications
for cosmology Stachowski:2016zio ; Szydlowski:2017uuy ; TeppaPannia:2018ale ;
Pinto:2018rfg ; Camera:2022myt ; Gialamas:2023flv  and astrophysical objects
Olmo:2019qsj ; Olmo:2019flu ; Wojnar:2020txr ; Wojnar:2020frr ; Benito:2021ywe
; Wojnar:2021xbr ; Kozak:2021ghd ; Sarmah:2021ule ; Kalita:2022trq ;
Kozak:2023axy , however in the vacuum and radiation dominated regions the
theory reduces to GR plus an effective cosmological constant, and easily
passes the Solar System tests Toniato:2019rrd ; Bonino:2020wps ; Hernandez-
Arboleda:2023abv . Gravitational waves in Palatini $f(\hat{R})$ carry only two
polarizations, like in the GR Alves:2009eg . Therefore if one is interested in
testing the Palatini $f(\hat{R})$ gravity with GWs, strong lensing effects are
of particular interest.

The paper is structured as follows. Section II introduces the essential
features of Palatini $f(\hat{R})$ gravity necessary for further discussion. In
the section III, we systematically derive the geometric optics limit in GR and
in Palatini $f(\hat{R})$ gravity. Then in section IV the propagation of GWs in
both theories is studied. In Section V we study the rotation of GW
polarization plane in GR and compare it with Palatini f(R) gravity. Then as an
example of lensing we discuss Singular Isothermal Sphere (SIS), which is a
reliable model for a galaxy scale lens, in the context of Palatini
$f(\hat{R})$ gravity. Differences between GR and Palatini formulation are
revealed and we conclude the paper in section VI.

##  II Palatini $f(\hat{R})$ gravity

Let us begin by briefly recalling some of the key features of Palatini
$f(\hat{R})$ gravity. The action of the theory can be written down as
Olmo:2011uz

$S[g,\hat{\Gamma},\psi_{m}]=\frac{1}{2\kappa^{2}}\int
d^{4}x\sqrt{|g|}f(\hat{R})+S_{\text{matter}}[g,\psi_{m}],$ (1)

where $\kappa^{2}=8\pi G/c^{4}$, $\hat{R}\equiv
g^{\mu\nu}\hat{R}_{\mu\nu}(\hat{\Gamma})$ is the Palatini curvature scalar,
and $|g|$ the absolute value of the determinant of the metric tensor. In
contrast with the usual metric approach, that is, instead of assuming that the
connection is the Levi-Civita one of the metric $g$, in the Palatini formalism
one considers an arbitrary torsion free (symmetric in lower indices)
connection $\hat{\Gamma}$ which is then used to construct the Riemann and
Ricci curvatures. The curvature scalar $\hat{R}$ is built from both geometric
structures, i.e., the metric and the independent of it, connection. On the
other hand, the matter action depends on the metric $g$ and matter fields
$\psi_{m}$ only,111That is, the matter fields are coupled to the metric only.
One may also couple them to the connection Iosifidis:2021nra , this broader
class is usually called metric-affine gravity. which guarantees that the
matter energy-momentum tensor $T_{\mu\nu}=-\frac{2}{\sqrt{|g|}}\frac{\delta
S_{m}}{\delta g_{\mu\nu}}$ obeys the usual conservation law Koivisto:2005yk ,

$\displaystyle\nabla_{\mu}T^{\mu\nu}=0,$ (2)

with respect to the Levi-Civita covariant derivative $\nabla$ of the metric
$g$.

The variation of the action (1) is taken with respect to both structures; the
metric one gives

$f^{\prime}(\hat{R})\hat{R}_{\mu\nu}-\frac{1}{2}f(\hat{R})g_{\mu\nu}=\kappa^{2}T_{\mu\nu},$
(3)

where prime is understood here as differentiating the function $f$ with
respect to its argument i.e. the curvature scalar $\hat{R}$. Contracting the
above equation with the metric $g_{\mu\nu}$ provides an algebraic relation
between the Palatini curvature scalar and trace of the energy-momentum tensor
$T\equiv g^{\mu\nu}T_{\mu\nu}$,

$f^{\prime}(\hat{R})\hat{R}-2f(\hat{R})=\kappa^{2}T.$ (4)

This feature allows to solve the above equations in some particular choices of
the function $f(\hat{R})$, providing $\hat{R}=\hat{R}(T)$.

On the other hand, the relation between the connection and the metric tensor
is given by the variation of (1) with respect to $\hat{\Gamma}$, which can be
written as

$\hat{\nabla}_{\beta}(\sqrt{|g|}f^{\prime}(\hat{R}(T))g^{\mu\nu})=0.$ (5)

The above covariant derivative is understood as the one defined by the
independent connection $\hat{\Gamma}$. Defining a new metric tensor

$\hat{g}_{\mu\nu}=f^{\prime}(\hat{R}(T))g_{\mu\nu}$ (6)

which is conformally related to the original metric $g_{\mu\nu}$, the equation
(5) can be expressed as

$\hat{\nabla}_{\beta}(\sqrt{|\hat{g}|}\hat{g}^{\mu\nu})=0.$ (7)

Therefore, the connection $\hat{\Gamma}$ happens to be the Levi-Civita one
with respect to the conformal metric $\hat{g}$ (6). It is related to the Levi-
Civita connection $\Gamma$ of the original metric $g$ by

$\hat{\Gamma}^{\alpha}_{\mu\nu}=\Gamma^{\alpha}_{\mu\nu}+\frac{1}{2f^{\prime}}[\delta^{\alpha}_{\nu}\partial_{\mu}f^{\prime}+\delta^{\alpha}_{\mu}\partial_{\nu}f^{\prime}-g_{\mu\nu}\partial^{\alpha}f^{\prime}].$
(8)

As $\hat{\Gamma}$ can be computed from the original metric and the trace of
the matter energy momentum, it may be considered an auxiliary field that can
be integrated out, with the implication that the degrees of freedom are
related to the metric tensor $g$ only. Therefore, in the end the connection is
not an independent dynamical quantity in the theory.

Equivalently, we can apply the conformal transformation (6) to the metric
field equation (3) to encapsulate the gravitational dynamics in terms of the
conformal metric $\hat{g}$,

$\hat{G}_{\mu\nu}=\kappa^{2}\hat{T}_{\mu\nu}-\frac{1}{2}\hat{g}_{\mu\nu}\hat{U}(f^{\prime})\,,$
(9)

where $\hat{G}_{\mu\nu}$ is the Einstein tensor computed from the Levi-Civita
connection $\hat{\Gamma}$ of $\hat{g}$, while
$\hat{T}_{\mu\nu}=(f^{\prime})^{-1}{T}_{\mu\nu}$ is the conformally
transformed matter energy-momentum tensor, and the effective “potential” is
given by

$\hat{U}=\frac{\hat{R}f^{\prime}(\hat{R})-f(\hat{R})}{f^{\prime}(\hat{R})^{2}}\,.$
(10)

The conformal energy-momentum does not obey the conformal conservation law,
but222To see the relation between the geometric objects in different frames,
see e.g. Stachowski:2016zio .

$\displaystyle\hat{\nabla}_{\mu}\hat{T}^{\mu\nu}$
$\displaystyle=\frac{1}{2\kappa^{2}}\hat{\nabla}^{\nu}\hat{U}(f^{\prime})\,.$
(11)

Let us notice that if $f^{\prime}(\hat{R}(T))$ is a constant, the connection
in (7) boils down to the Levi-Civita connection of the original metric $g$.
This happens when one deals with trace free matter energy-momentum, $T=0$, or
with a linear lagrangian. For the vanishing $T$ the field equations (3) and
(9) coincide and reduce to ferraris1993universal ; ferraris1994universality

$\hat{R}_{\mu\nu}-\Lambda(\hat{R}_{0})g_{\mu\nu}=0,$ (12)

where the “potential” $\hat{U}$ acts like an effective cosmological constant
that depends on the constant vacuum value of $\hat{R}$, that is

$\Lambda(\hat{R}_{0})\equiv f(\hat{R}_{0})/2f^{\prime}(\hat{R}_{0})$ (13)

with $\hat{R}_{0}=\hat{R}(T=0)$. Although its value does depend on a given
choice of the $f(\hat{R})$ model, it always leads to a unique solution for
each choice of $f(\hat{R})$. Therefore, in the vacuum or radiation environment
(for which also $T=0$), the Palatini $f(\hat{R})$ gravity reduces to GR with a
cosmological constant. One can argue that the theory passes Solar System tests
because of this property. Indeed, as examined in detail for analytic functions
$f(\hat{R})$ Toniato:2019rrd , the center-of-mass orbits are the same as in
GR, while modifications appearing in the terms related to the momentum and
energy in the Euler equation are not sensitive enough for the current
experiments which use the Solar System orbits.333However, let us mention that
these considerations assume Minkowski background. Moreover, when atomic level
experiments are in our reach schwartz2019post ; olmo2008hydrogen ;
olmo2007violation ; Kozak:2021ghd ; Wojnar:2022dvo ; Kozak:2023axy ;
Hernandez-Arboleda:2023abv , one may arrive at the required accuracy to
confront the theory against the data provided by the experiments performed in
the Solar System.

For the purposes of the GW physics one can neglect the cosmological constant
in GR naf2009gravitational . In this paper we assume that the Palatini
$f(\hat{R})$ effective cosmological constant $\Lambda(\hat{R}_{0})$ is
compliant with the observations of our universe and we can drop it from the GW
calculations, although in general there might exist specific choices of
$f(\hat{R})$ for which the quantity (13) is not negligible.

In the context of $f(\hat{R})=\hat{R}+\alpha\hat{R}^{2}$ gravity it has been
observed that the value of the parameter $\alpha$ is connected to the
curvature regime Olmo:2005zr . This relationship arises from the fact that the
Palatini curvature scalar is directly proportional to the trace of the energy-
momentum tensor, resulting in a similar dependence for its value. Furthermore,
when examining the weak-field limit analytically, it has been found that
$|\alpha|$ is approximately less than $2\times 10^{8}\rm m^{2}$ Olmo:2005zr .
However, due to uncertainties in microphysics, the experiments conducted
within the Solar System have not been able to establish any constraints on
these parameters Toniato:2019rrd . On the other hand, taking microphysical
properties into account, seismic data from Earth constrained the theory’s
parameter to $|\alpha|\lesssim 10^{9}\text{m}^{2}$ Kozak:2023axy ;
Kozak:2023ruu  while analysing neutron stars’ equations of state and combining
the results with the observational data put the bound to be $|\alpha|\lesssim
10^{6}$ m2 Lope-Oter:2023urz . The lower bound of the parameter
$\alpha\gtrapprox-7\times 10^{7}\text{ m}^{2}$ was obtained in Wojnar:2023bvv
ensuring microscopic stability of matter. Note that in the case of studies in
the non-relativistic limit, it was shown that only the quadratic term is
relevant; further terms entering the equations are of the sixth or higher
order Toniato:2019rrd , see the Poisson equation in the section IV. Similar to
the case of general relativity, none of the $f(\hat{R})$ models can adequately
explain the rotation curves of galaxies Hernandez-Arboleda:2022rim ;
Hernandez-Arboleda:2023abv . As a result, no constraints have been derived
from galaxy catalogues thus far. In contrast, the latest cosmological data
such SNIa and BAO data, provides bounds on the parameter $\alpha$ about forty
orders of magnitude bigger, that is, $|\alpha|\leq 10^{49}\,\text{m}^{2}$
Gomes:2023xzk .

In our investigation the most relevant feature of the family of the Palatini
theories is that the gravitational waves are propagating perturbations of the
background conformal metric $\hat{g}$. It is so since in the case of the
absence of the anisotropic stress, the tensor perturbations for the conformal
and the spacetime metrics coincide jimenez2015tensor ;
jimenez2017gravitational ; jimenez2018born . We will use that fact further in
this paper.

##  III Propagation of gravitational wave in the geometric optics regime

In order to be self-consistent and not to refer the reader to dispersed pieces
of literature, we present in this section basic formalism of GWs propagation
in the eikonal approximation. We start with the GR and then extend our
analysis to the Palatini gravity. Such an approach allows to trace
similarities and differences between these theories.

###  III.1 General relativity

In the study of gravitational waves over some arbitrary background, the metric
describing the background geometry (of larger amplitude) is considered to be
smooth over a typical characteristic length scale $L_{B}$. Gravitational waves
constitute perturbations of much smaller amplitude propagating on top of the
larger ones. The typical scale of spatial variation of this perturbation is
denoted by $\ell$, which obeys the condition $\ell\ll L_{B}$. Therefore, the
configuration can be described as

$g^{\text{total}}_{\mu\nu}=g_{\mu\nu}+{h}_{\mu\nu},$ (14)

where a natural distinction can be made in frequency space with $g_{\mu\nu}$
being a slowly varying background (or a low momenta component) while the
gravitational wave described by ${h}_{\mu\nu}$ is a high frequency
fluctuation. How fast the metric components vary is given by how their
derivatives change, that is, $\partial g\sim g/L_{B}$ and $\partial h\sim
h/\ell$. The separation of scales in frequency space is given by
$\mathfrak{f}\gg\mathfrak{f}_{B}$. The method of separation of scales and the
expansion based on that is called the short-wave expansion Maggiore:2007ulw .
Following this procedure we expand the Ricci tensor as

$R^{\text{total}}_{\mu\nu}=R_{\mu\nu}+R_{\mu\nu}^{[1]}+R_{\mu\nu}^{[2]}+\ldots$
(15)

where $R_{\mu\nu}$ is the background Ricci tensor containing only the slowly
varying component of the metric, while on the other hand, $R_{\mu\nu}^{[1]}$
and $R_{\mu\nu}^{[2]}$ contain terms which are linear and quadratic in the
metric perturbation ${h}_{\mu\nu}$, respectively. As $R_{\mu\nu}^{[1]}$ is
linear, it contains only higher frequency modes whereas $R_{\mu\nu}^{[2]}$ is
quadratic and contains both high and low frequency modes.

We can rewrite Einstein’s equations as

$R^{\text{total}}_{\mu\nu}=\kappa^{2}\left(T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T,\right)$
(16)

and the distinction of scales guarantees us the decomposition into two parts
Maggiore:2007ulw :

$R_{\mu\nu}=-\left(R_{\mu\nu}^{[2]}\right)^{\text{Low}}+\kappa^{2}\left(T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T\right)^{\text{Low}},$
(17)
$R_{\mu\nu}^{[1]}=-\left(R_{\mu\nu}^{[2]}\right)^{\text{High}}+\kappa^{2}\left(T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T\right)^{\text{High}}.$
(18)

It is the equation (17) that determines the dynamics of the background metric
$g{{}_{\mu\nu}}$ and the energy momentum tensor associated with the high
frequency mode $h_{\mu\nu}$. The order of magnitude of the high frequency
field equations suggests that $R_{\mu\nu}^{[1]}\sim\partial^{2}h\sim
h/\ell^{2}$ while $R_{\mu\nu}^{[2]}\sim\partial^{2}h^{2}\sim h^{2}/\ell^{2}$.
Therefore, $\left(R_{\mu\nu}^{[2]}\right)^{\text{High}}$ is quadratic in the
perturbation and negligible with respect to $R_{\mu\nu}^{[1]}$ and thus can be
neglected in the linear order equations. Similarly, the analogous estimation
applied to the right hand side of (18) provides that $T_{\mu\nu}$ is a smooth
macroscopic external matter and its high frequency contribution can arise only
via its dependence on $g_{\mu\nu}$ and via the term $g_{\mu\nu}T$. Since it
has the order $\mathcal{O}(h)$,

$\left(T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T\right)^{\text{High}}\sim\mathcal{O}(\frac{h}{L_{B}^{2}}).$
(19)

In comparison to $R_{\mu\nu}^{[1]}\sim h/\ell^{2}$ we see that (19) is smaller
than $R_{\mu\nu}^{[1]}$ by a factor of $\mathcal{O}(\ell^{2}/L_{B}^{2})$. In
summary, the low frequency part of the Einstein equations describes the effect
of external matter and GWs which affect the curvature of the background, while
the high frequency part reduces to the wave equation in this background

$\Box h_{\mu\nu}-2h_{\alpha\beta}R^{\alpha}{{}_{\mu\nu}}^{\beta}=0.$ (20)

where $\Box=\nabla_{\mu}\nabla^{\mu}$ as well as
$R^{\alpha}{{}_{\mu\nu}}^{\beta}$ are computed from the Levi-Civita connection
of the background metric $g_{\mu\nu}$. The above equation can be studied by
the Wentzel-Kramers-Brillouin (WKB) approximation444Also known as the eikonal
approximation or stationary phase approximation in the literature..

In GW theory we are mainly interested in plane waves. Although in a curved
background it is difficult to describe waves as exactly plane, we can still
consider them as nearly plane in comparison to the overall curvature for a
typical small wavelength. Therefore, the eikonal ansatz can be used to study
the effects of background on the propagation of GWs in geometric optics limit
and beyond. Then, if the information about the GW amplitude and polarization
is carried by the complex symmetric tensor denoted as $\xi_{\mu\nu}$, the high
frequency fluctuation can be written in the following form

$h_{\mu\nu}=\mathfrak{Re}\Big{\\{}\Big{[}\xi^{(0)}_{\mu\nu}+\epsilon\xi^{(1)}_{\mu\nu}+\epsilon^{2}\xi^{(2)}_{\mu\nu}+....\Big{]}e^{i\Phi/\epsilon}\Big{\\}},$
(21)

where $\Phi(x)$ is a real scalar function that depends on the coordinates and
defines the phase of the GWs while $\epsilon$ is an expansion book-keeping
parameter. The geometric optics limit corresponds to $\epsilon\rightarrow 0$.
As a consequence of this ansatz, independently of the systematic expansion in
${h}_{\mu\nu}$ there exists a derivative expansion of high frequency modes
given by the parameter $\epsilon=\ell/L_{B}$. The hierarchy of the terms are
as follows $\epsilon\gg h\epsilon\gg h\epsilon^{2}...\gg h^{2}$. Therefore,
each derivative of the high frequency fluctuation $h_{\mu\nu}$ of (21) will
collect a factor $1/\epsilon$, allowing to separate out geometric optics and
beyond geometric optics terms. To study the evolution of the fluctuations up
to a required order, one applies the eikonal ansatz (21) in the wave equation
(20)

$\displaystyle
e^{i\Phi/\epsilon}\Big{\\{}\frac{1}{\epsilon^{2}}[-k^{\beta}k_{\beta}\xi^{(0)}_{\mu\nu}]$
(22)
$\displaystyle+\frac{1}{\epsilon}[i(\nabla_{\beta}k^{\beta}\xi^{(0)}_{\mu\nu}+k^{\beta}\nabla_{\beta}\xi^{(0)}_{\mu\nu})-k^{\beta}k_{\beta}\xi^{(1)}_{\mu\nu}]+\mathcal{O}(\epsilon^{0})\Big{\\}}=0,$

where we have defined the wave vector $k^{\mu}$ as the gradient of $\Phi$,
that is, $k^{\mu}=g^{\mu\nu}\partial_{\nu}\Phi$. Applying the ansatz to the
Hilbert gauge condition $\nabla^{\mu}h_{\mu\nu}=0$ yields

$\displaystyle
e^{i\Phi/\epsilon}\Big{\\{}\frac{1}{\epsilon}[ik^{\mu}\xi^{(0)}_{\mu\nu}]+\mathcal{O}(\epsilon^{0})\Big{\\}}=0\,.$
(23)

Considering the leading $\mathcal{O}(\epsilon^{-1})$ term in (23), and the
leading order $\mathcal{O}(\epsilon^{-2})$ term with the next to the leading
order $\mathcal{O}(\epsilon^{-1})$ term in (22) defines the geometric optics
limit. Beyond the optical limit terms $\mathcal{O}(\epsilon^{0})$ were
considered in Dalang:2021qhu ; Cusin:2019rmt .

Therefore, we see from (23) that at the leading order the GW polarization is
transverse to the direction of propagation in the geometric optics limit,

$k^{\mu}\xi^{(0)}_{\mu\nu}=0.$ (24)

Similarly, the leading order term $\mathcal{O}(\epsilon^{-2})$ in the
expansion (22) provides us with the information that the wave vector is a null
vector field, $k_{\beta}k^{\beta}=0$, hence the gravitational waves travel at
the speed of light. Taking a covariant derivative of this relation and making
use of the fact that the wave vector is a gradient of the phase function
$\Phi$ gives us the geodesic equation

$k^{\mu}\nabla_{\mu}k_{\nu}=0,$ (25)

which tells us that the gravitational wave vector is parallel propagated along
the null geodesic of the background metric $g_{\mu\nu}$. It is easy to see
that by defining $k^{\mu}=dx^{\mu}/d\lambda$, where $\lambda$ is an affine
parameter along the geodesic, and substituting it in the equation (25) we get
the geodesic equation in the familiar form

$\frac{d^{2}x^{\beta}}{d\lambda^{2}}+\Gamma^{\beta}_{\alpha\mu}\frac{dx^{\alpha}}{d\lambda}\frac{dx^{\mu}}{d\lambda}=0.$
(26)

where $\Gamma^{\beta}_{\alpha\mu}$ are the coefficients of the Levi-Civita
connection of the background metric $g_{\mu\nu}$. The next to leading order
provides the following relation

$2k^{\beta}\nabla_{\beta}\xi^{(0)}_{\mu\nu}+\nabla_{\beta}k^{\beta}\xi^{(0)}_{\mu\nu}=0$
(27)

In order to understand this equation better, let us rewrite the tensor
$\xi^{(0)}_{\mu\nu}$ as

$\xi^{(0)}_{\mu\nu}=\mathcal{A}\mathcal{A}_{\mu\nu},$ (28)

where $\mathcal{A}$ is the amplitude, defined as
$\mathcal{A}=\sqrt{\xi^{*}_{\mu\nu}\xi^{\mu\nu}}$ where $\xi^{*}_{\mu\nu}$ is
the complex conjugate of $\xi^{\mu\nu}$ , and $\mathcal{A}_{\mu\nu}$ is the
normalized polarization tensor. Using this decomposition and the gauge
condition in the equation (27), we obtain the following relations for the
amplitude and the polarization tensor:

$\displaystyle\nabla_{\rho}(k^{\rho}\mathcal{A}^{2})$ $\displaystyle=0,$ (29)
$\displaystyle k^{\alpha}\nabla_{\alpha}\mathcal{A}_{\mu\nu}$
$\displaystyle=0.$ (30)

The equation (30) suggests that in the geometric optics limit the polarization
tensor is parallel transported along the gravitational wave direction
$k^{\mu}$. Moreover, the equation (29) can be written as a conservation
equation. Let us define the momentum of the gravitons as $P^{\mu}=\hbar
k^{\mu}$. Then, the graviton number density along the geodesic bundle is given
by $N^{\mu}=\frac{\mathcal{A}^{2}}{\hbar^{2}}P^{\mu}$. Applying it into (29)
leads to the conservation of graviton number density

$\nabla_{\mu}N^{\mu}=0$ (31)

in the geometric optics limit Gravitational_lenses1992 ; Maggiore:2007ulw  .
This arises as a consequence of absence of high frequency fields that
contribute to $\left(T_{\mu\nu}\right)^{\text{High}}$ in the Einstein field
equations Garoffolo:2019mna .

###  III.2 Palatini $f(\hat{R})$ gravity

To obtain all basic relations which govern the propagation of gravitational
waves in Palatini $f(\hat{R})$ gravity, we need to perturb the field equations
which we obtained in the section II. It is more convenient to start from the
conformal representation and use the perturbed form of the metric

$\hat{g}{}^{\text{total}}_{\mu\nu}=\hat{g}_{\mu\nu}+\hat{h}_{\mu\nu}$ (32)

where $\hat{g}{}_{\mu\nu}=f^{\prime}g{}_{\mu\nu}$ and
${\hat{h}}_{\mu\nu}=f^{\prime}h_{\mu\nu}$, cf. (6) and (14). The inverse
metric is defined by
$\hat{g}{}^{\mu\nu}_{\text{total}}=\hat{g}^{\mu\nu}-\hat{h}^{\mu\nu}$. As
already mentioned at the end of the section II, since we are dealing with no
anisotropic stress, the perturbations coincide, that is,
$\hat{h}_{\mu\nu}={h}_{\mu\nu}$ jimenez2015tensor ; jimenez2017gravitational ;
jimenez2018born .

Substituting the perturbation (32) in the equation (9), neglecting the
“potential” (effective cosmological constant) contribution and considering
only the high frequency part of the Palatini $f(\hat{R})$ gravity yields

$\hat{\Box}{}h_{\mu\nu}-2{h}_{\alpha\beta}\hat{R}^{\alpha}{}_{\mu\nu}{}^{\beta}=0\,,$
(33)

where $\hat{\Box}$ and $\hat{R}^{\alpha}{}_{\mu\nu}{}^{\beta}$ are computed
from the Levi-Civita connection $\hat{\Gamma}$ of the long wavelength
background metric $\hat{g}$. Notice that, in this case, the absence of matter
does not reduce Palatini theory to GR. The slowly varying part of the field
equations does not have a vanishing trace of the energy momentum tensor as it
happens in the GR. As before, applying the Eikonal ansatz one obtains the
following expansion similar to the GR one:

$\begin{split}e^{i\Phi/\epsilon}\Big{\\{}\frac{1}{\epsilon^{2}}[-\hat{k}^{\beta}\hat{k}_{\beta}\xi^{(0)}_{\mu\nu}]+\frac{1}{\epsilon}[i(\hat{\nabla}_{\beta}\hat{k}^{\beta}\xi^{(0)}_{\mu\nu}\\\
+\hat{k}^{\beta}\hat{\nabla}_{\beta}\xi^{(0)}_{\mu\nu})-\hat{k}^{\beta}\hat{k}_{\beta}\xi^{(1)}_{\mu\nu}]+\mathcal{O}(\epsilon^{0})\Big{\\}}=0\,,\end{split}$
(34)

where the covariant derivatives $\hat{\nabla}$ are defined with respect to the
connection $\hat{\Gamma}$. We see that the leading order provides that the GWs
in Palatini $f(\hat{R})$ gravity also follow null paths defined by the wave
vector $\hat{k}_{\mu}\hat{k}^{\mu}=0$. On the other hand, the geodesic
equation requires a bit of discussion. Taking the covariant derivative (with
respect to the connection $\hat{\Gamma}$) of this leading order term gives an
autoparallel

$\hat{k}^{\mu}\hat{\nabla}{{}_{\mu}}\hat{k}_{\nu}=0\,,$ (35)

which can be rewritten in a familiar form

$\frac{dx^{\beta}}{d\hat{\lambda}^{2}}+\hat{\Gamma}^{\beta}_{\alpha\mu}\frac{dx^{\alpha}}{d\hat{\lambda}}\frac{dx^{\mu}}{d\hat{\lambda}}=0$
(36)

by defining the wave vector as $\hat{k}^{\mu}=dx^{\mu}/d\hat{\lambda}$, where
$\hat{\lambda}$ is the affine parameter in the conformal frame. This is a
geodesic equation of the conformal metric $\hat{g}$ which we introduced as a
convenient calculation tool in the theory. In principle the geodesics of
$\hat{g}$ can be contrasted with the geodesics of the original metric $g$, and
it is the latter that are obeyed by the matter particles (since the matter
action (1) depends on the metric $g$). However, due to the conformal
invariance of the null directions, the null geodesics of $g$ and $\hat{g}$
coincide, and thus in the geometric optics limit the GWs in GR and Palatini
$f(\hat{R}$) gravity follow the same paths.555Consider the conformal relations
for the wave vectors $\hat{k}^{\mu}=(f^{\prime}(\hat{R}))^{-1}k^{\mu}$ and
$\hat{k}_{\mu}=k_{\mu}$ Stachowski:2016zio ; fleurry_thesis .

Although the next to leading order $\mathcal{O}(\epsilon^{-1})$ also seems to
be similar,

$2k^{\alpha}\hat{\nabla}_{\alpha}\xi^{(0)}_{\mu\nu}+\hat{\nabla}_{\alpha}\hat{k}^{\alpha}\xi^{(0)}_{\mu\nu}=0,$
(37)

let us decompose the wave tensor $\xi^{(0)}_{\mu\nu}$ as previously in order
to understand the physical consequences. The polarization tensor
$\mathcal{A}_{\mu\nu}$ in this theory is parallel propagated along the
direction of $\hat{k}^{\alpha}$,

$\hat{k}^{\alpha}\hat{\nabla}_{\alpha}\mathcal{A}_{\mu\nu}=0.$ (38)

However, since due to the conformal relationship the vectors
$\hat{k}^{\alpha}$ and $k^{\alpha}$ are parallel, the polarization tensor is
also parallel propagated along the light cone direction $k^{\alpha}$ of the
original metric $g$, like the light rays.

Moreover, we also obtain the information about the conservation of graviton
number density. It goes as follows: if we define the momentum of gravitons as
$\hat{P}=\hbar\hat{k}^{\mu}$, then graviton number density along a geodesic
ray bundle can be defined as
$\hat{N}^{\mu}=\frac{\mathcal{A}^{2}}{\hbar^{2}}\hat{P}^{\mu}$. Therefore, as
we have seen in the case GR (29), one can obtain a relation for the amplitude
as,

$2\hat{k}^{\alpha}\hat{\nabla}_{\alpha}\mathcal{A}+\mathcal{A}\hat{\nabla}_{\alpha}\hat{k}^{\alpha}=0\,.$
(39)

Rewriting the above equation as
$\hat{\nabla}_{\mu}(\hat{k}^{\mu}\mathcal{A}^{2})=0$ and making use of the
expression for graviton number density we obtain

$\displaystyle\hat{\nabla}_{\mu}\hat{N}^{\mu}=0\,.$ (40)

Note that the graviton number density conservation in Palatini $f(\hat{R})$
gravity holds for the quantity $\hat{N}$ and conformal metric $\hat{g}$ (in
terms of its Levi-Civita covariant derivative $\hat{\nabla}_{\mu}$). If we
rewrite these formulas for the quantity ${N}^{\mu}$ and the vector $k^{\mu}$
as in GR, the conservation of graviton number density turns out to be
violated. Namely, one has

$2{k}^{\alpha}{\nabla}_{\alpha}\mathcal{A}+\mathcal{A}{\nabla}_{\alpha}{k}^{\alpha}-\frac{1}{2}\mathcal{A}f^{\prime}(\hat{R})k^{\alpha}\partial_{\alpha}f^{\prime}(\hat{R})=0.$
(41)

which can be rewritten as

$\nabla_{\mu}N^{\mu}=-N^{\mu}\partial_{\mu}\ln f^{\prime}.$ (42)

This feature can be compared to the conservation of matter energy-momentum in
terms of $\nabla_{\mu}$, (2), but not in terms of the conformal
$\hat{\nabla}_{\mu}$, (11). For gravitational waves the roles of
$\nabla_{\mu}$ and $\hat{\nabla}_{\mu}$ have switched, which is one of the key
features of Palatini $f(\hat{R})$ gravity jimenez2015tensor ; jimenez2018born
. Later in Sec. IV.5 we see that the graviton number nonconservation can be
attributed to the change in amplitude. Such effects are also observed in
scalar-tensor theories Garoffolo:2019mna  due to the presence of a high
frequency scalar fields. Further, the non-conservation has consequences
related to the observable such as GW luminosity distance Deffayet:2007kf ;
Saltas:2014dha ; Nishizawa:2017nef ; Belgacem:2017ihm .

##  IV Gravitational wave amplitude evolution

Passing by distributions of matter can affect the amplitude and polarisation
of a gravitational wave, as well as deflect the direction of the rays. We will
focus on the amplitude first.

###  IV.1 Newman-Penrose tetrad in both theories

The analysis of the polarization states and the evolution of amplitude can be
made more simple by introducing a Newman-Penrose (NP) tetrad Newman:1961qr ;
Chandrasekhar:1985kt  of null co-vectors

$e^{A}_{\mu}=\\{k_{\mu},m_{\mu},l_{\mu},n_{\mu}\\}\,,$ (43)

which are parallel transported along the null geodesics. The vectors $k^{\mu}$
and $n^{\mu}$ are real while $m^{\mu}$ and $l^{\mu}$ are complex conjugate
pairs whose the only non-vanishing contractions satisfy the relations:

$m_{\mu}l^{\mu}=1,\hskip 14.22636ptk_{\mu}n^{\mu}=-1.$ (44)

It can be shown that they can be used to express the spacetime metric as

$g_{\mu\nu}=m_{\mu}l_{\nu}+m_{\nu}l_{\mu}-n_{\mu}k_{\nu}-n_{\nu}k_{\mu}.$ (45)

The GW is represented by a second rank symmetric tensor. Therefore, ten
independent components of the tensor can be rewritten in terms of a symmetric
basis constructed from the null tetrad Cusin:2019rmt

$\Theta^{AB}_{\mu\nu}=\frac{1}{2}(A_{\mu}B_{\nu}+A_{\nu}B_{\mu}),$ (46)

where $A_{\mu},B_{\mu}$ have the same structure as $e_{\mu}$ in (43). The
tensor $\xi_{\mu\nu}$ can be then expanded in this basis as follows

$\begin{split}&\xi_{\mu\nu}=\mathcal{C}_{kk}\Theta^{kk}_{\mu\nu}+\mathcal{C}_{ll}\Theta^{ll}_{\mu\nu}+\mathcal{C}_{mm}\Theta^{mm}_{\mu\nu}+\mathcal{C}_{nn}\Theta^{nn}_{\mu\nu}+\mathcal{C}_{kl}\Theta^{kl}_{\mu\nu}\\\
&\quad+\mathcal{C}_{km}\Theta^{km}_{\mu\nu}+\mathcal{C}_{kn}\Theta^{kn}_{\mu\nu}+\mathcal{C}_{ml}\Theta^{ml}_{\mu\nu}+\mathcal{C}_{nl}\Theta^{ml}_{\mu\nu}+\mathcal{C}_{mn}\Theta^{mn}_{\mu\nu},\end{split}$
(47)

where $\mathcal{C}_{AB}$ are complex coefficients. In the geometric optics
limit, after using all gauge freedom this tensor reduces to the following form

$\xi^{(0)}_{\mu\nu}=\mathcal{C}^{(0)}_{ll}l_{\mu}l_{\nu}+\mathcal{C}^{(0)}_{mm}m_{\mu}m_{\nu},$
(48)

which reveals the two polarizations. To see this let us substitute
$l^{\mu}=(x^{\mu}+iy^{\mu})/\sqrt{2}$ and
$m^{\mu}=(x^{\mu}-iy^{\mu})/\sqrt{2}$, whereby (48) can be rewritten as

$\xi^{(0)}_{\mu\nu}=H_{+}\mathcal{A}^{+}_{\mu\nu}+H_{\times}\mathcal{A}^{\times}_{\mu\nu}\,.$
(49)

Here $\mathcal{A}^{+}_{\mu\nu}=x_{\mu}x_{\nu}-y_{\mu}y_{\nu}$ and
$\mathcal{A}^{\times}_{\mu\nu}=x_{\mu}y_{\nu}+y_{\nu}y_{\nu}$ are the two
polarization tensors in GR and $H_{+},H_{\times}$ are their corresponding
amplitudes. The complex amplitudes $\mathcal{C}_{AB}$ are related to the
polarizations amplitudes as follows

$H_{+}=\frac{1}{2}(\mathcal{C}^{(0)}_{ll}+\mathcal{C}^{(0)}_{mm})\;;\;H_{\times}=\frac{i}{2}(\mathcal{C}^{(0)}_{ll}-\mathcal{C}^{(0)}_{mm}).$
(50)

###  IV.2 General relativity

Gravitational waves in the geometrical optics limit are rays of null
congruences. The congruence is a set of parametrized curves such that only one
curve passes through each point. The behaviour of one null geodesic with
respect to another reference ray is governed by the separation vector
$\chi^{\mu}$. The evolution of $\chi^{\mu}$ depends on the background Riemann
tensor $R_{\mu\nu\rho\sigma}$ and obeys the geodesic deviation equation

$\frac{D^{2}\chi^{\mu}}{d\lambda^{2}}=R^{\mu}{}_{\nu\rho\sigma}k^{\nu}k^{\rho}\chi^{\sigma},$
(51)

where $\lambda$ denotes the affine parameter and $D$ is the directional
covariant derivative, that is, $k^{\mu}\nabla_{\mu}$ with respect to the the
null geodesics $k^{\mu}$. It can be written in the matrix form (see Appendix
(A)) and it is also known as the Sachs equation when projected on a two-
dimensional spacelike surface (see the discussion below).

The cross-section of the beam has a nonzero extension except from the point of
observation where the rays converge and the separation vector
$\chi(\lambda_{O})=0$, with $\lambda_{O}$ being the value of the affine
parameter at the observation point. Therefore, a given observer with four
velocity $u^{\mu}$ can also project the separation vector into a 2-dimensional
surface which is called a screen space. The screen space is perpendicular to
the line of the sight direction $d^{\mu}$ of the observer, with
$d^{\mu}d_{\mu}=1$ and $d^{\mu}u_{\mu}=0$. The screen projector is defined as

$Q^{\mu\nu}=g^{\mu\nu}+u^{\mu}u^{\nu}-d^{\mu}d^{\nu}\,.$ (52)

To study the morphology of the beam in the screen space by an observer with
$u^{\mu}$, let us consider the physical area $A$ of the cross section of the
beam on the screen space. It is given by the following expression

$A=\int
d\chi^{1}d\chi^{2}=\int\text{det}\mathcal{D}d\dot{\chi}^{1}_{O}d\dot{\chi}^{2}_{O},$
(53)

where $\mathcal{D}$ is the Jacobi matrix satisfying the equation (101), where
$\dot{\chi}^{1}_{O}$ denotes the derivative with respect to the affine
parameter.

Let us now consider an infinitesimal light beam. Then, the Jacobi matrix
$\mathcal{D}$ can be regarded as a constant and the evolution rate of area of
cross section of the beam (53) with respect to the affine parameter $\lambda$
can be expressed with the use of the expansion scalar $\theta$ (see its
definition in the Appendix A)

$\frac{\dot{A}}{A}=\frac{1}{\det\mathcal{D}}\frac{d(\det\mathcal{D})}{d\lambda}=\text{tr}(\mathcal{S})=2\theta.$

Therefore, the expansion scalar $\theta$ can be interpreted as the evolution
rate of the beam’s area:

$\theta=\frac{1}{D_{A}}\frac{dD_{A}}{d\lambda},$ (54)

where the angular diameter distance $D_{A}$ is proportional to $\sqrt{A}$
Fluerry_2015 . Now on, using the above relations along with amplitude
evolution equation (27) rewritten in terms of two polarization states in GR
(48), we get

$2k^{\mu}\nabla_{\mu}\mathcal{C}_{AB}+\nabla_{\mu}k^{\mu}\mathcal{C}_{AB}=0$
(55)

and

$\frac{d\ln\mathcal{C}_{AB}}{d\lambda}+\frac{d\ln D_{A}}{d\lambda}=0\,.$ (56)

Integrating this equation from the source to the observation point we get the
amplitude of polarization at an arbitrary position in terms of the affine
parameter $\lambda$

$\mathcal{C}^{\text{GR}}_{AB}(\lambda)=\frac{\mathcal{C}_{AB}(\lambda_{s})D(\lambda_{s})}{D(\lambda)}\,,$
(57)

where $\lambda_{s}$ is the affine parameter at the source position.

###  IV.3 Palatini $f(\hat{R})$ gravity

As mentioned in earlier sections of Palatini $f(\hat{R})$ gravity, the metric
$\hat{g}$ that is convenient to describe the gravitational dynamics in this
theory is related to the physical metric $g$ via the conformal factor
$f^{\prime}(\hat{R})$. Therefore all quantities that describe the morphology
of the beam, including the Sachs vector equations and expansion, carry an
imprint of this conformal factor. The screen projector is given by
$\hat{Q}_{\mu\nu}=f^{\prime}Q_{\mu\nu}$, this is due to re-scaling of line of
sight direction and the four velocity of the observer666Notice that $\hat{u}$
is not a physical observer - the observer vector field is related to the
metric $g$, and in our case she co-moves with the fluid given by the energy-
momentum tensor $T_{\mu\nu}$.

$\hat{d}_{\mu}=\sqrt{f^{\prime}}d_{\mu}\,,\qquad\hat{u}_{\nu}=\sqrt{f^{\prime}}u_{\nu}\,.$
(58)

As a consequence of this conformal transformation the expansion term (54) gets
an extra term which is related to the conformal factor and hence to the trace
of the energy momentum tensor via the trace equation (4). Therefore the Ricci
and Weyl lensing terms get modified in Palatini $f(\hat{R})$ gravity. The
expansion scalar in Palatini $f(\hat{R})$ gravity is related to the expansion
scalar in general relativity as

$\hat{\theta}=\frac{1}{2}\hat{\nabla}_{\alpha}\hat{k}^{\alpha}=\frac{1}{2}\frac{d\ln
f^{\prime}}{d\lambda}+\frac{1}{f^{\prime}}\theta$ (59)

where $\theta$ is related to the angular diameter distance (54).

Using the relation for the evolution of amplitude in Palatini $f(\hat{R})$
gravity given by (37) and making use of the new basis introduced in the
section (IV.1) and the gauge freedom, we obtain the following equation for the
evolution of amplitude

$2\hat{k}^{\mu}\hat{\nabla}_{\mu}\mathcal{C}_{AB}+\hat{\nabla}_{\mu}\hat{k}^{\mu}\mathcal{C}_{AB}=0$
(60)

where $\mathcal{C}_{AB}$ denotes the amplitude of two polarizations. It is to
be remembered that the wave direction is given by
$\hat{k}^{\mu}=\frac{k^{\mu}}{f^{\prime}}$ and the covariant derivative
$\hat{\nabla}$ is given by the connection $\hat{\Gamma}$ in Palatini
$f(\hat{R})$ gravity. Integrating the above equation from source to the
observation point and making use of (59) one gets the evolution of the
gravitational wave amplitude in this theory of gravity,

$\displaystyle\mathcal{C}^{\text{Palatini}}_{AB}$ (61)
$\displaystyle=\mathcal{C}^{\text{GR}}_{AB}(\lambda)\exp\Big{[}-\int_{\lambda_{s}}^{\lambda_{o}}\Big{(}2\frac{d\ln
f^{\prime}}{d\lambda}+\frac{1}{2}f^{\prime}\frac{d\ln
f^{\prime}}{d\lambda}\Big{)}d\lambda\Big{]},$

where $\mathcal{C}^{\text{GR}}_{AB}(\lambda)$ is given by (57) where
$\lambda_{o}$ and $\lambda_{s}$ are the value of the affine parameter at the
observer and source respectively. This is our main result. We observe that in
the geometric optics limit the two polarizations evolve independently, in the
same way as in GR. However, in Palatini $f(\hat{R})$ gravity we are dealing
with the additional term making the amplitude to decrease with the distance.
The exponential term is related to the mass distribution of the lens through
the conformal factor. This extra exponential factor modifies the amplitude of
GW detected, hence modifying the luminosity distance inferred from the signal
based on the GR prediction. The corrected luminosity distance in Palatini
theory takes the following form

$d^{\text{P}}_{L}=d^{\text{GR}}_{L}\exp\Big{[}\int_{\lambda_{s}}^{\lambda_{o}}\Big{(}2\frac{d\ln
f^{\prime}}{d\lambda}+\frac{1}{2}f^{\prime}\frac{d\ln
f^{\prime}}{d\lambda}\Big{)}d\lambda\Big{]}$ (62)

As mentioned in Sec. III, this is a consequence of non-conservation of
graviton number density. Therefore in the event of GW lensing in Palatini
$f(\hat{R})$ one observes a larger luminosity. Further, in the case of multi-
messenger GW event with the electromagnetic follow up observation, luminosity
distance measured in the electromagnetic window would be different from
measured using GW standard sirens.

###  IV.4 Nonrelativistic limit of Palatini $f(\hat{R})$ gravity

In what follows, we will consider a non-relativistic object described by the
Palatini $f(\hat{R})$ gravity, such as a star or galaxy. To do so, we need to
briefly discuss the non-relativistic limit of the theory. For an analytic
function $f(\hat{R})=\sum_{i=0}\alpha_{i}\hat{R}^{i}$, the Poisson equation is
given by Toniato:2019rrd ; Hernandez-Arboleda:2023abv

$\nabla^{2}\Phi\approx 4\pi G(\rho+2\alpha\nabla^{2}\rho),$ (63)

where $\alpha$ comes from the quadratic term of the Lagrangian (we neglect
cosmological constant). Therefore, in the non-relativistic objects, the only
correction to the GR equations which has a negligible effect is given by the
polynomial term

$f(\hat{R})=\hat{R}+\alpha\hat{R}^{2}.$ (64)

Inserting this functional into the trace equation (3) provides the analogous
to GR relation between the trace of the energy-momentum tensor and the
Palatini curvature scalar:

$\hat{R}=-\frac{8\pi G}{c^{4}}T,$ (65)

where $T=3p-c^{2}\rho$ if we consider a perfect-fluid description of matter.
Since we are interested in the non-relativistic objects, pressure, when it
appears together with the energy density, is negligible in such objects, that
is, $p/c^{2}\ll\rho$. Therefore, any modification introduced by Palatini
gravity in the Poisson equation (63) and further in the gravitational
potential, is signalized by the functions of the energy density $\rho(r)$
(with the parameter $\alpha$).

###  IV.5 Gravitational wave amplitude evolution

The above set of equations can be used to evaluate gravitational wave
amplitude evolution and hence the luminosity distance correction in Palatini
$f(\hat{R})$ gravity. From Eqs. (64), (65) the conformal factor becomes
$f^{\prime}=1+2\alpha\hat{R}=1+A\rho$, where $A=2\alpha\kappa^{2}c^{2}$. As a
first approximation one can rewrite the amplitude evolution expression (61) in
the form

$\mathcal{C}^{P}_{AB}=\mathcal{C}^{GR}_{AB}\exp\Big{[}-\frac{5A}{2}\int_{\rho_{s}}^{\rho_{o}}d\rho\Big{]},$
(66)

It is quite clear from (62) the above equation can be expressed as the
correction to luminosity distance. Any deviation from GR comes from the
exponential factor present in the above equation.

In order to see the effect of Palatini gravity on the amplitude evolution, we
need to consider some density profile which describes our lens. As Palatini
$f(\hat{R})$ gravity differs from GR only in the presence of matter, it is
evident from above that the most popular lens model considered in the
literature, that is, a point like mass in the vacuum Gravitational_lenses1992
, will not give any additional effect compared to GR. Because of that fact, we
need to consider an extended body with a given density profile that resides
asymmetrically on the path of the gravitational ray, whereby the entry and
exit values of the density differ. For instance, we may consider a GW source
close to the centre of a galaxy cluster and estimate the accumulated extra
effect as the wave exits the cluster and passes into the cosmic intergalactic
medium. Computing the terms in the exponential for NFW density profile
Navarro:1996gj  for an XCOP galaxy cluster sample Ettori:2018tus  A85 with a
radial distance of $b=3\text{ Mpc}$ from the cluster centre one obtains the
magnitude of the correction factor to be of the order $10^{-42}$ for
$\alpha\sim 10^{9}$. Analogously, assuming a Singular Isothermal Sphere (SIS)
for the same radial distance one estimates the correction to be of the order
$10^{-43}$. Therefore it can be seen that the correction to the GW amplitude
as well as to the luminosity distance during the wave propagation is
negligible and clearly beyond the sensitivity limit of the GW detectors. Out
of curiosity, we may also estimate the accumulated extra effect for an ancient
wave that has been travelling since early universe. Integrating over the
cosmic mean density from a redshift $z_{s}=3000$ till now, will have the
luminosity distance corrected by a factor of $10^{-33}$ which is still
negligible. However, if one uses the constraint for $\alpha\sim 10^{49}m^{2}$
from SNIa and BAO Gomes:2023xzk  the corrections for the GW luminosity
distance is $10^{-2}$ which is non-negligible.

##  V Gravitational wave deflection angle and rotation of polarization plane

So far we have discussed propagation of GW through an extended distribution of
matter characterized by some energy density $\rho$ distribution profile. In
this section we discuss how GW propagation is influenced by a compact matter
distribution, i.e. gravitational lensing.

###  V.1 Gravitational lensing

The subject of gravitational lensing has became a mature field of astrophysics
Gravitational_lenses1992 ; Meneghetti  having its own rich mathematical
formalism Petters  and numerous practical applications Treu . Referring the
reader to the literature just cited for a detailed and systematic
introduction, we give here a short, concise review of basic concepts of strong
lensing and notation used in sections below. Strong gravitational lensing
occurs whenever the source located at a distance777In cosmological scales
these distances have a meaning of angular diameter distances. $D_{S}$, a
deflecting mass – the lens at a distance of $D_{L}$ from the observer
(distance between the source and the lens being $D_{LS}$) are almost perfectly
aligned. Then instead of a single image of the source, multiple images can be
formed. This is possible, because the weak perturbation of the background
metric by the lens endow the space with the effective refractive index
$n=1-\frac{2\Phi}{c^{2}}$, where $\Phi$ is the Newtonian potential of the
lens.

In practice the size of the lens is negligible in comparison of the
characteristic distances of this optical system and it is convenient to use
thin lens approximation i.e. to assume that all changes in propagation of null
geodesics from the source to the observer (with respect to propagation through
the background space undistorted by the lens) occurs at the locus of the
closest encounter with the lens. This suggests a natural choice of (local)
coordinate system: the line of sight to the center of the lens $z$ (optical
axis of the system), the lens plane - perpendicular to $z$ at lens location
and the source plane - perpendicular to $z$ at source location. Intersection
of $z$ axis with lens and source planes define the origin of coordinate
systems on these planes spanned by arbitrarily (but consistently) chosen
orthogonal vectors. Radius vectors (in physical [length] units) on the lens
and source planes are denoted $\vec{\xi}$ and $\vec{\eta}$, respectively.
Because observations on the sky (the celestial sphere) are expressed in terms
of angles, it is convenient to describe strong lensing in the following way.
The rays from the source at the intrinsic angular position $\vec{\beta}$, i.e.
lying on the source plane at a distance $\vec{\eta}=\vec{\beta}D_{S}$ from the
optical axis, impact the lens plane at the position
$\vec{\xi}=\vec{\theta}D_{L}$ and reach the observer because they are
deflected by the angle $\vec{\gamma}$. In consequence the observer registers
the source’s image at angular position $\vec{\theta}$ instead of intrinsic
$\vec{\beta}$, which is governed by the lens equation

$\vec{\beta}=\vec{\theta}-\vec{\gamma}(\vec{\theta})\,.$ (67)

It is convenient to choose some length scale $\xi_{0}$ on the lens plane (and
associated to it $\eta_{0}=\xi_{0}D_{S}/D_{L}$ on the source plane), and use
dimensionless vectors $\vec{x}=\vec{\xi}/\xi_{0}$,
$\vec{y}=\vec{\eta}/\eta_{0}$. In this case the lens equation reads
$\vec{y}=\vec{x}-\vec{\gamma}(\vec{x})$ and the deflection angle is
$\vec{\gamma}=\nabla_{\vec{x}}(\Psi(\vec{x}))$. The dimensionless 2D
deflection potential is defined as

$\Psi(\vec{\xi})=\frac{D_{L}D_{LS}}{\xi_{0}D_{S}}\frac{2}{c^{2}}\int\Phi(\vec{\xi},z)dz.$
(68)

Let us finally remark, that although choice of the length scale $\xi_{0}$ is
arbitrary, it is convenient to use for this purpose the so called Einstein
radius on the lens plane $\xi_{0}=\theta_{E}D_{L}$, where angular Einstein
radius $\theta_{E}$ satisfies $\theta_{E}=\vec{\gamma}(\theta_{E})$. It is
also useful to introduce the convergence

$\kappa(\vec{\xi})=\Sigma(\vec{\xi})/\Sigma_{cr}$ (69)

where

$\Sigma(\vec{\xi})=\int_{-\infty}^{\infty}\rho(\vec{\xi},z)dz$ (70)

is surface mass density of the lens and

$\Sigma_{cr}=\frac{c^{2}}{4\pi G}\frac{D_{S}}{D_{L}D_{LS}}$ (71)

is the so called critical surface density Gravitational_lenses1992 .
Dimensionless 2D deflection potential and convergence satisfy 2D Poisson
equation:

$\nabla^{2}_{\vec{x}}\Psi(\vec{x})=2\kappa(\vec{x}).$ (72)

###  V.2 General Relativity

To evaluate the corrections in polarization due to lensing, one needs to
perturb the tetrad $k^{\mu}=\bar{k}^{\mu}+\delta k^{\mu}$, where $\delta
k^{\mu}$ is the first order correction which satisfies the perturbed geodesic
equation fleurry_thesis :

$\bar{k}^{\nu}\nabla_{\nu}\delta
k^{\mu}=-\delta\Gamma^{\mu}_{\nu\rho}\bar{k}^{\nu}\bar{k}^{\rho}\,.$ (73)

The Christoffel symbols are evaluated for the following weak-field metric

$ds^{2}=-(1+2\Phi)dt^{2}+(1-2\Phi)(dx^{2}+dy^{2}+dz^{2}),$ (74)

where $\Phi$ is the gravitational potential of the lens888We have already used
the fact that gravitational potentials coincide in the case of GR; see
Appendix B for the modified gravity case.. Solving (73) for the non-vanishing
contributions of the Christoffel components arising from the (74), one obtains

$\delta
x^{\mu}=\frac{1}{2}\Big{(}\vec{\gamma}\cdot\vec{x},-(\vec{\gamma}\cdot\vec{x})\vec{k}\Big{)}\,,$
(75) $\delta
y^{\mu}=\frac{1}{2}\Big{(}\vec{\gamma}\cdot\vec{y},-(\vec{\gamma}\cdot\vec{y})\vec{k}\Big{)}\,,$
(76)

where $\vec{\gamma}=\nabla_{x}(\Phi(x))$ is the deflection angle due to
lensing. The vectors $\vec{x},\vec{y}$ are the spatial parts of
$x^{\mu},\,y^{\nu}$ introduced in IV.1, and $\vec{k}$ denotes the direction of
propagation of the GW in space. Now, one can use the basis described in IV.1
to write the wave tensor in a more familiar form. Therefore, the wave tensor
(49) can be rewritten by changing the wave vector as $x^{\mu}\rightarrow
x^{\mu}+\delta x^{\mu}$, $y^{\mu}\rightarrow y^{\mu}+\delta y^{\mu}$ arsing
from the presence of the lens, where $\delta x$ and $\delta y$ are given by
(75) and (76), respectively. The polarization tensor undergoes a rotation and
the wave tensor becomes

$\xi^{(0)}_{\mu\nu}=H_{+}\tilde{\mathcal{A}}^{+}_{\mu\nu}+H_{\times}\tilde{\mathcal{A}}^{\times}_{\mu\nu},$
(77)

where $\tilde{\mathcal{A}}^{+}$ and $\tilde{\mathcal{A}}^{\times}$ are the new
modified polarization tensors due to the presence of the lens. They are
defined as

$\tilde{\mathcal{A}}^{+}_{ij}=\mathcal{A}^{+}_{ij}-\frac{1}{2}(\vec{\gamma}\cdot\vec{x})\mathcal{A}^{x}_{ij}+\frac{1}{2}(\vec{\gamma}\cdot\vec{y})\mathcal{A}^{y}_{ij}$
(78)

and

$\tilde{\mathcal{A}}^{\times}_{ij}=\mathcal{A}^{+}_{ij}-\frac{1}{2}(\vec{\gamma}\cdot\vec{y})\mathcal{A}^{x}_{ij}-\frac{1}{2}(\vec{\gamma}\cdot\vec{x})\mathcal{A}^{y}_{ij}$
(79)

where $\mathcal{A}^{x}_{ij}=x_{i}k_{j}+k_{i}x_{j}$ and
$\mathcal{A}^{y}_{ij}=y_{i}k_{j}+k_{i}y_{j}$ are the polarization tensors for
the vector -x and vector -y polarizations, respectively. One can see that the
effect of the lens is that it introduces additional polarizations in GR, which
is an illusion (see e.g.Hou:2019wdg ). Namely, one can always choose
appropriate coordinates to remove these non-tensorial polarizations.
Therefore, we do not deal here with any additional physical effect.

###  V.3 Palatini $f(\hat{R})$ gravity

Similarly, to evaluate the perturbed quantities we can use the perturbed
geodesic equation similar to the perturbation of light rays (fleurry_thesis ).
Substituting the relation $k^{\mu}=\bar{k}^{\mu}+\delta k^{\mu}$ into the
geodesic equation and evaluating it to first order in perturbation we obtain
the following relation for the tetrad $\delta k^{\mu}$ (and similarly for
other tetrads),

$\bar{k}^{\nu}\hat{\nabla}_{\nu}\delta
k^{\mu}=-\delta\hat{\Gamma}^{\mu}_{\nu\rho}\bar{k}^{\nu}\bar{k}^{\rho}\,.$
(80)

To evaluate the above equation, one must calculate the non-vanishing
components of the affine connection corresponding to the background geometry
associated with the lensing mass (see Appendix B). Since we are interested in
a static case, the conformal factor $f^{\prime}(\hat{R})=f^{\prime}(T)$, where
$T$ is the trace of the energy-momentum tensor, is independent of time, and
hence we will deal with the spatial derivatives only in the above expression.
Analogously as in GR, we deal with the following solutions of (80):

$\delta
x^{\mu}=\frac{1}{2}\Big{(}\vec{\gamma}^{P}\cdot\vec{x},-(\vec{\gamma}^{P}\cdot\vec{x})\vec{k}\Big{)}$
(81) $\delta
y^{\mu}=\frac{1}{2}\Big{(}\vec{\gamma}^{P}\cdot\vec{y},-(\vec{\gamma}^{P}\cdot\vec{y})\vec{k}\Big{)},$
(82)

where

$\vec{\gamma}^{P}=\nabla_{\vec{x}}(f^{\prime}\Psi(\vec{x}))$ (83)

is the gravitational deflection due to the lens in this theory, where as usual
$\Psi(x)$ is the lens potential in a given theory of gravity and the
deflection angle $\gamma^{P}$ reduces to GR when $f^{\prime}=1$ ($\alpha=0$
for the considered cases). Recall that $\nabla_{\perp}$ is the projected
gradient on the plane perpendicular to the GW propagation. Therefore, the wave
tensor (49) can be rewritten by changing the wave vector as
$x^{\mu}\rightarrow x^{\mu}+\delta x^{\mu}$, $y^{\mu}\rightarrow
y^{\mu}+\delta y^{\mu}$. The polarization tensor undergoes a rotation and the
wave tensor becomes

$\xi^{(0)}_{\mu\nu}=H_{+}\tilde{\mathcal{A}}^{+}_{\mu\nu}+H_{\times}\tilde{\mathcal{A}}^{\times}_{\mu\nu}$
(84)

where $\tilde{\mathcal{A}}^{+}$ and $\tilde{\mathcal{A}}^{\times}$ are the new
modified polarization tensors defined as

$\tilde{\mathcal{A}}^{+}_{ij}=\mathcal{A}^{+}_{ij}-\frac{1}{2}(\vec{\gamma}^{P}\cdot\vec{x})\mathcal{A}^{x}_{ij}+\frac{1}{2}(\vec{\gamma}^{P}\cdot\vec{y})\mathcal{A}^{y}_{ij}$
(85)

and

$\tilde{\mathcal{A}}^{\times}_{ij}=\mathcal{A}^{+}_{ij}-\frac{1}{2}(\vec{\gamma}^{P}\cdot\vec{y})\mathcal{A}^{x}_{ij}-\frac{1}{2}(\vec{\gamma}^{P}\cdot\vec{x})\mathcal{A}^{y}_{ij}$
(86)

where $\mathcal{A}^{x}_{ij}=x_{i}k_{j}+k_{i}x_{j}$ and
$\mathcal{A}^{y}_{ij}=y_{i}k_{j}+k_{i}y_{j}$ are the polarization tensors for
the vector -x and vector- y polarizations respectively.

It is evident from the above expressions (85) and (86) that the polarization
plane is rotated in the geometric optics limit in the similar manner to what
we have observed in GR Hou:2019wdg . However, the rotation is different from
GR as a consequence of the extra conformal factor $f^{\prime}(\hat{R})$
appearing in $\vec{\gamma}^{P}$. Let us recall that it is customary in the
literature to neglect the spin 2 nature of gravitational waves and it is
acceptable as the polarization tensor is parallel transported along the
geodesic as we saw in section III.

###  V.4 Lensing by Singular Isothermal sphere

Singular Isothermal Sphere (SIS) Kormann1994  is one of the simplest
parametrizations of the axially symmetric lens model where the matter content
of the lens behaves like an ideal gas in thermal and hydrostatic equilibrium.
The SIS density profile is given by the density distribution

$\rho(r)=\frac{\sigma_{v}^{2}}{2\pi Gr^{2}}$ (87)

where $\sigma_{v}$ is the velocity distribution of the gas particles. Using
the above density distribution and the Poisson equation given in (63) it is
easy to calculate the convergence

$\displaystyle\kappa({\vec{\xi}})$ $\displaystyle=$
$\displaystyle\frac{1}{c^{2}}\frac{D_{L}D_{LS}}{D_{S}}\int_{-\infty}^{\infty}\nabla^{2}\Phi(\vec{\xi},z)dz$
(88) $\displaystyle=$ $\displaystyle\frac{4\pi
G}{c^{2}}\frac{D_{L}D_{LS}}{D_{S}}\int_{-\infty}^{\infty}(\rho+2\alpha\nabla^{2}\rho)dz$

which provides

$\kappa({\vec{\xi}})=\frac{2\pi\sigma_{v}^{2}}{c^{2}}\frac{D_{L}D_{LS}}{D_{S}}\Big{[}\frac{1}{\xi}+\frac{2\alpha}{\xi^{3}}\Big{]}.$
(89)

The above equation can be re-written in terms of the dimensionless variable
$x=\xi/\xi_{0}$,

$\kappa(x)=\frac{1}{2x}+\frac{\alpha}{x^{3}\xi_{0}^{2}}\,,$ (90)

where the length scale $\xi_{0}$ and is given by

$\xi_{0}=4\pi\frac{\sigma_{v}^{2}}{c^{2}}\frac{D_{L}D_{LS}}{D_{S}}.$ (91)

As expected, taking $\alpha=0$ recovers the convergence in GR for the SIS
density profile.

Lensing potential $\Psi(\vec{x})=\Psi(x)$ can be obtained from the 2D Poisson
equation $\nabla^{2}_{\vec{x}}\Psi(x)=2\kappa(x).$ For this purpose it is
convenient to introduce polar coordinates in the lens plane:
$\vec{x}=x(\cos{\varphi},\sin{\varphi})$. Then
$\nabla^{2}_{\vec{x}}=\frac{1}{x}\frac{\partial}{\partial
x}\left(x\frac{\partial}{\partial
x}\right)+\frac{1}{x^{2}}\frac{\partial^{2}}{\partial\varphi^{2}}$ and 2D
Poisson equation reads

$\frac{1}{x}\frac{d}{dx}\left(x\frac{d\Psi}{dx}\right)=\frac{1}{x}+\frac{2\alpha}{x^{3}\xi_{0}^{2}}.$
(92)

Its solution is

$\Psi(x)=x+\frac{2\alpha}{\xi_{0}^{2}x}\,.$ (93)

Taking the derivative of the lensing potential (83) gives us the corresponding
deflection angle $\vec{\gamma}$ for the Palatini polynomial $f(\hat{R})$
models,

$\vec{\gamma}^{P}(x)=\frac{\vec{x}}{x}-\frac{2\alpha}{\xi_{0}^{2}}\frac{\vec{x}}{x^{3}}\,.$
(94)

In addition to the convergence $\kappa(x)$ being informative about Ricci part
of curvature effect on the bundle of null geodesics, shear modulus999For more
details see Appendix C. $\gamma(x)$ quantifies the magnitude of the Weyl part
of the curvature tensor. In the case of Palatini theory one has, for the SIS
lens:

$\displaystyle\gamma=\frac{1}{2x^{3}}\left|\frac{6\alpha}{\xi_{0}^{2}}-x^{2}\right|$
(95)

and the magnification is

$\displaystyle\mu=x^{6}\left(x^{3}-\frac{4\alpha}{\xi_{0}^{2}}\right)^{-1}\left(x^{3}-x^{2}+\frac{2\alpha}{\xi_{0}^{2}}\right)^{-1}\,,$
(96)

where again GR formulae are recovered in the limit of vanishing $\alpha$.

##  VI Conclusions

The aim of the present research was to examine the effects of Palatini
$f(\hat{R})$ gravity on the propagation of gravitational waves and
gravitational lensing of them, in particular. To our knowledge, this is the
first work in which the effects of metric-affine gravity were studied in the
context of the lensing of gravitational waves. Firstly, we have confirmed the
results obtained by previous works in the literature jimenez2015tensor ;
jimenez2018born , demonstrating that gravitational waves propagate along
autoparallel curves in this theory which coincides with geodesics in general
relativity. This consistency reinforces the validity of our approach and
contributes to the overall understanding of gravitational waves in modified
theories of gravity. Following the standard procedure, we have studied using
WKB approximation the geometric optics limit in Palatini $f(\hat{R})$ gravity.
This allowed us to successfully compute the evolution of gravitational wave
amplitude and explore the effects of metric-affine gravity on the rotation of
the gravitational wave polarization plane. Furthermore, as a useful and
realistic example we discussed the SIS model, which is a robust model for a
galaxy acting as a lens, in the context of Palatini $f(\hat{R})$ gravity.

The examined model of gravity reduces to GR in vacuum, it does not provide any
changes in the gravitational wave amplitude and other lensing properties for a
point lens. Therefore, one has to consider the propagation of GWs through the
region of space-time where a non-zero energy density distribution is present
or in a close vicinity of it. Note that the non-relativistic gravitational
potential describing such a situation carries an additional term due to the
modification of the Poisson equation (63). An example of the first case is
when a GW signal emitted from the source travels not exactly through the
vacuum but rather through the dark matter halo of the host galaxy, through the
dark matter of the galaxy cluster to which the host belongs, and through the
dark matter halo of our Galaxy. This attenuates the GW amplitude in the
metric-affine gravity theory studied. Moreover, we have shown that the GW
amplitude evolution equation can be expressed as a correction to the measured
luminosity distance (62) of the source. One should note that the corrections
depend on the matter density distribution along the line of sight through the
conformal factor. As an example, we considered a galaxy halo modelled by NFW
(SIS) density profile. Taking the value of $\alpha\sim 10^{9}$ m2 we found
that the corrections are of the order of $10^{-42}\,(10^{-43})$, which is
negligible. However, the corrections could be significant – of order of
$10^{-2}$ if one uses the constraints for $\alpha\sim 10^{49}$ m2 obtained
from SNIa and BAO observations Gomes:2023xzk . Furthermore, we have derived
the corrections to the polarization tensor (85) and (86) in theory considered
and found that the polarization plane is rotated and a vector polarization
appears. However, just like in GR this is illusory and can be removed by
proper gauge choice. It is worth mentioning that these corrections depend on
the modified deflection angle (94) and hence they are theory dependent. As an
example of the second case, we studied the SIS lens model and we demonstrated
that Palatini $f(\hat{R})$ gravity indeed modifies the convergence (90),
lensing potential (93), and the deflection angle (94), reducing to the GR
limit when $\alpha$ vanishes. Let us stress that the exemplary value of
$\alpha$ used by us is not the firmly established one and can be updated in
the future. Hence the claims of negligible modified gravity effects are not
firmly established, as well. Therefore, the results we obtained could be
valuable for further metric-affine gravity vs. GR tests involving lensing of
GWs and comparison of luminosity distances measured from electromagnetic and
GW sources.

### Acknowledgements

This work was supported by the European Regional Development Fund through the
Center of Excellence TK133 “The Dark Side of the Universe.” SH acknowledges
the Dora Plus grant for visiting students supported by European Regional
Development Fund and Laboratory of Theoretical Physics, University of Tartu
where this research is carried out. AW acknowledges financial support from
MICINN (Spain) Ayuda Juan de la Cierva - incorporación 2020 No.
IJC2020-044751-I. LJ and MS were supported by the Estonian Research Council
grant PRG356 “Gauge Gravity”.

## References

  * (1) LIGO Scientific, Virgo Collaboration, B. P. Abbott et al., “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116 (2016) no. 6, 061102, arXiv:1602.03837 [gr-qc].
  * (2) LIGO Scientific, VIRGO, KAGRA Collaboration, R. Abbott et al., “GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run,” arXiv:2111.03606 [gr-qc].
  * (3) J. Antoniadis et al., “The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals,” arXiv:2306.16214 [astro-ph.HE].
  * (4) NANOGrav Collaboration, G. Agazie et al., “The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background,” Astrophys. J. Lett. 951 (2023) no. 1, L8, arXiv:2306.16213 [astro-ph.HE].
  * (5) H. Xu et al., “Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I,” Res. Astron. Astrophys. 23 (2023) no. 7, 075024, arXiv:2306.16216 [astro-ph.HE].
  * (6) D. J. Reardon et al., “Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array,” Astrophys. J. Lett. 951 (2023) no. 1, L6, arXiv:2306.16215 [astro-ph.HE].
  * (7) J. Sakstein and B. Jain, “Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories,” Phys. Rev. Lett. 119 (2017) no. 25, 251303, arXiv:1710.05893 [astro-ph.CO].
  * (8) J. M. Ezquiaga and M. Zumalacárregui, “Dark Energy After GW170817: Dead Ends and the Road Ahead,” Phys. Rev. Lett. 119 (2017) no. 25, 251304, arXiv:1710.05901 [astro-ph.CO].
  * (9) T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller, and I. Sawicki, “Strong constraints on cosmological gravity from GW170817 and GRB 170817A,” Phys. Rev. Lett. 119 (2017) no. 25, 251301, arXiv:1710.06394 [astro-ph.CO].
  * (10) P. Schneider, J. Ehlers, and E. E. Falco, Gravitational Lenses.  Astronomy and Astrophysics Library. Springer, 1992.
  * (11) O. A. Hannuksela, K. Haris, K. K. Y. Ng, S. Kumar, A. K. Mehta, D. Keitel, T. G. F. Li, and P. Ajith, “Search for gravitational lensing signatures in LIGO-Virgo binary black hole events,” Astrophys. J. Lett. 874 (2019) no. 1, L2, arXiv:1901.02674 [gr-qc].
  * (12) LIGO Scientific, VIRGO Collaboration, R. Abbott et al., “Search for Lensing Signatures in the Gravitational-Wave Observations from the First Half of LIGO–Virgo’s Third Observing Run,” Astrophys. J. 923 (2021) no. 1, 14, arXiv:2105.06384 [gr-qc].
  * (13) LIGO Scientific, VIRGO, KAGRA Collaboration, R. Abbott et al., “Search for gravitational-lensing signatures in the full third observing run of the LIGO-Virgo network,” arXiv:2304.08393 [gr-qc].
  * (14) M. Sereno, P. Jetzer, A. Sesana, and M. Volonteri, “Cosmography with strong lensing of LISA gravitational wave sources,” Mon. Not. Roy. Astron. Soc. 415 (2011) 2773, arXiv:1104.1977 [astro-ph.CO].
  * (15) S. Cao, J. Qi, Z. Cao, M. Biesiada, J. Li, Y. Pan, and Z.-H. Zhu, “Direct test of the FLRW metric from strongly lensed gravitational wave observations,” Sci. Rep. 9 (2019) no. 1, 11608, arXiv:1910.10365 [astro-ph.CO].
  * (16) K. Liao, X.-L. Fan, X.-H. Ding, M. Biesiada, and Z.-H. Zhu, “Precision cosmology from future lensed gravitational wave and electromagnetic signals,” Nature Commun. 8 (2017) no. 1, 1148, arXiv:1703.04151 [astro-ph.CO]. [Erratum: Nature Commun. 8, 2136 (2017)].
  * (17) Y. Li, X. Fan, and L. Gou, “Constraining Cosmological Parameters in the FLRW Metric with Lensed GW+EM Signals,” Astrophys. J. 873 (2019) no. 1, 37, arXiv:1901.10638 [astro-ph.CO].
  * (18) M. Biesiada and S. Harikumar, “Gravitational Lensing of Continuous Gravitational Waves,” Universe 7 (2021) no. 12, 502, arXiv:2111.05963 [gr-qc].
  * (19) M. Grespan and M. Biesiada, “Strong Gravitational Lensing of Gravitational Waves: A Review,” Universe 9 (2023) no. 5, 200.
  * (20) K.-H. Lai, O. A. Hannuksela, A. Herrera-Martín, J. M. Diego, T. Broadhurst, and T. G. F. Li, “Discovering intermediate-mass black hole lenses through gravitational wave lensing,” Phys. Rev. D 98 (2018) no. 8, 083005, arXiv:1801.07840 [gr-qc].
  * (21) A. K. Meena, “Gravitational lensing of gravitational waves: Probing intermediate mass black holes in galaxy lenses with global minima,” arXiv:2305.02880 [astro-ph.CO].
  * (22) J. M. Diego, “Constraining the abundance of primordial black holes with gravitational lensing of gravitational waves at LIGO frequencies,” Phys. Rev. D 101 (2020) no. 12, 123512, arXiv:1911.05736 [astro-ph.CO].
  * (23) M. Oguri and R. Takahashi, “Probing Dark Low-mass Halos and Primordial Black Holes with Frequency-dependent Gravitational Lensing Dispersions of Gravitational Waves,” Astrophys. J. 901 (2020) no. 1, 58, arXiv:2007.01936 [astro-ph.CO].
  * (24) J. M. Ezquiaga and M. Zumalacárregui, “Gravitational wave lensing beyond general relativity: birefringence, echoes and shadows,” Phys. Rev. D 102 (2020) no. 12, 124048, arXiv:2009.12187 [gr-qc].
  * (25) S. Goyal, K. Haris, A. K. Mehta, and P. Ajith, “Testing the nature of gravitational-wave polarizations using strongly lensed signals,” Phys. Rev. D 103 (2021) no. 2, 024038, arXiv:2008.07060 [gr-qc].
  * (26) X.-L. Fan, K. Liao, M. Biesiada, A. Piorkowska-Kurpas, and Z.-H. Zhu, “Speed of Gravitational Waves from Strongly Lensed Gravitational Waves and Electromagnetic Signals,” Phys. Rev. Lett. 118 (2017) no. 9, 091102, arXiv:1612.04095 [gr-qc].
  * (27) T. E. Collett and D. Bacon, “Testing the speed of gravitational waves over cosmological distances with strong gravitational lensing,” Phys. Rev. Lett. 118 (2017) no. 9, 091101, arXiv:1602.05882 [astro-ph.HE].
  * (28) V. K. Sharma, S. Harikumar, M. Grespan, M. Biesiada, and M. M. Verma, “Probing massive gravitons in $f(R)$ with lensed gravitational waves,” arXiv:2310.10346 [gr-qc].
  * (29) A. Mishra, N. V. Krishnendu, and A. Ganguly, “Unveiling Microlensing Biases in Testing General Relativity with Gravitational Waves,” arXiv:2311.08446 [gr-qc].
  * (30) T. Baker and M. Trodden, “Multimessenger time delays from lensed gravitational waves,” Phys. Rev. D 95 (2017) no. 6, 063512, arXiv:1612.02004 [astro-ph.CO].
  * (31) M. Maggiore et al., “Science Case for the Einstein Telescope,” JCAP 03 (2020) 050, arXiv:1912.02622 [astro-ph.CO].
  * (32) M. Evans et al., “A Horizon Study for Cosmic Explorer: Science, Observatories, and Community,” arXiv:2109.09882 [astro-ph.IM].
  * (33) S. Kawamura et al., “Current status of space gravitational wave antenna DECIGO and B-DECIGO,” PTEP 2021 (2021) no. 5, 05A105, arXiv:2006.13545 [gr-qc].
  * (34) P. A.-S. et al., “Laser interferometer space antenna,” arXiv:1702.00786 [astro-ph.IM].
  * (35) S.-S. Li, S. Mao, Y. Zhao, and Y. Lu, “Gravitational lensing of gravitational waves: A statistical perspective,” Mon. Not. Roy. Astron. Soc. 476 (2018) no. 2, 2220–2229, arXiv:1802.05089 [astro-ph.CO].
  * (36) M. Biesiada, X. Ding, A. Piorkowska, and Z.-H. Zhu, “Strong gravitational lensing of gravitational waves from double compact binaries - perspectives for the Einstein Telescope,” JCAP 10 (2014) 080, arXiv:1409.8360 [astro-ph.HE].
  * (37) L. Yang, X. Ding, M. Biesiada, K. Liao, and Z.-H. Zhu, “How Does the Earth’s Rotation Affect Predictions of Gravitational Wave Strong Lensing Rates?,” Astrophys. J. 874 (2019) no. 2, 139, arXiv:1903.11079 [astro-ph.GA].
  * (38) L. Yang, S. Wu, K. Liao, X. Ding, Z. You, Z. Cao, M. Biesiada, and Z.-H. Zhu, “Event rate predictions of strongly lensed gravitational waves with detector networks and more realistic templates,” Mon. Not. Roy. Astron. Soc. 509 (2021) no. 3, 3772–3778, arXiv:2105.07011 [astro-ph.GA].
  * (39) A. Piórkowska-Kurpas, S. Hou, M. Biesiada, X. Ding, S. Cao, X. Fan, S. Kawamura, and Z.-H. Zhu, “Inspiraling Double Compact Object Detection and Lensing Rate: Forecast for DECIGO and B-DECIGO,” Astrophys. J. 908 (2021) no. 2, 196, arXiv:2005.08727 [astro-ph.HE].
  * (40) T. T. Nakamura and S. Deguchi, “Wave Optics in Gravitational Lensing,” Prog. Theor. Phys. Suppl. 133 (1999) 137–153.
  * (41) R. Takahashi and T. Nakamura, “Wave effects in gravitational lensing of gravitational waves from chirping binaries,” Astrophys. J. 595 (2003) 1039–1051, arXiv:astro-ph/0305055.
  * (42) C. M. Will, Theory and Experiment in Gravitational Physics.  Cambridge University Press, 1993.
  * (43) D. M. Eardley, D. L. Lee, A. P. Lightman, R. V. Wagoner, and C. M. Will, “Gravitational-wave observations as a tool for testing relativistic gravity,” Phys. Rev. Lett. 30 (1973) 884–886.
  * (44) K. S. Thorne, D. L. Lee, and A. P. Lightman, “Foundations for a Theory of Gravitation Theories,” Phys. Rev. D 7 (1973) 3563–3578.
  * (45) S. Hou, X.-L. Fan, and Z.-H. Zhu, “Gravitational Lensing of Gravitational Waves: Rotation of Polarization Plane,” Phys. Rev. D 100 (2019) no. 6, 064028, arXiv:1907.07486 [gr-qc].
  * (46) C. Dalang, G. Cusin, and M. Lagos, “Polarization distortions of lensed gravitational waves,” Phys. Rev. D 105 (2022) no. 2, 024005, arXiv:2104.10119 [gr-qc].
  * (47) G. Cusin and M. Lagos, “Gravitational wave propagation beyond geometric optics,” Phys. Rev. D 101 (2020) no. 4, 044041, arXiv:1910.13326 [gr-qc].
  * (48) G. J. Olmo, “Palatini Approach to Modified Gravity: f(R) Theories and Beyond,” Int. J. Mod. Phys. D 20 (2011) 413–462, arXiv:1101.3864 [gr-qc].
  * (49) T. P. Sotiriou, “Curvature scalar instability in f(R) gravity,” Phys. Lett. B 645 (2007) 389–392, arXiv:gr-qc/0611107.
  * (50) A. Stachowski, M. Szydłowski, and A. Borowiec, “Starobinsky cosmological model in Palatini formalism,” Eur. Phys. J. C 77 (2017) no. 6, 406, arXiv:1608.03196 [gr-qc].
  * (51) M. Szydłowski, A. Stachowski, and A. Borowiec, “Emergence of running dark energy from polynomial f(R) theory in Palatini formalism,” Eur. Phys. J. C 77 (2017) no. 9, 603, arXiv:1707.01948 [gr-qc].
  * (52) F. A. Teppa Pannia, S. E. Perez Bergliaffa, and N. Manske, “Cosmography and the redshift drift in Palatini $f({\cal R})$ theories,” Eur. Phys. J. C 79 (2019) no. 3, 267, arXiv:1811.08176 [gr-qc].
  * (53) P. Pinto, L. Del Vecchio, L. Fatibene, and M. Ferraris, “Extended cosmology in Palatini f(R)-theories,” JCAP 11 (2018) 044, arXiv:1807.00397 [gr-qc].
  * (54) S. Camera, S. Capozziello, L. Fatibene, and A. Orizzonte, “The effective equation of state in Palatini $f({{\mathcal{R}}})$ cosmology,” Eur. Phys. J. Plus 138 (2023) no. 2, 180, arXiv:2212.13825 [gr-qc].
  * (55) I. D. Gialamas, A. Karam, T. D. Pappas, and E. Tomberg, “Implications of Palatini gravity for inflation and beyond,” arXiv:2303.14148 [gr-qc].
  * (56) G. J. Olmo, D. Rubiera-Garcia, and A. Wojnar, “Minimum main sequence mass in quadratic Palatini $f(R)$ gravity,” Phys. Rev. D 100 (2019) no. 4, 044020, arXiv:1906.04629 [gr-qc].
  * (57) G. J. Olmo, D. Rubiera-Garcia, and A. Wojnar, “Stellar structure models in modified theories of gravity: Lessons and challenges,” Phys. Rept. 876 (2020) 1–75, arXiv:1912.05202 [gr-qc].
  * (58) A. Wojnar, “Early evolutionary tracks of low-mass stellar objects in modified gravity,” Phys. Rev. D 102 (2020) no. 12, 124045, arXiv:2007.13451 [gr-qc].
  * (59) A. Wojnar, “Lithium abundance is a gravitational model dependent quantity,” Phys. Rev. D 103 (2021) no. 4, 044037, arXiv:2009.10983 [gr-qc].
  * (60) M. Benito and A. Wojnar, “Cooling process of brown dwarfs in Palatini f(R) gravity,” Phys. Rev. D 103 (2021) no. 6, 064032, arXiv:2101.02146 [gr-qc].
  * (61) A. Wojnar, “Jupiter and jovian exoplanets in Palatini f(R¯) gravity,” Phys. Rev. D 104 (2021) no. 10, 104058, arXiv:2108.13528 [gr-qc].
  * (62) A. Kozak and A. Wojnar, “Metric-affine gravity effects on terrestrial exoplanet profiles,” Phys. Rev. D 104 (2021) no. 8, 084097, arXiv:2106.14219 [gr-qc].
  * (63) L. Sarmah, S. Kalita, and A. Wojnar, “Stability criterion for white dwarfs in Palatini f(R) gravity,” Phys. Rev. D 105 (2022) no. 2, 024028, arXiv:2111.08029 [gr-qc].
  * (64) S. Kalita, L. Sarmah, and A. Wojnar, “Metric-affine effects in crystallization processes of white dwarfs,” Phys. Rev. D 107 (2023) no. 4, 044072, arXiv:2212.04918 [gr-qc].
  * (65) A. Kozak and A. Wojnar, “Planetary seismology as a test of modified gravity proposals,” Phys. Rev. D 108 (2023) no. 4, 044055, arXiv:2303.17213 [gr-qc].
  * (66) J. D. Toniato, D. C. Rodrigues, and A. Wojnar, “Palatini $f(R)$ gravity in the solar system: post-Newtonian equations of motion and complete PPN parameters,” Phys. Rev. D 101 (2020) no. 6, 064050, arXiv:1912.12234 [gr-qc].
  * (67) A. Bonino, S. Camera, L. Fatibene, and A. Orizzonte, “Solar System tests in Brans–Dicke and Palatini $f({\mathcal{R}})$-theories,” Eur. Phys. J. Plus 135 (2020) no. 12, 951, arXiv:2011.06303 [gr-qc].
  * (68) A. Hernandez-Arboleda, D. C. Rodrigues, J. D. Toniato, and A. Wojnar, “Palatini $f(R)$ gravity tests in the weak field limit: Solar System, seismology and galaxies,” arXiv:2306.04475 [gr-qc].
  * (69) M. E. S. Alves, O. D. Miranda, and J. C. N. de Araujo, “Probing the f(R) formalism through gravitational wave polarizations,” Phys. Lett. B 679 (2009) 401–406, arXiv:0908.0861 [gr-qc].
  * (70) D. Iosifidis, “The Perfect Hyperfluid of Metric-Affine Gravity: The Foundation,” JCAP 04 (2021) 072, arXiv:2101.07289 [gr-qc].
  * (71) T. Koivisto, “Covariant conservation of energy momentum in modified gravities,” Class. Quant. Grav. 23 (2006) 4289–4296, arXiv:gr-qc/0505128.
  * (72) M. Ferraris, M. Francaviglia, and I. Volovich, “Universal gravitational equations,” Nuovo Cim. B 108 (1993) 1313–1317.
  * (73) M. Ferraris, M. Francaviglia, and I. Volovich, “The Universality of vacuum Einstein equations with cosmological constant,” Class. Quant. Grav. 11 (1994) 1505–1517, arXiv:gr-qc/9303007.
  * (74) P. K. Schwartz and D. Giulini, “Post-Newtonian Hamiltonian description of an atom in a weak gravitational field,” Phys. Rev. A 100 (2019) no. 5, 052116, arXiv:1908.06929 [quant-ph].
  * (75) G. J. Olmo, “Hydrogen atom in Palatini theories of gravity,” Phys. Rev. D 77 (2008) 084021, arXiv:0802.4038 [gr-qc].
  * (76) G. J. Olmo, “Violation of the Equivalence Principle in Modified Theories of Gravity,” Phys. Rev. Lett. 98 (2007) 061101, arXiv:gr-qc/0612002.
  * (77) A. Wojnar, “Fermi gas and modified gravity,” Phys. Rev. D 107 (2023) no. 4, 044025, arXiv:2208.04023 [gr-qc].
  * (78) J. Naf, P. Jetzer, and M. Sereno, “On Gravitational Waves in Spacetimes with a Nonvanishing Cosmological Constant,” Phys. Rev. D 79 (2009) 024014, arXiv:0810.5426 [astro-ph].
  * (79) G. J. Olmo, “The Gravity Lagrangian according to solar system experiments,” Phys. Rev. Lett. 95 (2005) 261102, arXiv:gr-qc/0505101.
  * (80) A. Kozak and A. Wojnar, “Earthquakes as probing tools for gravity theories,” arXiv:2308.01784 [gr-qc].
  * (81) E. Lope-Oter and A. Wojnar, “Constraining Palatini gravity with GR-independent equations of state,” arXiv:2306.00870 [gr-qc].
  * (82) A. Wojnar, “Unveiling Phase Space Modifications: A Clash of Modified Gravity and the Generalized Uncertainty Principle,” arXiv:2311.14066 [gr-qc].
  * (83) D. C. Rodrigues, A. Hernandez-Arboleda, and A. Wojnar, “Normalized additional velocity distribution: Testing the radial profile of dark matter halos and MOND,” Phys. Dark Univ. 41 (2023) 101230, arXiv:2204.03762 [astro-ph.GA].
  * (84) D. A. Gomes, R. Briffa, A. Kozak, J. Levi Said, M. Saal, and A. Wojnar, “Cosmological constraints of Palatini $f(\mathcal{R})$ gravity,” arXiv:2310.17339 [gr-qc].
  * (85) J. Beltran Jimenez, L. Heisenberg, and G. J. Olmo, “Tensor perturbations in a general class of Palatini theories,” JCAP 06 (2015) 026, arXiv:1504.00295 [gr-qc].
  * (86) J. Beltran Jimenez, L. Heisenberg, G. J. Olmo, and D. Rubiera-Garcia, “On gravitational waves in Born-Infeld inspired non-singular cosmologies,” JCAP 10 (2017) 029, arXiv:1707.08953 [hep-th]. [Erratum: JCAP 08, E01 (2018)].
  * (87) J. Beltran Jimenez, L. Heisenberg, G. J. Olmo, and D. Rubiera-Garcia, “Born–Infeld inspired modifications of gravity,” Phys. Rept. 727 (2018) 1–129, arXiv:1704.03351 [gr-qc].
  * (88) M. Maggiore, Gravitational Waves. Vol. 1: Theory and Experiments.  Oxford University Press, 2007.
  * (89) A. Garoffolo, G. Tasinato, C. Carbone, D. Bertacca, and S. Matarrese, “Gravitational waves and geometrical optics in scalar-tensor theories,” JCAP 11 (2020) 040, arXiv:1912.08093 [gr-qc].
  * (90) P. Fleury, “Light propagation in inhomogeneous and anisotropic cosmologies,” arXiv:1511.03702 [gr-qc].
  * (91) C. Deffayet and K. Menou, “Probing Gravity with Spacetime Sirens,” Astrophys. J. Lett. 668 (2007) L143–L146, arXiv:0709.0003 [astro-ph].
  * (92) I. D. Saltas, I. Sawicki, L. Amendola, and M. Kunz, “Anisotropic Stress as a Signature of Nonstandard Propagation of Gravitational Waves,” Phys. Rev. Lett. 113 (2014) no. 19, 191101, arXiv:1406.7139 [astro-ph.CO].
  * (93) A. Nishizawa, “Generalized framework for testing gravity with gravitational-wave propagation. I. Formulation,” Phys. Rev. D 97 (2018) no. 10, 104037, arXiv:1710.04825 [gr-qc].
  * (94) E. Belgacem, Y. Dirian, S. Foffa, and M. Maggiore, “Gravitational-wave luminosity distance in modified gravity theories,” Phys. Rev. D 97 (2018) no. 10, 104066, arXiv:1712.08108 [astro-ph.CO].
  * (95) E. Newman and R. Penrose, “An Approach to gravitational radiation by a method of spin coefficients,” J. Math. Phys. 3 (1962) 566–578.
  * (96) S. Chandrasekhar, The mathematical theory of black holes.  Clarendon Press, 1985.
  * (97) P. Fleury, C. Pitrou, and J.-P. Uzan, “Light propagation in a homogeneous and anisotropic universe,” Phys. Rev. D 91 (2015) no. 4, 043511, arXiv:1410.8473 [gr-qc].
  * (98) J. F. Navarro, C. S. Frenk, and S. D. M. White, “A Universal density profile from hierarchical clustering,” Astrophys. J. 490 (1997) 493–508, arXiv:astro-ph/9611107.
  * (99) S. Ettori, V. Ghirardini, D. Eckert, E. Pointecouteau, F. Gastaldello, M. Sereno, M. Gaspari, S. Ghizzardi, M. Roncarelli, and M. Rossetti, “Hydrostatic mass profiles in X-COP galaxy clusters,” Astron. Astrophys. 621 (2019) A39, arXiv:1805.00035 [astro-ph.CO].
  * (100) M. Meneghetti, Introduction to Gravitational Lensing; With Python Examples, vol. 956.  2021\.
  * (101) A. O. Petters, H. Levine, and J. Wambsganss, Singularity theory and gravitational lensing.  2001\.
  * (102) T. Treu, “Strong lensing by galaxies,” Annual Review of Astronomy and Astrophysics 48 (2010) no. 1, 87–125. https://doi.org/10.1146/annurev-astro-081309-130924.
  * (103) R. Kormann, P. Schneider, and M. Bartelmann, “Isothermal elliptical gravitational lens models,” Astronomy and Astrophysics 284 (1994) 285–299.

##  Appendix A Sachs formalism

With the notations of Sec.(IV.2), let us introduce an orthonormal basis (the
Sachs basis) $(s^{\mu}_{A})$, $A=\\{1,2\\}$, satisfying the conditions

$s^{\mu}_{A}u_{\mu}=s^{\mu}_{A}d_{\mu}=0,\qquad
s^{\mu}_{A}s_{B\mu}=\delta_{AB}.$ (97)

In the further part, we follow fleurry_thesis .

One of the ways to study the beam pattern is provided by the solutions of the
deviation equation (51) provided in the Sachs basis

$\frac{d^{2}\chi^{A}}{d\lambda^{2}}=\mathcal{R}^{A}_{B}\chi^{B}$ (98)

where $\mathcal{R}^{A}_{B}$ is the optical tidal matrix which is symmetric
under the exchange of their indices. Moreover, it can be decomposed into pure
trace and trace-free parts as

$\mathcal{R}=\begin{pmatrix}\mathscr{R}&0\\\
0&\mathscr{R}\end{pmatrix}+\begin{pmatrix}-\text{Re}\mathscr{W}&\text{Im}\mathscr{W}\\\
\text{Im}\mathscr{W}&\text{Re}\mathscr{W}\end{pmatrix}$ (99)

where $\mathscr{R}=-\frac{1}{2}R_{\mu\nu}k^{\mu}k^{\nu}$ is the Ricci lensing
term, responsible for an homothetic transformation of the beam pattern. On the
other hand,

$\mathscr{W}=-\frac{1}{2}C_{\mu\rho\nu\sigma}(s^{\mu}_{1}-is^{\mu}_{2})k^{\rho}k^{\nu}(s^{\sigma}_{1}-is^{\sigma}_{2}),$
(100)

is the Weyl lensing matrix responsible for the elongation and contraction of
the beam pattern in different directions. The Weyl tensor
$C_{\mu\rho\nu\sigma}$ has the same symmetries as the Riemann tensor and it is
trace free. As we will see, these quantities are useful in lensing studies to
express the wave’s polarization amplitude.

Let us now denote the optical tidal matrix as $\mathcal{R}$. Then, the Sachs
equation (98) can be written in terms of the Jacobi matrix $\mathcal{D}$ which
relates position of images on the observers’ celestial sphere to the physical
separation of the sources:

$\ddot{\mathcal{D}}=\mathcal{R}\mathcal{D},$ (101)

where dot denotes derivative with respect to the affine parameter $\lambda$.
Let us now define the deformation rate of the beam as

$\mathcal{S}=\frac{\dot{\mathcal{D}}}{\mathcal{D}}.$ (102)

Note that the deformation matrix (102) can be expressed in terms of Sach’s
basis as

$\mathcal{S}_{AB}=s^{\mu}_{A}s^{\nu}_{B}\nabla_{\mu}k_{\nu}.$ (103)

Now, we can decompose the matrix $\mathcal{S}$ into anti-symmetric, trace, and
traceless symmetric parts

$\mathcal{S}=\begin{pmatrix}0&\omega\\\
-\omega&0\end{pmatrix}+\begin{pmatrix}\theta&0\\\
0&\theta\end{pmatrix}+\begin{pmatrix}-\sigma_{1}&\sigma_{2}\\\
\sigma_{2}&\sigma_{1}\end{pmatrix},$ (104)

where $\omega$, $\theta$ and $\sigma=\sigma_{1}+i\sigma_{2}$ are optical
scalars defined as

$\omega^{2}=\frac{1}{2}\omega^{\alpha\beta}\omega_{\alpha\beta}\,,\;\;\theta=\frac{1}{2}k^{\alpha}_{\,;\alpha}\,,\;\;\sigma^{2}=\sigma^{\alpha\beta}\sigma_{\alpha\beta}$
(105)

with

$\displaystyle\omega_{\alpha\beta}$ $\displaystyle=k_{[\alpha;\beta]}$ (106)
$\displaystyle\sigma_{\alpha\beta}$
$\displaystyle=k_{(\alpha;\beta)}-\frac{1}{2}Q_{\alpha\beta}\theta\,.$ (107)

The semicolon ;β means the application of $\nabla_{\beta}$ while the round and
rectangular brackets denote symmetrization and antisymetrization in the
indices, respectively. Since we are dealing with a null congruence, the
vorticity tensor $\omega_{\alpha\beta}$ is zero since the null vector was
defined as $k_{\mu}=\Phi,_{\mu}$ after the expression (21). Therefore, the
matrix’s components (104) consist of the expansion $\theta$ and and shear
scalars $\sigma_{1}$ and $\sigma_{2}$ only.

##  Appendix B Weak-field limit equations for Palatini gravity

In the weak-field limit the metric in modified gravity can be written as

$ds^{2}=-(1+2\Phi)dt^{2}+(1-2\Psi)(dx^{2}+dy^{2}+dz^{2}),$ (108)

where $\Phi$ and $\Psi$ are Bardeen potentials. We can directly obtain the
non-vanishing components of the affine connection corresponding to the
background geometry associated with the lensing mass:

$\displaystyle\delta\hat{\Gamma}^{0}_{i0}$
$\displaystyle=-\partial_{i}\Phi+\Phi\partial_{i}f^{\prime}$ (109)
$\displaystyle\delta\hat{\Gamma}^{i}_{00}$ $\displaystyle=\partial_{i}\Psi$
(110) $\displaystyle\delta\hat{\Gamma}^{0}_{ij}$ $\displaystyle=0$ (111)
$\displaystyle\delta\hat{\Gamma}^{i}_{j0}$ $\displaystyle=0$ (112)
$\displaystyle\delta\Gamma^{i}_{jk}$
$\displaystyle=-\partial_{k}\Psi\delta^{i}_{j}-\delta^{i}_{k}\partial_{j}\Psi+\delta_{jk}\partial^{i}\Psi$
$\displaystyle-\partial_{j}\ln f^{\prime}\delta^{i}_{k}\Psi-\partial_{k}\ln
f^{\prime}\delta^{i}_{j}\Psi+\partial^{i}\ln f^{\prime}\delta_{jk}\Psi.$ (113)

##  Appendix C lensing potential and convergence

First order strong lensing effects, i.e. image distortion and magnification
can be derived from the Jacobi matrix of lens equation:

$[A]_{ij}:=\left(\frac{\partial\vec{y}}{\partial\vec{x}}\right)_{ij}=\left(\delta_{ij}-\frac{\partial^{2}\Psi}{\partial
x_{i}\partial x_{j}}\right)$ (114)

With a shorthand notation $\Psi_{ij}$ for partial derivatives in the above
expression, one can express the shear tensor
$\Gamma=\begin{pmatrix}\gamma_{1}&\gamma_{2}\\\
\gamma_{2}&-\gamma_{1}\end{pmatrix}$ where:
$\gamma_{1}=\frac{1}{2}\left(\Psi_{11}-\Psi_{22}\right)=:\gamma\cos{2\varphi}$,
$\gamma_{2}=\Psi_{12}=\Psi_{21}=:\gamma\sin{2\varphi}$ For the SIS lens in
Palatini theory, one has

$\displaystyle\Psi_{11}$ $\displaystyle=$
$\displaystyle\frac{1}{x^{5}}\left[x^{2}\left(x^{2}-\frac{2\alpha}{\xi_{0}^{2}}\right)+\left(\frac{6\alpha}{\xi_{0}^{2}}-x^{2}\right)x_{1}^{2}\right]$
$\displaystyle\Psi_{12}$ $\displaystyle=$
$\displaystyle\Psi_{21}=\frac{x_{1}x_{2}}{x^{5}}\left(\frac{6\alpha}{\xi_{0}^{2}}-x^{2}\right)$
$\displaystyle\Psi_{22}$ $\displaystyle=$
$\displaystyle\frac{1}{x^{5}}\left[x^{2}\left(x^{2}-\frac{2\alpha}{\xi_{0}^{2}}\right)+\left(\frac{6\alpha}{\xi_{0}^{2}}-x^{2}\right)x_{2}^{2}\right]\,,$

hence the shear components are

$\displaystyle\gamma_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2x^{5}}\left(\frac{6\alpha}{\xi_{0}^{2}}-x^{2}\right)\left(x_{1}^{2}-x_{2}^{2}\right)$
$\displaystyle\gamma_{2}$ $\displaystyle=$
$\displaystyle\left(\frac{6\alpha}{\xi_{0}^{2}}-x^{2}\right)\frac{x_{1}x_{2}}{x^{5}}\,,$
(115)

where the notation $\vec{x}=(x_{1},x_{2})$ and $x=\sqrt{x_{1}^{2}+x_{2}^{2}}$
was used.