ƒ(R) Gravity, Relic Coherent Gravitons and Optical Chaos

Abstract

We discuss the production of massive relic coherent gravitons in a particular class of $f\left(R\right)$ gravity, which arises from string theory, and their possible imprint in the Cosmic Microwave Background. In fact, in the very early Universe, these relic gravitons could have acted as slow gravity waves. They may have then acted to focus the geodesics of radiation and matter. Therefore, their imprint on the later evolution of the Universe could appear as filaments and a domain wall in the Universe today. In that case, the effect on the Cosmic Microwave Background should be analogous to the effect of water waves, which, in focusing light, create optical caustics, which are commonly seen on the bottom of swimming pools. We analyze this important issue by showing how relic massive gravity waves (GWs) perturb the trajectories of the Cosmic Microwave Background photons (gravitational lensing by relic GWs). The consequence of the type of physics discussed is outlined by illustrating an amplification of what might be called optical chaos.

1. Introduction

Modified gravity currently obtains a lot of attention from the scientific community. The main reason is the remarkable issue that it enables a description of early-time inflation, as well as late-time acceleration epoch (dark energy) in a unified way.

In recent years, superstring/Mtheory caused a lot of interest about higher order gravity in more than four dimensions [1]. These models work in the effective low-energy action of superstring theory [1,2]. Within the classical framework, they have to be inserted among the class of the so-called $f\left(R\right)$ theories of gravity (for a recent review, see [3]).

Motivations for a potential extension of Einstein’s general relativity (GR) [4] are various. First of all, as distinct from other field theories, like the electromagnetic theory, GR is very difficult to quantize. This fact rules out the possibility of treating gravitation like other quantum theories and precludes the unification of gravity with other interactions. At the present time, it is not possible to realize a consistent quantum gravity theory that leads to the unification of gravitation with the other forces. One of the most important goals of modern physics is to obtain a unified theory, which could, in principle, show the fundamental interactions as different forms of the same symmetry. Considering this point of view, today, one observes and tests the results of one or more breaks of symmetry. In this way, it is possible to say that we live in an unsymmetrical world [5]. In the last 60 years, the dominant idea has been that a fundamental description of physical interactions arises from quantum field theory [6]. In this approach, different states of a physical system are represented by vectors in a Hilbert space defined in a spacetime, while physical fields are represented by operators (i.e., linear transformations) on such a Hilbert space. The greatest problem is that this quantum mechanical framework is not consistent with gravitation, because this particular field, i.e., the metric, $g_{\mu \nu }$, describes both the dynamical aspects of gravity and the spacetime background [5]. In other words, one says that the quantization of dynamical degrees of freedom of the gravitational field is meant to give a quantum-mechanical description of the spacetime. This is an unequaled problem in the context of quantum field theories, because the other theories are founded on a fixed spacetime background, which is treated like a classical continuum. Thus, at the present time, an absolute quantum gravity theory, which implies a total unification of various interactions, has not been obtained [5]. In addition, GR assumes a classical description of the matter, which is totally inappropriate at subatomic scales, which are the scales of the early Universe [3,5].

In the general context of cosmological evidence, there are also other considerations that suggest an extension of GR [3,7]. As a matter of fact, the accelerated expansion of the Universe, which is observed today, implies that cosmological dynamics is dominated by the so-called dark energy, which gives a large negative pressure. This is the standard picture, in which this new ingredient should be some form of unclustered, non-zero vacuum energy, which, together with the clustered dark matter, drives the global dynamics. This is the so-called “concordance model” (ΛCDM), which gives, in agreement with the Cosmic Microwave Background Radiation, Large Scale Structure and Supernovae Ia data, a good picture of the observed Universe today, but presents several shortcomings, such as the well-known “coincidence” and “Cosmological Constant” problems [8].

An alternative approach is seeing if the observed cosmic dynamics can be achieved through an extension of GR [3,7]. In this different context, it is not required to find candidates for dark energy and dark matter that, till now, have not been found; only the “observed” ingredients, which are curvature and baryon matter, have to be taken into account. Then, dark energy and dark matter have to be considered like pure effects of the presence of an intrinsic curvature in the Universe. Considering this point of view, one can think that gravity is different at various scales, and there is room for alternative theories.

Note that we are not claiming that GR is wrong. It is well known that, even in the context of extended theories of gravity, GR remains the most important part of the structure [7]. We are only trying to understand if weak modifications of such a structure could be needed to solve some theoretical and current observational problems. In this picture, we also recall that even Einstein tried to modify the framework of GR by adding the “Cosmological Constant” [9]. In any case, cosmology and Solar System tests show that modifications of GR in the sense of extended theories of gravity have to be very weak [3,7].

In principle, the most popular dark energy and dark matter models can be achieved in the framework of extended theories of gravity, i.e., $f\left(R\right)$ theories of gravity [3] and scalar tensor theories of gravity [7], which are generalizations of the Jordan–Fierz–Brans–Dicke Theory [10,11,12]. One assumes that geometry (for example, the Ricci curvature scalar, R) interacts with material quantum fields, generating back-reactions, which modify the gravitational action, adding interaction terms (examples are high-order terms in the Ricci scalar and/or in the Ricci tensor and non-minimal coupling between matter and gravity). This approach enables the modification of the Lagrangian, with respect to the standard Einstein–Hilbert gravitational Lagrangian [13], through the addition of high-order terms in the curvature invariants (terms like $R^2$, $R^{\alpha \beta }{R}_{\alpha \beta }$, $R^{\alpha \beta \gamma \delta }{R}_{\alpha \beta \gamma \delta }$, $R\square R$ and $R{\square }^{k}R$, in the sense of $f\left(R\right)$ theories [3,7]) and/or terms with scalar fields non-minimally coupled to geometry (terms like ${\varphi }^{2}R$) in the sense of scalar-tensor theories [7].

In the tapestry of $f\left(R\right)$ theories, the higher order terms are physically a type of back reaction from geometry acting upon matter, which further modifies geometry. This is a topological massive gravity, which represents a form of intrinsic curvature to spacetime. These terms are related to the Bel–Robinson tensor [14]:

$${T^{\mu }}_{\nu \sigma \rho }={R^{\mu \alpha \beta }}_{\sigma }{R}_{\nu \alpha \beta \rho }+R^{\mu \alpha \beta }{R}_{\nu \alpha \beta \rho }-\frac{1}{2}{{\delta }^{\mu }}_{\nu }{R^{\alpha \beta \gamma }}_{\sigma }{R}_{\alpha \beta \gamma \rho }$$

(1)

Contraction over indices gives the result,

$${{\delta }^{\sigma }}_{\mu }{g}^{\nu \rho }{T^{\mu }}_{\nu \sigma \rho }=R^{\mu \alpha }{R}_{\mu \beta }+R^{\mu \alpha \beta }{R^{\nu }}_{\alpha \beta \mu }-\frac{1}{2}{R}^{\alpha \beta \gamma \nu }{R}_{\alpha \beta \gamma \nu }$$

(2)

The physical consequences of this extension to curvature are fairly remarkable. The Bel–Robinson tensor is a vacuum curvature $\nabla T=0$, and it predicts gravity waves (GWs).

2. Gravity Waves in $f\left(R\right)$ Theories

In $f\left(R\right)$ gravity, the GWs have longitudinal structure [7,15,16], which makes them comparable to acoustical waves in a media. The linearized theory of weak GWs with a metric perturbation [15,16,17]:

$$g_{mu\nu }={\eta }_{\mu \nu }+h_{\mu \nu }$$

(3)

gives a traceless solution in standard GR [17]:

$$\square \tilde{h}=0$$

(4)

The modified gravity results in a terms that acts as a mass, where the wave equation is [15,16]:

$$\square \tilde{h}=m^2\tilde{h}$$

(5)

The decomposition of the solution gives the standard $h_{++}$ and $h_{\times \times }$ polarization modes; the mass introduces a third polarization, which is a longitudinal mode [15,16].

Let us consider a string theory setting [1]. The gravitational action is expanded in powers of ${\alpha }^n{R}^{2n}$ [2], for α, the string parameter. The action is [2]:

$$S=\int \left[\frac{1}{2\kappa }R+{\alpha }^{\prime }{R}^{\mu \nu \sigma \rho }{R}_{\mu \nu \sigma \rho }+L\right]\sqrt{-g}{d}^{4}x$$

(6)

being L the Lagrangian for everything else compactified on a $Dp$-brane. This action may be trivially rewritten as standard $R^2$ gravity [18]. Following the advice of [18], by using the Gauss Bonnet identity (its invariant) [19,20,21,22,23,24,25,26,27,28,29,30,31], one can indeed express a Riemann tensor squared term as a combination of pure $R^2$ and a Ricci tensor squared term. In such form, the action (Equation (6)) has been already studied by number of researchers; see [3,32,33,34,35] and the references within. This is a key point. In fact, although in the form of Equation (6), this action looks to be neither a renormalizable nor ghost-free theory [18], by using the Gauss Bonnet invariant, one can reduce it to the simpler form of $R^2$ gravity, which is the simplest one among the class of viable models with $R^m$ terms in addition to the Einstein–Hilbert theory. In [35], it has been shown that such models may lead to the (Cosmological Constant or quintessence) acceleration of the Universe, as well as an early-time era of inflation. Moreover, they seem to pass the Solar System tests, i.e., they have the acceptable Newtonian limit, no instabilities and no Brans–Dicke problem (decoupling of the scalar) in the scalar-tensor version.

The extremization of Equation (6) gives:

$$\begin{array}{cc}\hfill \delta S=0=& \int [\frac{1}{2\kappa }\frac{\delta R}{\delta g^{\mu \nu }}+2\alpha R^{\mu \nu \sigma \rho }\frac{\delta R_{\mu \nu \sigma \rho }}{\delta g^{\mu \nu }}+\frac{\delta L}{\delta g^{\mu \nu }}]\sqrt{-g}{d}^{4}x\hfill \\ \hfill & +\int \left[\frac{1}{2\kappa }R+\alpha R^{\mu \nu \sigma \rho }{R}_{\mu \nu \sigma \rho }+L\right]\frac{2}{\sqrt{-g}}\frac{\delta (-g)}{\delta g^{\mu \nu }}{d}^{4}x\hfill \end{array}$$

(7)

This action derives a modified form of the Einstein field equation:

$$R_{\mu \nu }+\frac{R}{2}{g}_{\mu \nu }+\alpha R_{\mu \nu \sigma \rho }{g}^{\sigma \rho }=\kappa T_{\mu \nu }$$

(8)

Now, let us consider a metric with the form of Equation (3), which is a classical expression. The quadratic term corresponds to quantum corrections on the order of the parameter, α. We consider this correction due to fields ${\varphi }_{\nu }^{\mu }$, so that the quantum correction to the metric is:

$$g_{\mu \nu }={\eta }_{\mu \nu }+h_{\mu \nu }+{\varphi }_{\mu }^{\sigma }{\varphi }_{\nu \sigma }$$

(9)

where we can regard ${\varphi }_{\mu }^{\sigma }{\varphi }_{\nu \sigma }=\delta h_{\mu \nu }$. These fields are physically a quantum correction to the classical gravitational radiation, $h_{\mu \nu }$. In general, these fields are quantized fields. In a string theory framework [1], we may define operators of the form:

$${\varphi }^{\mu \nu }=\sum _{m,n=1}^{\infty }{\alpha }_{m-n}^{\mu }{\alpha }_n^{\nu }$$

(10)

which is a harmonic oscillator quantization condition compatible with a string theory interpretation [1]. The graviton fields are given by the $n=m-n=-1$ states:

$${\varphi }_{\mu }^{\sigma }{\varphi }_{\nu \sigma }={\alpha }_{-1}^{\mu }V\left(x\right)e^{ikx}{\alpha }_{-1}^{\nu }V\left(x\right)e^{ik{x}^{\prime }}$$

(11)

such that ${\alpha }_{-1}^{\mu }{\alpha }_{-1}^{\nu }|0\rangle =|{\omega }_{\mu \nu }\rangle$ constructs the elementary states.

The connection terms are computed as:

$${\omega }_{\nu \sigma }^{\mu }={\partial }_{\nu }{\varphi }_{\rho }^{\mu }{\varphi }_{\sigma }^{\rho }$$

(12)

where the field is treated as a vierbein. Now, let us compute the curvature as:

$${R^{\mu }}_{nu\sigma \rho }={\partial }_{\sigma }{\omega }_{\nu \rho }^{\mu }-{\partial }_{\rho }{\omega }_{\nu \sigma }^{\mu }+{[{\omega }_{\sigma },{\omega }_{\rho }]}_{\nu }^{\mu }$$

(13)

The connection terms are on the order of the fundamental length, α, which is small enough to ignore the second order term. The quantized graviton field may then be written in this linearized fashion as:

$$\begin{array}{ccc}\hfill R_{\mu \nu \sigma \rho }& =\hfill & {\partial }_{\sigma }{\omega }_{\nu \rho }^{\mu }-{\partial }_{\rho }{\omega }_{\nu \sigma }^{\mu }={\partial }_{\sigma }\left({\partial }_{\nu }{\varphi }_{\gamma }^{\mu }{\varphi }_{\rho }^{\gamma }\right)-{\partial }_{\rho }\left({\partial }_{\nu }{\varphi }_{\gamma }^{\mu }{\varphi }_{\sigma }^{\gamma }\right)\hfill \\ \hfill & =\hfill & \left({\partial }_{\sigma }{\partial }_{\nu }{\varphi }_{\gamma }^{\mu }\right){\varphi }_{\rho }^{\gamma }-\left({\partial }_{\rho }{\partial }_{\nu }{\varphi }_{\gamma }^{\mu }\right){\varphi }_{\sigma }^{\gamma }+{\partial }_{\nu }{\varphi }_{\gamma }^{\mu }{\partial }_{\sigma }{\varphi }_{\rho }^{\gamma }-{\partial }_{\nu }{\varphi }_{\gamma }^{\mu }{\partial }_{\rho }{\varphi }_{\sigma }^{\gamma }\hfill \end{array}$$

(14)

where the last term is zero in a linearized approximation.

The linearized approximation occurs for long wavelength gravitons. Assume the connection term ${\omega }_{\nu \sigma }^{\mu }={\partial }_{\nu }{\varphi }_{\rho }^{\mu }{\varphi }_{\sigma }^d$ is eigenvalued with a wave number, $k^a$:

$${\omega }_{\nu \sigma }^{\mu }=k_{\nu }{\varphi }_{\rho }^{\mu }{\varphi }_{\sigma }^d$$

(15)

so the curvature tensor is:

$$\begin{array}{ccc}\hfill R_{\mu \nu \sigma \rho }& \simeq \hfill & {\partial }_{\sigma }{\omega }_{\nu \rho }^{\mu }-{\partial }_{\rho }{\omega }_{\nu \sigma }^{\mu }\hfill \\ \hfill & =\hfill & \left(k_{\sigma }{k}_{\nu }\right){\varphi }_{\gamma }^{\mu }{\varphi }_{\rho }^{\gamma }-k_{\rho }{k}_{\nu }{\varphi }_{\gamma }^{\mu }{\varphi }_{\sigma }^{\gamma }=\left(k_{\sigma }{k}_{\nu }\right)\delta h_{\rho }^{\mu }-k_{\rho }{k}_{\nu }\delta h_{\sigma }^{\mu }\hfill \end{array}$$

(16)

The second order term in the action is then:

$$R_{\mu \nu \sigma \rho }{R}^{\mu \nu \sigma \rho }\simeq [\left(k_{\sigma }{k}_{\nu }\right)\delta h_{\rho }^{\mu }-\left(k_{\rho }{k}_{\nu }\right)\delta h_{\sigma }^{\mu }][\left(k^{\sigma }{k}^{\nu }\right)\delta h_{\sigma \nu }-\left(k^{\rho }{k}^{\nu }\right)\delta h_{\mu \sigma }]=6k^4$$

(17)

The term, $\alpha k^4$, is a quartic term in mass, where the string coupling constant, $\sim G_N$, has naturalized units of area.

This is an intrinsic curvature in spacetime. The string coupling constant is about $\alpha \sim {10}^{-60}$ cm${}^2$, which is a small number. This also guarantees the viability of the action (Equation (6)), because the theory can pass Solar System and cosmology tests [7].

Is it possible that this mass effect should then become apparent in the laboratory? The question is, what is the laboratory? The obvious laboratory is the Cosmic Microwave Background (CMB). In fact, we recall that relic gravitons should have been produced in the Inflationary Era. This is a consequence of general assumptions. Essentially, it derives from a mixing between basic principles of classical theories of gravity and of quantum field theory [36,37,38]. The strong variations of the gravitational field in the early Universe amplify the zero-point quantum oscillations and produce relic GWs. It is well known that the detection of relic GWs is the only way to learn about the evolution of the very early Universe, up to the bounds of the Planck epoch and the initial singularity [36,37,38]. It is very important to stress the unavoidable and fundamental character of this mechanism. The model derives from the inflationary scenario for the early Universe [37], which is tuned in a good way with the WMAPdata on the CMB (in particular, exponential inflation and spectral index $\approx 1$) [39,40]. Inflationary models of the early Universe were analyzed in the early and mid-1980s [37]. These are cosmological models in which the Universe undergoes a brief phase of a very rapid expansion in early times. In this context, the expansion could be power-law or exponential in time. Inflationary models provide solutions to the horizon and flatness problems [37] and contain a mechanism that creates perturbations in all fields [36,38]. Important for our goals is that this mechanism also provides a distinctive spectrum of relic GWs. The GWs perturbations arise from the uncertainty principle, and the spectrum of relic GWs is generated from the adiabatically-amplified zero-point fluctuations [36,38].

Relic gravitons can be characterized by a dimensionless spectrum [36,38]:

$${\Omega }_{gw}\left(f\right)\equiv \frac{1}{{\rho }_c}\frac{d{\rho }_{gw}}{dlnf}$$

(18)

where:

$${\rho }_c\equiv \frac{3H_0^2}{8G}$$

(19)

is the (actual) critical density energy, ${\rho }_c$, of the Universe, $H_0$ the actual value of the Hubble expansion rate and $d{\rho }_{gw}$ the energy density of relic GWs in the frequency range, f to $f+df$.

In the standard inflationary model, the spectrum is flat over a wide range of frequencies; see [36,38] and Figure 1. The more recent value for the flat part of the spectrum that arises from the WMAP data can be found in [38],

$${\Omega }_{gw}\left(f\right)\le 9\times {10}^{-13}$$

(20)

Based on the weakness of the signal, it will be very difficult to detect relic gravitons on Earth, but a potential detection could be, in principle, realized with LISA [15]. However, the presence of relic gravitons may have perturbed the early Universe in ways that might be observable in the fine details of the CMB background. These gravitons would introduce a small dispersion in GWs, which might then leave an imprint on the CMB. We will discuss the potential presence of such an imprint in next Section.

Now, let us expand the field, ${\varphi }_{\nu }^{\mu }$, according to harmonic oscillator operators, $b,b^{†}$, as a simple model of a string. The fields are expanded as:

$${\varphi }_{\nu }^{\mu }=\left(\frac{1}{\sqrt{2}}\right)\sum _k{E}_{\nu }^{\mu }b\left(k\right)e^{i\theta \left(k\right)}+b^{†}{e}^{-i\theta \left(k\right)}$$

(21)

where $E_{\nu }^{\mu }$ is a tetrad, which is discussed more below. The summation runs from $\{-\infty ,\infty \}$. The product ${\varphi }_a^c{\varphi }_{\nu \sigma }=\delta h_{\mu \nu }$ is a harmonic oscillator operator:

$$\begin{array}{cc}\hfill {\varphi }_{\nu }^{\mu }{\varphi }_{\sigma \mu }=& \left(\frac{1}{2}\right)\sum _{k{k}^{\prime }}{E}_{\nu \sigma }^{2}b\left(k\right)b^{†}\left(k^{\prime }\right)e^{i\theta \left(k\right)-i\theta \left(k^{\prime }\right)}+b^{†}\left(k\right)b\left(k^{\prime }\right)e^{-i\theta \left(k^{\prime }\right)-i\theta \left(k\right)}\hfill \\ \hfill & +\left(\frac{1}{2}\right)\sum _{k{k}^{\prime }}{E}_{\nu \sigma }^{2}b\left(k\right)b\left(k^{\prime }\right)e^{i\theta \left(k\right)+i\theta \left(k^{\prime }\right)}+b^{†}\left(k\right)b^{†}\left(k^{\prime }\right)e^{-i\theta \left(k\right)-i\theta \left(k^{\prime }\right)}\hfill \end{array}$$

(22)

Figure 1. The spectrum of relic scalar GWs in inflationary models is flat over a wide range of frequencies. The horizontal axis is ${\log }_{10}$ of frequency, in hertz. The vertical axis is ${\log }_{10}{\Omega }_{gsw}$. The inflationary spectrum rises quickly at low frequencies (the wave that re-entered in the Hubble sphere after the Universe became matter dominated) and falls off above the (appropriately redshifted) frequency scale, $f_{max}$, associated with the fastest characteristic time of the phase transition at the end of inflation. The amplitude of the flat region depends only on the energy density during the inflationary stage; we have chosen the largest amplitude consistent with the WMAP constraints on scalar perturbations. This means that at LIGOand LISAfrequencies, ${\Omega }_{gw}\left(f\right)h_{100}^2<9\times {10}^{-13}$. Adapted from [41].

Figure 1. The spectrum of relic scalar GWs in inflationary models is flat over a wide range of frequencies. The horizontal axis is ${\log }_{10}$ of frequency, in hertz. The vertical axis is ${\log }_{10}{\Omega }_{gsw}$. The inflationary spectrum rises quickly at low frequencies (the wave that re-entered in the Hubble sphere after the Universe became matter dominated) and falls off above the (appropriately redshifted) frequency scale, $f_{max}$, associated with the fastest characteristic time of the phase transition at the end of inflation. The amplitude of the flat region depends only on the energy density during the inflationary stage; we have chosen the largest amplitude consistent with the WMAP constraints on scalar perturbations. This means that at LIGOand LISAfrequencies, ${\Omega }_{gw}\left(f\right)h_{100}^2<9\times {10}^{-13}$. Adapted from [41].

The sum gives a delta function on k and $k^{\prime }$, and the first term is the Hamiltonian, which after the use of a commutator the RHSterm is:

$${\varphi }_{\nu }^{\mu }{\varphi }_{\sigma \mu }=\left(\frac{1}{2}\right)\sum _k{E}_{\nu \sigma }^2{b}^{†}\left(k\right)b\left(k\right)+\left(\frac{1}{2}\right)\sum _k{E}_{\nu \sigma }^{2}b\left(k\right)b(-k)+b^{†}\left(k\right)b^{†}(-k)$$

(23)

where the zeta point energy (ZPE) term has been dropped. The first RHS term is a familiar Hamiltonian type of term, while the second term is similar to a squeeze operator in quantum optics [42].

The tetrad, $E_{\nu }^{\mu }$, is the amplitude of the field. This plays a role similar to the minimal electric field $E=\sqrt{\hslash \omega /V{ϵ}_0}$ in box normalization [43]. A plausible choice for the tetrad is then $E_{\nu }^{\mu }=\sqrt{\alpha \omega }{\delta }_{\nu }^{\mu }$, where α is the string parameter and ω the frequency. For $\alpha \ll 1/\omega$, this is a small term.

The curvature in quantum modes is then:

$$R_{\mu \nu \sigma \rho }\simeq \left(\frac{1}{2}\right)\sum _k\left(k_{\sigma }{k}_{\nu }{\left(E^2\right)}_{\mu }^{\beta }{E}_{\beta \rho }^2-k_{\rho }{k}_{\nu }{E}_{\mu \beta }^2{\left(E^2\right)}_{\sigma }^{\beta }\left)\right(b^{†}\left(k\right)b\left(k\right)+b\left(k\right)b(-k)+b^{†}\left(k\right)b^{†}(-k)\right)$$

(24)

which is $O\left(\alpha \right)$ in the scale parameter. In fluctuations of the curvature, the metric is $g\sim \delta L/L$, for $\delta L>L_p$. The connection terms are of order $\Gamma \sim \delta L/L^2$ and curvatures are $R\sim \delta L/L^3$. The wave vectors are $k\sim 1/L$, and the scaling parameter is $\alpha \omega \sim \delta L/L$. From a dimensional and scaling perspective, this answer appears at least proximal.

For the sake of simplicity, let us write the curvature tensor as:

$$R_{\mu \nu \sigma \rho }\simeq \left(\frac{1}{2}\right)\sum _k{\Pi }_{\mu \nu \sigma \rho }\left(k\right)\left(b^{†}\left(k\right)b\left(k\right)+b\left(k\right)b(-k)+b^{†}\left(k\right)b^{†}(-k)\right)$$

(25)

The second order term is formed from the total contraction on the Riemann tensor:

$$\begin{array}{c}{R}_{\mu \nu \sigma \rho }{R}^{\mu \nu \sigma \rho }\simeq \left(\frac{1}{4}\right)\sum _{k{k}^{\prime }}{\Pi }_{\mu \nu \sigma \rho }\left(k\right){\Pi }^{\mu \nu \sigma \rho }\left(k^{\prime }\right)x\\ (b^{†}\left(k\right)b\left(k\right)b^{†}\left(k^{\prime }\right)b\left(k^{\prime }\right)+b^{†}\left(k\right)b\left(k\right)(b\left(k^{\prime }\right)b(-k^{\prime }+b^{†}\left(k^{\prime }\right)b^{†}(-k^{\prime }))+(b\left(k^{\prime }\right)b(-k^{\prime })+\\ b^{†}\left(k^{\prime }\right)b^{†}(-k^{\prime }))b^{†}\left(k\right)b\left(k\right)+\left(b\left(k\right)b(-k)+b^{†}\left(k\right)b^{†}(-k)\left)\right(b\left(k^{\prime }\right)b(-k^{\prime })+b^{†}\left(k^{\prime }\right)b^{†}(-k^{\prime })\right))\end{array}$$

(26)

This term is to $O\left({\alpha }^2\right)$ and contributes a term $O\left({\alpha }^3\right)$ to the Lagrangian.

Consider the operator matrix operation $R_{\mu \nu \sigma \rho }{R}^{\mu \nu \sigma \rho }|m\rangle$. The first term has the operator matrix elements:

$$\begin{array}{ccc}\hfill b^{†}\left(k\right)b\left(k\right)b^{†}\left(k^{\prime }\right)b\left(k^{\prime }\right)|m\rangle & =\hfill & b^{†}\left(k\right)\sum _n|n\rangle \langle n|b\left(k\right)b^{†}\left(k^{\prime }\right)b\left(k^{\prime }\right)|m\rangle \hfill \\ \hfill & =\hfill & m\left(k^{\prime }\right)n\left(k\right){\delta }_{mn}{\delta }_{k{k}^{\prime }}\hfill \end{array}$$

(27)

where ${\sum }_n|n\rangle \langle n|$ is a completeness sum and the momentum values assumed in the states, $|m\rangle$ and $|n\rangle$. This contributes an energy-squared. A similar analysis for $\langle m|R_{\mu \nu \sigma \rho }{R}^{\mu \nu \sigma \rho }$ gives:

$$\langle m|b^{†}\left(k\right)b\left(k\right)(b\left(k^{\prime }\right)b(-k^{\prime })+b^{†}\left(k^{\prime }\right)b^{†}(-k^{\prime }))=m\left(k\right)(b\left(k^{\prime }\right)b(-k^{\prime })+b^{†}\left(k^{\prime }\right)b^{†}(-k^{\prime }))$$

(28)

and for $R_{\mu \nu \sigma \rho }{R}^{\mu \nu \sigma \rho }|m\rangle$,

$$(b\left(k\right)b(-k)+b^{†}\left(k\right)b^{†}(-k))b^{†}\left(k^{\prime }\right)b\left(k^{\prime }\right)|m\rangle =m(b\left(k\right)b(-k)+b^{†}\left(k\right)b^{†}(-k))|m\rangle$$

(29)

The operators, $b\left(k\right)b(-k)+b^{†}\left(k\right)b^{†}(-k)$, form the squeeze operator [44]:

$$S=exp\left(\left(\frac{1}{2}\right)(z^{*}b\left(k\right)-z{b}^{†}\left(k\right))\right)$$

(30)

where $z^{*}=z=i\left((1/4){\Pi }_{\mu \nu \sigma \rho }\left(k\right){\Pi }^{\mu \nu \sigma \rho }\right)$. Hence, the action phase due to the action, $e^{iS}$, contains a squeeze operator.

The final operator term is more complicated. The operator terms, $b\left(k\right)b(-k)b\left(k\right)b(-k)$ and $b^{†}\left(k\right)b^{†}(-k)b\left(k\right)b(-k)$, are evaluated by commuting operators, and this leads to the square of number operators, $n\left(k\right)n(-k)$. The terms, $b\left(k\right)b(-k)b\left(k\right)b(-k)$ and $b^{†}\left(k\right)b^{†}(-k)b^{†}\left(k\right)b^{†}(-k)$, are then a product of terms that represent a squeezed state operator.

The squeeze operator, $S\left(z\right)$, acts upon the displacement operator $D\left(\alpha \right)=exp(\alpha b^{†}-{\alpha }^{*}b)$, so that $S\left(z\right)D\left(\alpha \right)\ne D\left(\alpha \right)S\left(z\right)$,

$$\begin{array}{ccc}\hfill S\left(z\right)D\left(\alpha \right)& =\hfill & exp[(z^{*}{b}^2-z{\left(b^{†}\right)}^2)/2]exp(\alpha b^{†}-{\alpha }^{*}b)\hfill \\ \hfill & =\hfill & exp[(z^{*}{b}^2-z{\left(b^{†}\right)}^2)/2+\alpha b^{†}-{\alpha }^{*}b)exp[-\left(\frac{1}{4}\right)(z^{*}\alpha b^{†}-z{\alpha }^{*}b)]\hfill \end{array}$$

(31)

which effectively creates a modified displacement operator,

$$exp[(z^{*}{b}^2-z{\left(b^{†}\right)}^2)/2+\alpha b^{†}-{\alpha }^{*}b)S\left(z\right)D\left(\alpha \right)=exp[-\left(\frac{1}{4}\right)(z^{*}\alpha b^{†}-z{\alpha }^{*}b)]=D\left(z^{*}\alpha \right)$$

(32)

The action of the squeezed state operator on b is $Sb{S}^{†}=bcosh\left(\right|z\left|\right)+b^{†}sinh\left(\right|z\left|\right)$, which is a Bogoliubov transformed operator [45]. For a set of bosons, here, linear gravitons, with the same state, there is then ${\sum }_n{\alpha }^n/\sqrt{n}|n\rangle$ states with the operator acting on this ${\sum }_n{\alpha }^n{\left(a^{†}\right)}^n/\sqrt{n}$ acting on the vacuum. This operator may be formed from the $S\left(z\right)D\left(\alpha \right)S^{†}\left(z\right)$ for α small and $\left|z\right|\gg \left|\alpha \right|$ with,

$$S\left(z\right)D\left(\alpha \right)S^{†}\left(z\right)\simeq exp[-\left(\frac{1}{4}\right)(z^{*}\alpha b^{†}-z{\alpha }^{*}b)]$$

(33)

and where we may then define $z^{*}\alpha /4\to \alpha$, and the Bogoliubov transformation of the operator, $b-b^{†}$, constructs a displacement operator [46].

In this way, the $R^2$ term in the action describes the squeezed state operator, which acts on the field raising and lowering operators to define a displacement operator for coherent states, which, in the case of photons, are laser states of light.

We then evaluate a Wilson loop [47] $W\left({\varphi }_{\nu }^{\mu }\right)=exp(\int i{\varphi }_{\mu \nu }{e}^{\mu }d{x}^{\nu }$). In the path integral,

$$Z[\varphi ,W]=\int \mathcal{D}\left[\varphi \right]W{e}^{iS\left[\varphi \right]}$$

(34)

The infinitesimal shift in the field $\varphi \to \varphi +\delta \varphi$ adjusts $Z[\varphi ,W]\to Z[\varphi +\delta \varphi ,W]=\langle W\rangle$, and the expansion is:

$$\begin{array}{ccc}\hfill \langle W\rangle & =\hfill & \int D\left[\varphi \right]W(\varphi +\delta \varphi )e^{iS[\varphi +\delta \varphi ]}\hfill \\ \hfill & =\hfill & \langle W\rangle +\int D\left[\varphi \right]\delta \varphi \left(\frac{\delta W}{\delta \varphi }+i\frac{W\delta S}{\delta \varphi }\right)W{e}^{iS\left[\varphi \right]})\hfill \end{array}$$

(35)

where the invariance of the expectation gives:

$$\frac{\delta W}{\delta \varphi }+i\frac{W\delta S}{\delta \varphi }=0\to \frac{\delta ln\left(W\right)}{\delta \varphi }+i\frac{\delta S}{\delta \varphi }=0$$

(36)

This formula is only well defined for a polynomial function. Therefore, we make the following approximation. The loop is considered to be very small, and in that way, we can approximate the Wilson loop with:

$$W\left(\varphi \right)=1+i{\varphi }_{\mu \nu }{e}^{\mu }\delta x^{\nu }$$

(37)

so that the functional derivative of $W\left(\varphi \right)$ is:

$$\frac{\delta W}{\delta {\varphi }^{\mu \nu }}\simeq i{ϵ}_{\sigma \nu }{\delta }_{\mu }^{\nu \sigma }\delta (x-x^{\prime })$$

(38)

for ${ϵ}_{\mu }^{\nu }$, a unit area. The solution is then,

$$〈\frac{\delta W}{\delta \varphi }+iW\frac{\delta S}{\delta \varphi }〉=0\to 〈W\frac{\delta S}{\delta {\varphi }^{\mu \nu }}〉\simeq {ϵ}_{\sigma \nu }{\delta }_{\mu }^{\sigma }\delta (x-x^{\prime })$$

(39)

Now, let us consider the second order expansion,

$$W\left(\varphi \right)=1+i{\varphi }_{\mu \nu }{e}^{\mu }\delta x^{\nu }+\frac{1}{2}{\varphi }_{\mu \nu }{\varphi }_{\sigma \rho }{e}^{\mu }{e}^{\sigma }\delta x^{\nu }\delta x^{\rho }=W^0\left(\varphi \right)+W^1\left(\varphi \right)+W^2\left(\varphi \right)$$

(40)

which gives the result:

$$〈W^2\frac{\delta S}{\delta {\varphi }^{\mu \nu }}〉\simeq -\frac{i}{2}{ϵ}_{\sigma \nu }\langle {\varphi }_{\mu }^{\sigma }\rangle \delta (x-x^{\prime })$$

(41)

where, by continuing the series, this leads to:

$$〈W\frac{\delta S}{\delta {\varphi }^{\mu \nu }}〉\simeq i{e}^{i\langle {\varphi }_{\mu \nu }\rangle e^{\mu }\delta x^{\nu }}\delta (x-x^{\prime })=i\langle e^{i{\varphi }_{\mu \nu }{e}^{\mu }\delta x^{\nu }}\rangle \delta (x-x^{\prime })$$

(42)

The input of an expansion of the field, ϕ, results in the expectation of an operator with the form of the displacement operator.

It is now important to understand the form the fields in the expansion in $D\left(\alpha \right)$. The Wilson loop is a form of the Stokes’ law [48] and:

$$-iln\left(W\right)={\int }_{ϵ}{\partial }_{\alpha }{\varphi }_{\mu \nu }{e}^{\mu }d{ϵ}^{\alpha \nu }$$

(43)

In vacuum, the canonical $h_{++}$ and $h_{\times \times }$ polarizations obey $\square h_{++}=\square h_{\times \times }=0$ [17]. The longitudinal modes due to $R^2$ terms obeys [15,16]:

$$\square h_c=m^2{h}_c$$

(44)

where the mass is a topologically-induced mass. The longitudinal $h_c={\varphi }^2$ then defines the equation [15,16,38]:

$$\square {\varphi }_{mu\nu }=m^2{\varphi }_{mu\nu }$$

(45)

where a Lorenz gauge sets terms with $\square \varphi =0$ [15,16,38]. This term plays a role similar to the Helmholtz potential in electromagnetism:

$$\Phi =\frac{1}{4\pi {ϵ}_0}{\int }_V{d}^{3}r\frac{\rho \left(\overrightarrow{r}\right)}{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}$$

(46)

but in the case of $f\left(R\right)$ theories, it results an effective potential through the identifications [15,38]:

$$\begin{array}{ccccc}\Phi \to f^{\prime }\left(R\right)& & \mathrm{and}& & \frac{dV}{d\Phi }\to \frac{2f\left(R\right)-R{f}^{\prime }\left(R\right)}{3}\end{array}$$

(47)

which give a Klein–Gordon equation for the effective Φ scalar field [15,38]:

$$\square \Phi =\frac{dV}{d\Phi }$$

(48)

The ${\varphi }_{\mu \nu }$, which physically contributes to the Wilson integral, has a source term, which is the topological mass.

3. Potential Imprint in Cosmic Microwave Background

We recall that the CMB is thermal radiation filling the observable Universe almost uniformly [39,40]. Precise measurements of the CMB are fundamental for cosmology, because any viable proposed model of the Universe must explain this radiation. The CMB has a thermal black body spectrum at a temperature of $\sim 2.7$ K [39,40]. At the present time, the best available data on the CMB arise from the Planck satellite [39,40], which has produced detailed all-sky observations over nine frequency bands between 30 and 857 GHz. According to the data, subtle fluctuations in the CMB temperature were imprinted on the deep sky during the recombination era, i.e., when the Universe was about $370,000$ years old. That imprint reflects ripples that arose from the early era, at about ${10}^{-30}$ s after the initial singularity. It is a common opinion that such ripples should give rise to the current cosmic structure of galactic clusters and dark matter.

The Planck satellite works within the Solar System, and to take into account weak potential effects on the CMB by relic massive GWs, we can use the weak field approximation (the linearized theory). In the linearized theory, the standard expansion $g_{\mu \nu }={\eta }_{\mu \nu }+h_{\mu \nu }$ with “small” $h_{\mu \nu }$ is performed in an asymptotically Cartesian coordinate system. This frame is the proper reference frame of a local observer, which we assume to be located in A within the Solar System. In other words, we assume that the spacetime within the Solar System is locally flat with respect to the global distribution of the CMB. Our goal is to understand how relic massive GWs perturb the trajectories of CMB photons between A and B. The global effect results in a particular gravitational lensing [49], due to relic massive GWs. Some clarifications are needed concerning this issue. In our linearized approach, gravitational lensing can be described in the local Lorentz frame perturbed by the first order post-Newtonian potential. Hence, one can define a refractive index [32,49]:

$$n\equiv 1+2\left|V\right|$$

(49)

In the usual geometrical optics, the condition $n>1$ implies that the light in a medium is slower than in vacuum [50]. Then, the effective speed of light in a gravitational field is expressed by [32,49,50]:

$$v=\frac{1}{n}\approx 1-2\left|V\right|$$

(50)

Thus, one can obtain the Shapiro delay [51] by integrating over the optical path between the source and the observer:

$${\int }_{source}^{observer}2\left|V\right|dl$$

(51)

The situation is analogous to the prism [50].

3.1. Gravitational Lensing in the Direction of the Propagating Gravity Wave

For the sake of simplicity, we assume that A and B are both located in the direction of the propagating massive GW, which we assume to be the z direction.

By using the proper reference frame of a local observer, the time coordinate, $x_0$, is the proper time of observer A, and the spatial axes are centered on A. In the special case of zero acceleration and zero rotation, the spatial coordinates, $x_j$, are the proper distances along the axes, and the frame of the local observer reduces to a local Lorentz frame [17]. The line element is [17]:

$$d{s}^2=-{\left(d{x}^0\right)}^2+{\delta }_{ij}d{x}^{i}d{x}^j+O\left(\right|x^j{|}^2)d{x}^{\alpha }d{x}^{\beta }$$

(52)

The connection between Newtonian theory and linearized gravity is well known [13]:

$$g_{00}=1+2V$$

(53)

V being the Newtonian potential. Let us consider the interval for photons propagating along the z-axis:

$$d{s}^2=g_{00}d{t}^2+d{z}^2$$

(54)

The condition for a null trajectory ($ds=0$) gives the coordinate velocity of the photons:

$$v_p^2\equiv {\left(\frac{dz}{dt}\right)}^2=1+2V(t,z)$$

(55)

which to the first order is well approximated by:

$$v_p\approx [1+V(t,z)]$$

(56)

Knowing the coordinate velocity of the photon, the propagation time for its traveling between A and B, which corresponds to the proper distance, AB, in the presence of the graviton, can be defined:

$$T_1\left(t\right)={\int }_{z_A}^{z_B}\frac{dz}{v_p}\approx T-{\int }_0^{T}V(t^{\prime },z)dz$$

(57)

where T represents the uniform propagation time of the photon between A and B (i.e., the proper distance between A and B in natural units), as if it were moving in a flat spacetime, i.e., in the absence of GW, and $t^{\prime }$ is the delay time, which corresponds to the unperturbed photon trajectory:

$$t^{\prime }=t-(T-z)$$

(58)

(i.e., t is the time at which the photon arrives in the position, T; so $T-z=t-t^{\prime }$).

In order to compute $T_1$, we need to know the Newtonian potential $V(t,z)$, which is generated by the massive GW. We recall that the effect of the gravitational force on test masses is described by the equation:

$$\ddot{x^i}=-{\tilde{R}}_{0k0}^i{x}^k$$

(59)

which is the equation for geodesic deviation in this frame [17]. ${\tilde{R}}_{0k0}^i$ is the linearized Riemann tensor [17].

On the other hand, with an opportune choice of the Lorenz gauge, the linearization process of $f\left(R\right)$ theories, which generates the third longitudinal mode $h_c=h_c(t-v_{G}z)$, enables a conformally flat line element [15,16,38]:

$$d{s}^2=[1+h_c(t-v_{G}z)](-d{t}^2+d{z}^2+d{x}^2+d{y}^2)$$

(60)

$v_G$ represents the group velocity of the massive GW. In fact, the velocity of every standard massless tensorial mode, ${\overline{h}}_{\mu \nu }$, is the light speed, c, but the dispersion law for the modes of $h_c$ is that of a massive field, which can be discussed like a wave-packet [15,16,38]. Furthermore, the group-velocity of a wave-packet of $h_c$ centered on $\overrightarrow{p}$ is [15,16,38]:

$$\overrightarrow{v_G}=\frac{\overrightarrow{p}}{\omega }$$

(61)

which is exactly the velocity of a massive particle with mass m (see Equation (44)) and momentum $\overrightarrow{p}$.

This group-velocity is a function of both of the mass and frequency of the wave-packet [15,16,38]:

$$v_G=\frac{\sqrt{{\omega }^2-m^2}}{\omega }$$

(62)

Even if the coordinates Equations (52) are different from the coordinates Equation (60), we recall that the linearized Riemann tensor is gauge invariant [17]. Hence, we can calculate it directly from Equation (60). Following [16] it is:

$${\tilde{R}}_{\mu \nu \alpha \beta }=\frac{1}{2}\{{\partial }_{\mu }{\partial }_{\beta }{h}_{\alpha \nu }+{\partial }_{\nu }{\partial }_{\alpha }{h}_{\mu \beta }-{\partial }_{\alpha }{\partial }_{\beta }{h}_{\mu \nu }-{\partial }_{\mu }{\partial }_{\nu }{h}_{\alpha \beta }\}$$

(63)

that, in the case of Equation (60), begins [16]:

$${\tilde{R}}_{0\gamma 0}^{\alpha }=\frac{1}{2}\{{\partial }^{\alpha }{\partial }_0{h}_c{\eta }_{0\gamma }+{\partial }_0{\partial }_{\gamma }{h}_c{\delta }_0^{\alpha }-{\partial }^{\alpha }{\partial }_{\gamma }{h}_c{\eta }_{00}-{\partial }_0{\partial }_0{h}_c{\delta }_{\gamma }^{\alpha }\}$$

(64)

the different elements are (only the non-zero ones will be written) [16]:

$${\partial }^{\alpha }{\partial }_0{h}_c{\eta }_{0\gamma }=\left\{\begin{array}{ccc}{\partial }_t^2{h}_c& for& \alpha =\gamma =0\\ -{\partial }_z{\partial }_t{h}_c& for& \alpha =3;\gamma =0\end{array}\right\}$$

(65)

$${\partial }_0{\partial }_{\gamma }{h}_c{\delta }_0^{\alpha }=\left\{\begin{array}{ccc}{\partial }_t^2{h}_c& for& \alpha =\gamma =0\\ {\partial }_t{\partial }_z{h}_c& for& \alpha =0;\gamma =3\end{array}\right\}$$

(66)

$$-{\partial }^{\alpha }{\partial }_{\gamma }{h}_c{\eta }_{00}={\partial }^{\alpha }{\partial }_{\gamma }{h}_c=\left\{\begin{array}{ccc}-{\partial }_t^2{h}_c& for& \alpha =\gamma =0\\ {\partial }_z^2{h}_c& for& \alpha =\gamma =3\\ -{\partial }_t{\partial }_z{h}_c& for& \alpha =0;\gamma =3\\ {\partial }_z{\partial }_t{h}_c& for& \alpha =3;\gamma =0\end{array}\right\}$$

(67)

$$-{\partial }_0{\partial }_0{h}_c{\delta }_{\gamma }^{\alpha }=\begin{array}{ccc}-{\partial }_z^2{h}_c& for& \alpha =\gamma \end{array}$$

(68)

By putting these results in Equation (64), one gets [16]:

$$\begin{array}{c}{\tilde{R}}_{010}^1=-\frac{1}{2}{\ddot{h}}_c\\ {\tilde{R}}_{010}^2=-\frac{1}{2}{\ddot{h}}_c\\ {\tilde{R}}_{030}^3=\frac{1}{2}\square h_c\end{array}$$

(69)

Let us put Equation (44) in the third of Equations (69). We obtain [16]:

$${\tilde{R}}_{030}^3=\frac{1}{2}{m}^2{h}_c$$

(70)

which shows that the field is not transversal.

In fact, Equation (59) implies [16]:

$$\ddot{x}=\frac{1}{2}{\ddot{h}}_c(t-v_{G}z)x$$

(71)

$$\ddot{y}=\frac{1}{2}{\ddot{h}}_c(t-v_{G}z)y$$

(72)

and,

$$\ddot{z}=-\frac{1}{2}{m}^2{h}_c(t-v_{G}z)z$$

(73)

Therefore, the effect of the mass is exactly the generation of a longitudinal force (in addition to the transverse one). Note that in the limit $m\to 0$, the longitudinal force vanishes.

Equivalently, we can say that there is a gravitational potential [16,17]:

$$V(\overrightarrow{r},t)=-\frac{1}{4}{\ddot{h}}_c(t-v_{G}z)[x^2+y^2]+\frac{1}{2}{m}^2{\int }_0^z{h}_c(t-v_{G}a)ada$$

(74)

which generates the tidal forces, and that the motion of the test mass is governed by the Newtonian equation [16,17]:

$$\ddot{\overrightarrow{r}}=-▽V$$

(75)

Now, we can use Equation (74) to compute $T_1$ in Equation (57). We get,

$$T_1\left(t\right)\approx T-{\int }_0^{T}V(t^{\prime },z)dz=T-\frac{1}{2}{m}^2{\int }_0^{T}dz{\int }_0^z{h}_c(t^{\prime }-v_{G}a)ada$$

(76)

Thus, the variation of the proper distance between A and B from its unperturbed value, T, which is due to the presence of the massive GW, $h_c$, is:

$$\begin{array}{ccc}\hfill \delta T_1\left(t\right)& \approx \hfill & \frac{1}{2}{m}^2{\int }_0^{T}dz{\int }_0^z{h}_c(t-T+a-v_{G}a)ada\hfill \\ \hfill & =\hfill & \frac{1}{4}{m}^2{\int }_0^T{h}_c(t-v_{G}z-T+z)dz-\frac{1}{4}{m}^2{\int }_0^T{\int }_0^z{h}_c^{\prime }(t-T+a-v_{G}a)z^{2}dadz\hfill \end{array}$$

(77)

Introducing the Fourier transform of $h_c$ defined by:

$${\tilde{h}}_c\left(\omega \right)={\int }_{-\infty }^{\infty }dt{h}_c\left(t\right)exp\left(i\omega t\right)$$

(78)

Equation (77) can be integrated in the frequency domain by using the Fourier translation and derivation theorems:

$$\frac{\delta {\tilde{T}}_1\left(\omega \right)}{T}=\mathrm{{\rm Y}}\left(\omega \right){\tilde{h}}_c\left(\omega \right)$$

(79)

where:

$$\begin{array}{c}\mathrm{{\rm Y}}\left(\omega \right)=\frac{1}{4}{m}^2\frac{exp\left(i\omega T\right)}{i\omega T(v_G-1)}\{expi\omega T(v_G-1)-1\\ +\frac{1}{i\omega (v_G-1)}\left[T^{2}expi\omega T(v_G-1)-2Texpi\omega T(v_G-1)+2expi\omega T(v_G-1)-1\right]-\frac{T^3}{3}\}\end{array}$$

(80)

is the longitudinal response function for relic gravitons.

In order to use Equations (79) and (80), we recall that relic gravitons represent a stochastic background [36,38]. Hence, one has to use average quantities [36,38]. The well-known equation for the characteristic amplitude [36], adapted for the third component of GWs, can be used [38]:

$$h_{cc}\left(f\right)\simeq 1.26\times {10}^{-18}\left(\frac{1Hz}{f}\right)\sqrt{h_{100}^2{\Omega }_{gw}\left(f\right)},$$

(81)

obtaining, for example, at 100 Hz and taking into account the bound (Equation (20)):

$$h_{cc}\left(100Hz\right)\simeq 1.7\times {10}^{-26}$$

(82)

Considering a graviton propagating with a speed of $v_G=0.99$ (ultra-relativistic case), if we insert these values in Equations (79) and (80), we get $\mathrm{{\rm Y}}\left(\omega \right)\approx 0.02$ and $\delta {\tilde{T}}_1\approx 3.4\times {10}^{-25}$ m for a proper distance between A and B of unperturbed value $T=1$ km. The situation is different for a speed of $0.9$ (relativistic case). In that case, one has $\mathrm{{\rm Y}}\left(\omega \right)\approx 0.19$ and $\delta {\tilde{T}}_1\approx 3.4\times {10}^{-24}$ m. For a speed of $0.1c$ (non relativistic case), we have $\mathrm{{\rm Y}}\left(\omega \right)\approx 0.99$ and $\delta {\tilde{T}}_1\approx 1.6\times {10}^{-23}$ m. The situation is better at lower frequencies. For $f=10$ Hz, Equation (81) gives $h_{cc}\simeq 1.7\times {10}^{-25}$. The response functions result in practically being unchanged; therefore, we gain an order of magnitude, i.e., $\delta {\tilde{T}}_1\approx 3.4\times {10}^{-24}$ m for $v_G=0.999$, $\delta {\tilde{T}}_1\approx 3.4\times {10}^{-23}$ m for $v_G=0.9$ and $\delta {\tilde{T}}_1\approx 1.6\times {10}^{-22}$ m for $v_G=0.1$.

Here, we discussed the variation of the photons’ paths in the z direction, which is the direction of the propagating relic GW. Clearly, analogous effects, which are due to the transverse effect of the GW (Equations (71) and (72)), are present in the x and y directions. Thus, Equations (74) and (50) can be used to discuss the general gravitational lensing in our model. We developed the complete computation in the z direction; the extension to the x and y directions is similar.

The global effect of these variations of the photons’ paths in the CMB should be analogous to the effect of water waves, which, in focusing light, create optical caustics, which are commonly seen on the bottom of swimming pools.

We stress that there are indications in the literature (see for instance [52]) that there is no amplification for $f\left(R\right)$ if compared with general relativity, while in this paper, we claim the amplification [18]. The key point here is the following. The ordinary transverse strain due to the scalar field in $f\left(R\right)$ theories is, in general, even lower with respect to the standard transverse strain in general relativity. On the other hand, due to the presence of the mass, in $f\left(R\right)$ theories, the third scalar polarization admits also a longitudinal strain. In this case, the correspondent longitudinal response function, i.e., Equation (80) in this paper, is frequency dependent. Thus, at high frequencies, the total signal can, in principle, be higher in $f\left(R\right)$ theories with respect to general relativity. This is also in agreement with the results in [7,15,16].

4. Chaos and Relativity in Orbital and Optical Systems

The consequences of GWs form $f\left(R\right)$ theories are observable fingerprints on the structure of the Universe. Massive GWs will act as lenses, which generate caustics in the motion of light and other particle fields. These caustics will then have measurable influences on the CMB or upon the distribution of galaxies in the Universe out to $z=1$ and beyond. The following looks at the issue of how general relativity can amplify chaotic dynamics and, further, can amplify optical chaos. This is illustrated in a three-body problem and in an elementary optical model. This digression into another aspect of relativity is meant as a way to set up analysis for the phenomenology of massive gravity waves. This illustrates how to proceed through the examination of elementary systems. The extension to more complex structures, such as a many-body problem of galaxies and dark matter, will require numerical methods.

One of the early tests of general relativity was that it predicted the perihelion precession in the orbit of Mercury [17]. This is a departure from Newtonian gravity that is largely post-Newtonian, or first order or to $O(1/c^2)$. These general relativistic corrections are completely integrable, and there is no chaotic dynamics associated with them. In a three-body problem, with a large central mass, a larger distant mass, which is treated as Newtonian, and a smaller satellite with $O(1/c^2)$ relativistic departures will exhibit chaotic dynamics in the small body. The additional relativistic corrections will interplay with the irregular chaotic dynamics and are shown below to contribute to a Lyapunov exponent [53]. In effect, a Lyapunov exponent $\lambda =log(\Lambda )$ will have a relativistic correction $\Lambda ={\Lambda }_0+\Lambda \left(O\left(c^{-2}\right)\right)$, and this correction then amplifies the chaotic behavior of the system. This is extended to optical systems. Einstein lenses [54] are a Newtonian gravitational phenomenon, and general relativistic corrections to $O(1/c^2)$ are minor, for the impact parameter on such a gravitating body is too small to be observationally significant. Yet, for a complex Einstein lens, say analogous to a compound lens, due to the smaller scale clumping of matter, a light ray may have a succession of small angular deviations. These angles of deviation will have a compounding effect similar to the angle deviations of a particle in an arena. This will result in increasingly complex optical caustics, which in analogue with chaos, are difficult to predict. This is further compounded if the gravitating clumps of matter are difficult to observe directly, such as with dark matter [55]. In a manner similar to the case with orbital dynamics, general relativistic corrections may also enhance this optical chaos or turbulence. This section connects two different aspects of chaos and relativity to present issues with the analysis of three-body systems with parameterized post-Newtonian parameters. Subtle enhancements of chaotic dynamics or the increase in a Lyapunov exponent might be documented in such a system. This should then be an observable characteristic of complex relativistic systems. The optical analogue illustrates how fine-detailed structure in a distribution of matter, which is an Einstein lens, could influence the complexity of optical caustics. Localized regions of large gravity fields could then further amplify this complexity, as well. This might lead to methods for mapping any local density variation in dark matter.

We stress that the numerical values of the Lyapunov exponent, λ, in general relativity are not gauge invariant, that is, they depend on the chosen coordinate system [56]. Therefore, for the same dynamical system, chaotic behavior may appear in some frames, but not in others [56]. Following [57], we find three different problems when one uses the Lyapunov exponent in general relativity:

• The reference systems have no unified time.
• The separation of space and time in the four-dimensional spacetime varies for different observers.
• Time and space coordinates works only for events and sometimes have no physical meaning.

Consequently, we could get different values of the Lyapunov exponent in different coordinate systems. The problem can be solved if one uses proper time and proper distances instead [57]. In that case, it is indeed possible to consider a particle, called the “observer”, moving along an orbit in the spacetime [57]. That particle can understand if its motion is or is not chaotic, observing if the proper distances from neighbor particles are increasing exponentially or not with its proper time [57]. Hence, the point in [56] that the Lyapunov exponent is not gauge-invariant in general relativity is correct. However, the point of this is to examine the possible role of general relativity in the amplification of chaotic dynamics. In effect, general relativity applies to a body close to the star or large mass, where these gravitationally interact by Newtonian gravity to a third body. The purpose is to illustrate how chaos in Newtonian mechanics may be amplified if the system interacts with a semi-relativistic system in a stronger gravity field. Within this approximation, the question concerning the invariance of the Lyapunov exponent in general relativity for the Newtonian dynamical body is a small effect. The Lyapunov exponent applies strictly to the Newtonian part of the problem.

4.1. General Relativity to $O(1/c^2)$

In general relativity, the equation of motion for a test mass particle around a fixed central mass is [17]:

$$\frac{d^{2}u}{d{\theta }^2}+u=\frac{GM}{l^2}+\frac{3GM{u}^2}{c^2}$$

(83)

Here, l is the constant specific angular momentum. We recognize this differential as the harmonic oscillator equation of Newtonian mechanics with a constant force, $GM/l^2$, plus the term, $\sim {(u/c)}^2$. The anomaly angle, θ, obeys the dynamical equation [17]:

$$\frac{d\theta }{ds}=\frac{l}{r^2}=l{u}^2$$

(84)

and,

$$\frac{dt}{ds}=\frac{E}{1-2GMu/c^2}$$

(85)

for E, the potential energy per unit mass of the particle “at infinity” $=constant$. For $GM/c^2\ll 1$, we may solve this problem by perturbation methods. The solution of interest is $O\left(1\right)$ plus $O\left(c^{-2}\right)$, which would be Newton plus first order GR correction. The expansion is carried out with the variables $u,\theta$ according to:

$$u=u_0+ϵu_1+O\left({ϵ}^2\right)$$

(86)

$$\theta ={\theta }_0+ϵ{\theta }_1+O\left({ϵ}^2\right)$$

(87)

Here, the term $ϵ=1/c^2$ gives the order of the expansion. The differential with respect to θ to the first order in ϵ is taken as:

$$\frac{d}{d\theta }\simeq \frac{d}{d{\theta }_0}+ϵ\frac{d}{d{\theta }_1}$$

(88)

If we input the expansion for u in Equation (86) into the differential equation of motion (Equation (83)) the following two equations are obtained:

$$O\left(1\right):\frac{d^2{u}_0}{d{\theta }_0^2}+u_0=\frac{\kappa }{l^2}$$

(89)

$$O\left(ϵ\right):\frac{d^2{u}_1}{d{\theta }_0^2}+u_1-3\kappa u_0^2=0$$

(90)

The term $d^2{u}_0/d{\theta }_{0}d{\theta }_1=0$, since $u_0$ is not a function of ${\theta }_1$. Further, the term $\kappa =GM$. The $O\left(1\right)$ differential Equation (89) has the solution:

$$u_0=\frac{\kappa }{l^2}(1+{ϵ}^{\prime }cos({\theta }_0+\alpha ))$$

(91)

which is the standard Newtonian solution for the radial velocity for a particle with orbital eccentricity ${ϵ}^{\prime }$ and anomaly angle α [17]. Now, let us consider on the expansion of θ. We set $E=1$ and insert this into the equation for the angular velocity equation:

$$\frac{d\theta }{dt}=(1-2\kappa ϵu)l{u}^2$$

(92)

where $1-2\kappa ϵu$ is the Schwarzschild transformation between proper and standard time coordinates [17]. This differential equation has the two contributing parts:

$$O\left(1\right):\frac{d{\theta }_0}{dt}=l{u}_0^2$$

(93)

$$O\left(ϵ\right):\frac{d{\theta }_1}{dt}=2l{u}_0-2\kappa l{u}_0^3$$

(94)

We are primarily concerned at this point in the solution to order $O\left(ϵ\right)$ for the orbit of a test mass in a GR orbit,

$$\frac{d^2{u}_1}{d{\theta }_0^2}+u_1-3\kappa u_0^2=0$$

(95)

where the Newtonian solution, $u_0$, is given by Equation (91). The square of $u_0$ in the non-homogenous term is:

$$u_0^2={\left(\frac{\kappa }{l^2}\right)}^2\left(1+2{ϵ}^{\prime }cos({\theta }_0+\alpha )+{{ϵ}^{\prime }}^{2}co{s}^2({\theta }_0+\alpha )\right)$$

(96)

which by elementary trigonometric identities is:

$$u_0^2={\left(\frac{\kappa }{l^2}\right)}^2\left((1-{ϵ}^{\prime })+2{ϵ}^{\prime }co{s}^2(({\theta }_0+\alpha )/2)+{{ϵ}^{\prime }}^{2}co{s}^2(\theta +\alpha )\right)$$

(97)

The reason for doing this is that the solution is elementary at this point. The first non-homogeneous term is going to give a solution:

$$\sim \frac{{\kappa }^3}{l^4}(1-{ϵ}^{\prime })(1+{ϵ}^{\prime }cos({\theta }_0+\alpha ))$$

(98)

and the quadratic trigonometric functions determine the solution:

$$u_1=\frac{{\kappa }^3}{l^4}\left(\left(1+{ϵ}^{\prime }cos({\theta }_0+\alpha )\right)+\frac{2{ϵ}^{\prime }}{3}\left(cos({\theta }_0+\alpha )-3\right)+\frac{{{ϵ}^{\prime }}^2}{2}\left(cos\left(2({\theta }_0+\alpha )\right)-3\right)\right)$$

(99)

The hard part is the perturbation of the third planet. The Jovian planet obeys a similar dynamical equation, but where $c\to \infty$ and Newtonian dynamics is recovered as:

$$\frac{d^{2}v}{d{\theta }^{\prime 2}}+v=\frac{\kappa }{L^2}$$

(100)

Here $v=1/r_2$ for this additional planet, and we define $u=u_0+ϵu_1=1/r_1$. Similarly, the angular momentum is defined by [17]:

$$\frac{d{\theta }^{\prime }}{dt}=\frac{L}{r_2^2}=L{v}^2$$

(101)

The angle, ${\theta }^{\prime }$, may exist in a different plane than θ, yet, as an approximation, we put both angles in the same plane of motion. Now, we need the coupling between the two bodies. We assume they are Newtonian as:

$$\overrightarrow{F}=GMm\frac{{\overrightarrow{r}}_1-{\overrightarrow{r}}_2}{|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2{|}^3}$$

(102)

which is approximately,

$$\overrightarrow{F}=\frac{GMm}{r_2^3}(1+\frac{3}{2}\frac{{\overrightarrow{r}}_1·{\overrightarrow{r}}_2)}{r_2^2}({\overrightarrow{r}}_1-{\overrightarrow{r}}_2)$$

(103)

To find the distance $|{\mathbf{r}}_1-{\mathbf{r}}_2|$, we consider the plane of the two orbits as complex valued and that the positions of the test mass and the larger mass are give by ${\mathbf{r}}_1=r_1{e}^{i{\theta }_1}$ and ${\mathbf{r}}_2=r_2{e}^{i{\theta }_2}$, and so, the distance between the two masses is given by:

$$|{\mathbf{r}}_1-{\mathbf{r}}_2{|}^2=r_1^2+r_2^2-2r_2{r}_{2}cos({\theta }_1+{\theta }_2)$$

(104)

The potential energy,

$$U(r_1,r_2)=-\frac{GMm}{|{\mathbf{r}}_1-{\mathbf{r}}_2|}$$

(105)

defines the force in Equation (102) by $\mathbf{F}=-\nabla U$. For $r_2\gg r_1$, the denominator in the potential may then be cast in the $u,v$ variables:

$$U(u_1,u_2)\simeq GMmv\left(1-{\left(\frac{v}{u}\right)}^2+2\frac{v}{u}cos\left(({\omega }_1+{\omega }_2)t\right)\right)$$

(106)

Here, ${\omega }_i=d{\theta }_i/dt$, for $i=1,2$ for the two bodies. This is the perturbing potential for the two orbits of the bodies in the same plane.

The total Hamiltonian is then,

$$\begin{array}{ccc}\hfill H& =\hfill & +ϵH_1^{ho}+H_v^{ho}+{\kappa }^{\prime }v\left\{(1-{\left(\frac{v}{u_0}\right)}^2+\left(\frac{v}{u_0}\right)cos(\theta +{\theta }^{\prime })\right\}\hfill \\ \hfill & \hfill & -3ϵu_0^2{u}_1-ϵ{\kappa }^{\prime }\left(1-\frac{v}{u_0}\right){\left(\frac{v}{u_0}\right)}^2{u}_1\hfill \end{array}$$

(107)

where the first three terms are harmonic oscillator Hamiltonians,

$$H_0^{ho}=\frac{1}{2}{p}_0^2+\frac{1}{2}\frac{\kappa }{l^2}{u}_0^2,H_1^{ho}=\frac{1}{2}{p}_1^2+\frac{1}{2}\frac{\kappa }{l^2}{u}_1^2,H_v^{ho}=\frac{1}{2}{p}_v^2+\frac{1}{2}\frac{{\kappa }^{\prime }}{L^2}{v}^2$$

(108)

We now have two order parameters $ϵ=1/c^2$ and another ${\kappa }^{\prime }=Gm$, where the mass, m, is the mass of the “Jovian” planet. The Hamiltonian term that scales according to $ϵ{\kappa }^{\prime }$ for $\theta ={\theta }_0+ϵ{\theta }_1$ is:

$$H_{ϵ\delta }\simeq -\kappa {}^{\prime }\left(1-\frac{v}{u_0}\right){\left(\frac{v}{u_0}\right)}^2{u}_1\left(1-\frac{v}{u_0}\right){\left(\frac{v}{u_0}\right)}^2{u}_1$$

(109)

To compute the Lyapunov exponent explicitly, the gradients of the Hamiltonian with $p_0,p_1,p_v$ and $u_0,u_1,v$ are first found. With $v/u_0\ll 1$, $u_1\ll u_0$, these are then to order ${(v/u_0)}^2$:

$${\nabla }_{p_0}H={\dot{u}}_0,{\nabla }_{p_1}H=ϵ{\dot{u}}_1,{\nabla }_{p_v}H={\dot{u}}_v$$

(110)

$${\nabla }_{u_0}H=\frac{\kappa }{l^2}{u}_0-6ϵu_0{u}_1-{\kappa }^{\prime }\frac{v^2}{u_0^2}cos(\theta +{\theta }^{\prime })$$

(111)

$${\nabla }_{u_1}H=\frac{ϵ\kappa }{l^2}{u}_1-3ϵu_0^2-ϵ{\kappa }^{\prime }\frac{v^2}{u_0^2}$$

(112)

$${\nabla }_{v}H=\frac{{\kappa }^{\prime }}{L^2}v+{\kappa }^{\prime }\left(1-2\frac{v^2}{u_0^2}+2\frac{v}{u_0}cos(\theta +{\theta }^{\prime })\right)$$

(113)

These are the forces $F=-\nabla H$, due to the three configuration variables, $u_0$, $u_1$ and v. The last right-hand side terms in ${\nabla }_{u_1}H$ are dependent upon both the general relativistic correction, $O(1/c^2)$, and the gravitational coupling with the Jovian planet, ${\kappa }^{\prime }$.

We consider the change in the phase space flow:

$$Z+\Delta Z=(u+\Delta u,p+\Delta p)$$

(114)

The change in momenta, due to the perturbation from the Jovian planet, is;

$$\Delta p\simeq \Delta t\left[{\kappa }^{\prime }\left(1-2\frac{v^2}{u_0^2}+2\frac{v}{u_0}cos(\theta +{\theta }^{\prime })\right)-{\kappa }^{\prime }\frac{v^2}{u_0^2}cos(\theta +{\theta }^{\prime })-ϵ{\kappa }^{\prime }\frac{v^2}{u_0^2}\right]$$

(115)

where the last term is a coupling of general relativistic $O(1/c^2)$ effects and planetary perturbation; to $O({\kappa }^{\prime }/c^2)$ $\Delta u\propto \Delta p$. Define $\Delta p\left(t\right)$ to be the deviation in momentum due to planetary perturbation, and let $\delta p\left(t\right)$ be the deviation due to the $O(1/c^2)$ coupling term. The Lyapunov exponent is then,

$$\lambda \simeq \underset{t\to \infty }{lim}\frac{1}{t}ln\left(\frac{\Delta p\left(t\right)+\delta p\left(t\right)}{\Delta p\left(t_0\right)}\right)\simeq \underset{t\to \infty }{lim}\frac{1}{t}\left[ln\left(\frac{\Delta p\left(t\right)}{\Delta p\left(t_0\right)}\right)+\frac{\Delta p\left(t_0\right)\delta p\left(t\right)}{\Delta p\left(t\right)}\right]$$

(116)

so that,

$$\lambda \simeq {\lambda }_0+ϵ{\left(\frac{v\left(t_0\right)}{u_0\left(t_0\right)}\right)}^2$$

(117)

with ${\lambda }_0$ defined for $ϵ=0$. The exponential divergence in phase space between nearby trajectories has then contribution in addition to ${\lambda }_0$ with $Z\left(t\right)\sim Z\left(t_0\right)e^{{\lambda }_{0}t}{e}^{ϵ{(v/u_0)}^{2}t}$. Thus, general relativity will bring about the onset of chaotic behavior, or the breakdown of numerical unpredictability, earlier.

This should then manifest itself in semi-relativistic systems with three bodies. A system, such as two neutron stars in a mutual relativistic orbit with a third companion further away and executing Newtonian dynamics, of this sort will then have more chaotic behavior, which is amplified by general relativistic effects. This simplified model suggests that a general parameterized post-Newtonian-multibody perturbative theory is needed. Such a model will then be more suited for the examination of complex general relativistic systems that include several bodies.

4.2. $O(1/c^2)$ Optical Corrections in Einstein Lensing

The Einstein lensing of light is now a common observational feature of deep space astronomy, since the launch and repair of the Hubble Space Telescope; see [58] and references within. Here, a complex optical gravitational lensing system is discussed with some analogues to the mechanics above. A large elliptical galaxy will have an overall gravitational lensing effect, but there may be sub-lensing, as well, if the density of dark matter that exists has some variation. This results in a type of optical turbulence, analogous to chaos. Further, this may also be amplified by general relativistic effects. Unknown configurations might exist with dark matter density increasing in the vicinity of a large black hole.

The general theory of gravitational lensing [32,49] shows that a light ray that approaches within a radius, $r\gg 2GM/c^2$, will be deflected approximately by an angle, $\theta =GM/r{c}^2$. In a more general setting, the deflection of light is given by the Einstein angular radius:

$${\theta }_E=\sqrt{\frac{4GM}{c^2}\frac{d_{ls}}{d_l{d}_s}}$$

(118)

where $d_{ls},d_l,d_s$ are the angular diameters to the gravitational lens, the source and the distance between the gravitational lens and the source; for $d_{ls},d_l,d_s$, the angular diameters to the gravitational lens, the source and the distance between the gravitational lens and the source. The condition $d_s=d_l+d_{sl}$ obtains locally where cosmological frame dragging is small. This theory is the weak gravitational lensing approximation, where the deflection of light is essentially a Newtonian result [49]. The distance relationships are determined by $\theta d_s=\beta d_s+{\alpha }^{\prime }{d}_{ls}$. The reduced angle of deflection $\alpha \left(\theta \right)=(d_{ls}/d_s){\alpha }^{\prime }\left(\theta \right)$ gives a relationship between the angles of importance $\alpha \left(\theta \right)+\beta =\theta$.

Complex distributions of matter can act similar to a compound lens in a weak gravitational limit. However, light rays that pass close to clumps of matter to exhibit $O(1/c^2)$ deviations will exhibit deviations from this linear summation. A light ray that passes through a set of random lenses will display caustics, which are similar to the caustics seen on the bottom of a swimming pool. The fine-grained structure in an Einstein lens can exhibit caustics, which occur due to nonlinear perturbation in the density profile of matter. This nonlinearity in the symmetry of the lens will produce caustics, which are analogous to chaos. The occurrence of a caustic has its connections with catastrophe theory [59] and the onset of a fold, which is also a mechanism for the bifurcation of vector fields in Hamiltonian chaos.

For the position of a source, $\overrightarrow{x}$, the propagation of light along the z-axis from this source then reduces the visual appearance of the object to $\overrightarrow{\xi }=({\xi }_x,{\xi }_y)$ along the axis of optical propagation. The weak gravitational lensing of light [49] then indicates that the deflection of the appearance of this object along the axis of optical propagation is given by:

$$\Delta \overrightarrow{\xi }=\nabla \Phi \left(\xi \right)$$

(119)

for ξ, the position of the image with the mass present, and $\Phi \left(\xi \right)$, the gravitational potential. The difference in the vector position of the image ${\overrightarrow{\xi }}_i-{\overrightarrow{\xi }}_s$ is the difference between the position with the mass present and without it being present. The potential term obeys the Poisson equation [13], so that:

$${\nabla }^2\Phi =2\frac{\Sigma \left(\overrightarrow{\xi }\right)}{{\Sigma }_c}$$

(120)

The integration over the direction of propagation then gives the mass density in the plane of the image, often called the surface mass density, $\Sigma \left(\overrightarrow{\xi }\right)$. The angle of deflection, α, is then determined by the Poisson equation and the potential as:

$${\overrightarrow{\alpha }}^{\prime }\left(\overrightarrow{\xi }\right)=\frac{4G}{c^2}\int \frac{(\overrightarrow{\xi }-{\overrightarrow{\xi }}^{\prime })\Sigma \left({\overrightarrow{\xi }}^{\prime }\right)}{|\overrightarrow{\xi }-{\overrightarrow{\xi }}^{\prime }{|}^2}{d}^2{\xi }^{\prime }$$

(121)

for $\Sigma \left(\overrightarrow{\xi }\right)$, a mass/area density distribution in the image. The function, $\Sigma \left(\overrightarrow{\xi }\right)$, plays the role of an index of refraction based upon a mass distribution, which, for a thin lens, will give the angle of deviation. For a gravitational thin lens, a weak field that is very small compared to the optical path length, and $\Sigma \left(\overrightarrow{\xi }\right)$ is a constant. The deflection angle is simply:

$$\alpha \left(\xi \right)=\frac{4\pi G}{c^2}\frac{\Sigma \left(\xi \right)d_{ls}\xi }{d_s}$$

(122)

where for small angles $|\overrightarrow{\xi }|=\xi =d_l\theta$ and,

$$\alpha \left(\xi \right)=\frac{4\pi G\Sigma }{c^2}\frac{d_{ls}{d}_l}{d_s}=\frac{\Sigma }{{\Sigma }_c}\theta$$

(123)

for the critical mass density ${\Sigma }_c=(c^2/4\pi G)(d_s/d_{ls}{d}_l)$. This is the minimal mass density that might be distributed in the area of an Einstein ring [60]. For a more complex arrangement of gravitational lenses, such as large density nonlinearities, the mass density, Σ, has a general form:

$$\Sigma \left(\xi \right)=\int dz\sigma (d_l{\xi }_x,d_l{\xi }_l,z)$$

(124)

The position of the image then plays the role of the vector, $\overrightarrow{r}$, in Equations (83) and (84) and beyond in the above discussion. Let the reciprocal of the vector norm $|\overrightarrow{\xi }|=1/\nu$ play the role of u. The analogue of the Newtonian equation of motion in Equation (1) with $c\to \infty$ is then,

$$\Delta \nu =\frac{\kappa }{{\left(c\xi \right)}^2}$$

(125)

for ξ, the impact parameter. Now, the Newtonian description of gravitational lens deflection has the effective photon angular momentum per mass term $j={\left(c\xi \right)}^{-1}$. The general relativistic extension of this equation is then,

$$\Delta \nu =\frac{\kappa }{{\left(c\xi \right)}^2}+\frac{3\kappa {\nu }^2}{c^2}$$

(126)

To order expansions, the analogue of Equations (89) and (90) are:

$$\Delta {\nu }_0=\frac{\kappa }{{\left(c\xi \right)}^2}$$

(127)

$$\Delta {\nu }_1=\frac{3\kappa \Delta {\nu }_0^2}{c^2}$$

(128)

for $u_0$, the reciprocal of ξ. The term $\kappa =GM$ for a general distribution in the plane of the Einstein ring is $\kappa /\xi \simeq 4\pi G\Sigma d_{ls}{d}_l/d_s$, which reproduces the Einstein ring case in the first approximation. The second order term is $\Delta \nu =\Delta \xi /{\xi }^2$ and $\alpha \simeq 2\Delta \xi /\xi$, which reproduces the weak field gravity lens result. The $O(1/c^2)$ correction gives an effective general relativistic correction term,

$$\Delta \alpha \simeq \frac{3\kappa d_l^2{\theta }^2}{c^2}{\left(\frac{\Sigma }{{\Sigma }_c}\right)}^2$$

(129)

This correction term is not likely to be detected directly by extra-galactic sources, such as the dark matter lensing of light by the Abell galaxy cluster [61].

4.3. Optical Chaos

Just as the $O(1/c^2)$ correction to Newtonian dynamics enhanced chaotic dynamics or contributed to a Lyapunov exponent, we might expect a similar amplification of optical chaos or the statistical appearance of caustics by the clumping of matter. Small local region where gravitating mass is clumped together will result in the deviation of the light ray by some small angle, $\delta \theta$, which is an error in computing the subsequent tracing of the ray. With a succession of n such small deviations, the first angle deviation is amplified by $\simeq 2^n\delta {\theta }_1$, the next by $\simeq 2^{n+1}\delta {\theta }_2$, where, for large n and ${\theta }_i\simeq \theta$$\forall i$, the total angular error in computing a ray trace will be approximately $2^{n+1}\delta \theta$. This is analogous to the arena problem of computing the trajectory of a ball.

The vector, $\overrightarrow{\xi }$, describes the visual appearance of a distant object along the axis of propagation. This vector describes the deformation of a wave front by the lensing action of the intervening gravitating body. The gravitational lens is usually considered as a symmetric lens [49], but nature may provide local clumping of material, which introduces some chaos in the ray tracing. Further, the overall gravitating lens may be sufficient enough to produce small $O(1/c^2)$ relativistic deviations from a purely Newtonian lensing. Above, the formula for this relativistic deviation is given. What is then needed is an analogue to the Lyapunov exponent for the classical unpredictability of a ray trace due to Newtonian gravitational sources. A multiple set of ray tracings is then a description of the deformation of an electromagnetic wave front and perturbations on the vector, $\overrightarrow{\xi }$. In what follows, such a development is presented to describe the chaotic perturbation of this vector.

The propagation of a plane electromagnetic wave front is given by $\psi \left(\overrightarrow{r}\right)={\psi }_0{e}^{i\overrightarrow{k}·\overrightarrow{r}-\omega t}$. The occurrence of a gravitational lens perturbs the the wave front according to:

$${\psi }^{\prime }\left(\overrightarrow{r}\right)=\left(\chi \left(\overrightarrow{r}\right)e^{i\varphi \left(\overrightarrow{r}\right)}\right)\psi \left(\overrightarrow{r}\right)$$

(130)

Here, the $\varphi \left(\overrightarrow{r}\right)$ is the change in the wave front phase and $\chi \left(\overrightarrow{r}\right)$ is the change in the wave front amplitude. The vectors describing the visual appearance of the image are $\overrightarrow{\xi }={\nabla }_{r_{\left|\right|}}\varphi \left(\overrightarrow{r}\right)$ for $r_{\left|\right|}$ coordinate directions along the wave front. This means that $\Delta \varphi \left(\xi \right)=\Phi \left(\xi \right)$, which is a Poisson equation [13]. The phase deviations are caused by an effective index of refraction [32,49] in the Newtonian limit, and the gravitational potential is the source in the Poisson equation. A Gaussian random distribution of sources results in the second order structure function:

$$D_{\varphi }\left(\overrightarrow{\rho }\right)=\langle |\varphi \left(\overrightarrow{r}\right)-\varphi (\overrightarrow{r}-\overrightarrow{\rho }){|}^2\rangle$$

(131)

The vector $\overrightarrow{\rho }=\overrightarrow{\xi }+\overrightarrow{z}$, where $\overrightarrow{\xi }=\overrightarrow{r}-\overrightarrow{\rho }$, so $D_{\varphi }\left(\overrightarrow{\rho }\right)$ is a phase variance between two different direction in the aperture plane.

The phase terms obey a Poisson equation, where some distribution of the sources is present. For optical perturbations compatible with the second order structure function, the gravitating perturbations are in a Gaussian distribution, where Gaussian distributions of perturbing sources means that the Equation (120) becomes:

$${\nabla }^2\Phi =2\frac{\Sigma \left(\overrightarrow{\xi }\right)}{{\Sigma }_c}+\frac{1}{4\pi }\frac{\mu }{{\left(\sqrt{2\pi }\sigma \right)}^3}\prod _i{e}^{-{\xi }_i^2/2{\sigma }^2}$$

(132)

where each ${\xi }_i\ll \xi$. Each one of these sources gives a solution:

$$\varphi =(\mu /4\pi )(1/r)Erf(r/\sqrt{2\sigma })$$

(133)

for the variable ${\xi }_i=r$, and the solution converges to a point course in the limit $\sigma \to 0$. Each of these perturbing changes on the aperture vector is due to a succession of matter clumps. A photon that passes close to each clump is modeled as having its angle deviated, and its path is then stochastically deviated away from a path given by Equation (132). The small angle of deviation for $\overrightarrow{\alpha }\left(\overrightarrow{\xi }\right)\to \overrightarrow{\alpha }\left(\overrightarrow{\xi }\right)+\delta \overrightarrow{\alpha }\left(\overrightarrow{\xi }\right)$ is determined by the Gaussian distribution as:

$$\delta \overrightarrow{\alpha }\left(\overrightarrow{\xi }\right)=\frac{4G}{c^2}\int \prod _i\frac{({\overrightarrow{\xi }}_i-{\overrightarrow{\xi }}_i^{\prime })\rho \left({\overrightarrow{\xi }}_i^{\prime }\right)}{|\overrightarrow{\xi }-{\overrightarrow{\xi }}_i^{\prime }{|}^2}{d}^2{\xi }_i$$

(134)

for $\rho \left({\overrightarrow{\xi }}_i\right)=(1/4\pi )(\mu /{\left(\sqrt{2\pi }\sigma \right)}^3)e^{-{\xi }_i^2/\left(2{\sigma }^2\right)}.$ For simplicity, the angle of deviation, $\delta \overrightarrow{\alpha }\left(\overrightarrow{\xi }\right)$, will now be treated as a scalar, and with $\xi \gg {\xi }_i$, the angle deviation is:

$$\delta \alpha \simeq \kappa \int \prod _i{\xi }^{-1}\rho \left({\xi }_i^{\prime }\right)d{\xi }_i^{\prime }$$

(135)

such that $\alpha \simeq \langle {\xi }^{-1}\rangle$. This is a partition function analogous to that in the Ising model [62], but here, instead of a set of spins that exist in space, there are stochastic angle changes in a ray trace of light. In this particular model, these stochastic angle changes are assumed to be, on average, the same.

This partition function can be demonstrated to be similar to the Ising model. For the variation in the stochastic variable $\delta {\xi }_j={\xi }_j-{\xi }_{j-1}$ in the exponent, the product of any two variations vanish $\delta {\xi }_j\delta {\xi }_j\simeq 0$, so that:

$${\xi }_{i-1}{\xi }_{j-1}+{\xi }_i{\xi }_j=2{\xi }_{i-1}{\xi }_j$$

(136)

for $i=j$, the sum of these stochastic variables is:

$$\frac{1}{2}\sum _{j=0}^{n-1}{\xi }_j{\xi }_j=\sum _{j=1}^{(n-1)/2}{\xi }_{j-1}{\xi }_j-\frac{1}{2}({\xi }_0^2+{\xi }_n^2)$$

(137)

This means there exist additional endpoint terms, which do not conform to the Ising type of construction. However, for a large enough n, this error should be minimal. The expectation is approximately:

$$\langle \xi \rangle \simeq \frac{1}{\sqrt{2\pi \sigma }}\int \prod _{j=1}^{(n-1)/2}d{\xi }_j{\xi }^{-1}exp\left(-{\xi }_{j-1}{\xi }_j\beta \right)$$

(138)

for $\beta =1/{\sigma }^2$. β is analogous to the Boltzmann factor, which is determined by the scale at which matter is lumped together. A correlation length scale is:

$${\lambda }^2\simeq 1/log\left(tanh\beta \right)$$

(139)

which for $\beta \ll 1$, or equivalently for large σ, is $\lambda \simeq \sigma$. This approximate formula is a ray trace path analogue of the Lyapunov exponent in time, which determines a length, λ, where the prediction of a ray trace breaks down. This also illustrates that this breakdown of ray tracing occurs on a scale comparable to the length scale of the perturbing. This loss of ray trace prediction is manifested in deformations of the angle deviation across the aperture distance or deviations in the symmetry of an Einstein ring.

The goal now is to determine if there are enhancements of optical chaos, analogous to optical turbulence in the Earth’s atmosphere, due to $O(1/c^2)$ corrections. To examine this, we consider the angle of deviation, due to Newtonian gravity, optical path length turbulence and relativistic corrections as a conjugate to action variables $J,J^{\prime },J"$, and the path is given in a classical setting by the action:

$$(H+H^{\prime }+H^{\prime \prime })dt=(Jd\alpha \left(\xi \right)+J^{\prime }d{\alpha }^{\prime }\left(\delta \xi \right)+J^{\prime \prime }d{\alpha }^{\prime \prime }\left({\xi }_1\right)$$

(140)

For the angular momentum variables $J\simeq J^{\prime }\simeq J^{\prime \prime }$ in this equation, then the action is entirely governed by the angle deviation, which for $d\alpha =(d\alpha /dt)dt$ expresses this as a principle of least time, just as in the case of planetary motion. The angle of deviation due to the clumpiness of matter is approximated as:

$$\delta {\alpha }_c\simeq \kappa exp\left({\xi }^{2}log\left(n\right)\beta \right)$$

(141)

for n regions of matter or dark matter clumping. The region where the light ray is the most distorted by gravitating bodies is of a distance $\sim n\sigma =\sqrt{\beta /2}$, which then gives an approximate relativistic $O(1/c^2)$ correction:

$$\delta {\alpha }_g\simeq \frac{3\kappa n{\sigma }^2({\theta }^2+2\theta \delta {\alpha }_c)}{c^2}{\left(\frac{\Sigma }{{\Sigma }_c}\right)}^2$$

(142)

where $\delta \theta \simeq \delta {\alpha }_c$. An approximate Lyapunov exponent is then:

$$\lambda \simeq \underset{n\to \infty }{lim}\frac{1}{n}log\left(1+\kappa \left(\frac{{\Sigma }_c}{\theta {\Sigma }_c}\right)\left[e^{{\xi }^{2}log\left(n\right)\beta }+3\kappa n{\sigma }^2({\theta }^2+2\theta e^{{\xi }^{2}log\left(n\right)\beta }){\left(\frac{\Sigma }{{\Sigma }_c}\right)}^2\right]\right)$$

(143)

Here, there is an amplification of the ray trace uncertainty, or chaos, by the introduction of the $O(1/c^2)$ term, as seen in the term, $2\theta e^{{\xi }^{2}log\left(n\right)\beta }{(\Sigma /{\Sigma }_c)}^2$. For $\delta {\alpha }_c\simeq \theta$, the contributions to the chaotic ray traced path from relativistic corrections and chaos are comparable and will contribute equally to the randomness of the caustic gravitational lens.

The difference in the perturbed aperture vectors $\Delta \delta \overrightarrow{\xi }=\delta {\overrightarrow{\xi }}_i-\delta {\overrightarrow{\xi }}_s{\nabla }_{\delta \xi }\Phi$ determines the magnification $M=d\left(\delta {\xi }_i\right)/d\left(\delta {\xi }_s\right)$. From Hamilton’s equations, this is generalized to:

$$\Delta \delta \overrightarrow{\xi }={\nabla }_{\xi }H\simeq 2\kappa \xi log\left(n\right)d_l^2{\alpha }_c\left(1+\frac{3n{d}_l^2\theta }{c^2}{\left(\frac{\Sigma }{{\sigma }_c}\right)}^2\right)$$

(144)

with the deviation magnification computed accordingly. For this, written according to the radius of curvature R of a surface for a ray curve along a line of sight, we have that:

$$\frac{2\pi A\sqrt{2\left|R\right|}}{{\Sigma }_c}\Theta \left(\delta {\overrightarrow{\xi }}_i\right)=\frac{3n{d}_l^2\theta }{c^2}{\left(\frac{\Sigma }{{\Sigma }_c}\right)}^2$$

(145)

The magnification for $\eta =2\pi \sqrt{2\left|R\right|}/\left({\Sigma }_c\right)$ is $\mathcal{M}=1+\eta \Theta +O\left({\eta }^2\right)$. The curvature, R, defines a tangent for the ray trace, which defines a caustic when the line of sight is along this tangent. The caustics along lines of site occur at swallowtail folds in the magnification map.

For the orbits in different planes, the above must be generalized some. Similarly, the angular components are given by the tangent vector parallel to $p_i$, and from there, we may find the vector to each dark matter clump in its plane of motion with coordinates $\{u,\theta \}$ and $\{v,{\theta }^{\prime }\}.$ Further, as the angular momentum vector is given by $\overrightarrow{l}={\overrightarrow{r}}_1\wedge {\overrightarrow{p}}_1/m$ (similar for $\overrightarrow{L}$), then $\overrightarrow{L}$ is rotated relative to $\overrightarrow{l}$ by the Euler angles, α, β and γ. The rotation matrix is then:

$$\left[\mathbf{R}\right]=\left[\begin{array}{c}cos\left(\gamma \right)\end{array}\right]\left[\begin{array}{c}1\end{array}\right]\left[\begin{array}{c}cos\left(\alpha \right)\end{array}\right]$$

(146)

With these, we may be able to put the problem in a general setting. This part is yet to be worked, and their may be resources to aid in this effort.

The formation of filaments and domain walls of galaxies is then proposed to occur by this mechanism. The massive gravity waves in the very early Universe, such as in the post-inflationary period, deviate the motion of relativistic particles in a manner similar to the optical focusing of light. These focal points of matter then set up their own gravity fields, which persist through the subsequent expansion of the Universe. A mesh of caustics with swallowtail cusps heuristically may be seen to produce a web of regions where mass-energy is concentrated. The distribution of dark matter may then be established by caustics of gravitons and gravity waves in the early Universe.

5. Conclusions and Remarks

In this paper, the production of massive relic coherent gravitons in in a particular class of $f\left(R\right)$ gravity, which arises from string theory and their possible imprint in the CMB, have been discussed. The key point is that in the very early Universe, these relic gravitons could have acted as slow gravity waves. They may have then acted to focus the geodesics of radiation and matter. Therefore, their imprint on the later evolution of the Universe could appear as filaments and the domain wall in the Universe today. In that case, the effect on the CMB should be analogous to the effect of water waves, which, in focusing light, create optical caustics, which are commonly seen on the bottom of swimming pools. This issue has been carefully analyzed by showing gravitational lensing by relic GWs, i.e., how relic massive GWs perturb the trajectories of CMB photons.

The consequence of the type of physics discussed has been outlined from the point of view of an amplification of what could be called optical chaos.

For the sake of completeness, we stress that multiple imaging by gravitational waves and the associated caustic structure have been studied by other authors in frameworks different with respect to the approach of this paper; see for example [63] and the references within.