f(R) Gravity, relic coherent gravitons and optical chaos

We discuss the production of massive relic coherent gravitons in a particular class of f(R) gravity which arises from string theory and their possible imprint in 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 domain wall in the Universe today. In that case, the effect on 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 GWs perturb the trajectories of 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.


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/M theory caused a lot of interest about higher order gravity in more than 4 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 (R) 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 which leads to the unification of gravitation with the other forces. One of the most important goals of modern physics is to obtain an 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 space-time, 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 µν , describes both the dynamical aspects of gravity and the space-time 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 space-time. This is an unequaled problem in the context of quantum field theories, because the other theories are founded on a fixed space-time 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 which 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 un-clustered, 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, untill 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 modification of 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 (R) 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 modify of the Lagrangian, with respect to the standard Einstein-Hilbert gravitational Lagrangian [13], through the addition of highorder terms in the curvature invariants (terms like R 2 , R αβ R αβ , R αβγδ R αβγδ , R R, R k R, in the sense of f (R) Theories [3,7]) and/or terms with scalar fields non-minimally coupled to geometry (terms like φ 2 R) in the sense of Scalar-Tensor Theories [7].
In the tapestry of f (R) 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 Contraction over indices gives the result The physical consequences of this extension to curvature are fairly remarkable. The Bel-Robinson tensor is a vacuum curvature ∇T = 0, and it predicts gravity waves (GWs).

Gravity waves in f (R) theories
In f (R) 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] gives a traceless solution in standard GR [17] h = 0.
The modified gravity results in a terms that acts as a mass, where the wave equation is [15,16] h = m 2h .
The decomposition of the solution gives the standard h ++ and h ×× 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 α ′n R 2n [2], for α the string parameter. The action to being L the Lagrangian for everything else compactified on a Dp-brane. The extremization of this gives This action derives a modified form of the Einstein field equation Now, let us consider a metric with the form (3), which is a classical expression. The quadratic term corresponds to quantum corrections on the order of the parameter α. We consider this correction as due to fields φ µ ν so that the quantum correction to the metric is where we can regard φ σ µ φ νσ = δh µν . These fields are physically a quantum correction to the classical gravitational radiation h µν . In general, these fields are quantized fields. In a string theory framework [1] we may define operators of the form which is a harmonic oscillator quantization condition compatible with a string theory interpretation [1]. The graviton fields are given by the level matching condition on a closed string m = 0, which for the first mode n = 1 are states of the form such that α µ 1 α ν −1 |0 = |ω µν constructs the elementary states. The connection terms are computed as where the field is treated as a vierbein. Now, let us compute the curvature as 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 where the last term is zero in a linearized approximation. The linearized approximation occurs for long wavelength gravitons. Assume the connection term ω µ νσ = ∂ ν φ µ ρ φ d σ is eigenvalued with a wave number k a so the curvature tensor is The second order term in the action is then The term αk 4 is a quartic term in mass, where the string coupling constant ∼ G N has naturalized units of area. This is an intrinsic curvature in spacetime. The string coupling constant is about α ′ ∼ 10 −60 cm 2 , which is a small number. This also guarantees the viability of the action (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 [18,19,20]. The strong variations of the gravitational field in the early universe amplify the zeropoint 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 [18,19,20]. 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 [19], which is tuned in a good way with the WMAP data on the CMB (in particular exponential inflation and spectral index ≈ 1) [21,22]. Inflationary models of the early Universe were analyzed in the early and middles 1980's [19]. 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 [19] and contain a mechanism which creates perturbations in all fields [18,20]. 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 [18,20].
Relic gravitons can be characterized by a dimensionless spectrum [18,20] where is the (actual) critical density energy, ρ c of the Universe, H 0 the actual value of the Hubble expansion rate and dρ gw the energy density of relic GWs in the frequency range f to f + df . In standard inflationary model the spectrum is flat over a wide range of frequencies, see [18,20] and figure 1. The more recent value for the flat part of the spectrum that arises from the WMAP data can be found in [20], [20] 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 Hz. The vertical axis is log 10 Ω gsw . The inflationary spectrum rises quickly at low frequencies (wave which 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 constrains on scalar perturbations. This means that at LIGO and LISA frequencies, Ω gw (f )h 2 100 < 9 * 10 −13 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 φ µ ν according to harmonic oscillator operators b, b † , as a simple model of a string. The fields are expanded as where E µ ν is a tetrad, which is discussed more below. The summation runs from {−∞, ∞}. The product φ c a φ νσ = δh µν is a harmonic oscillator operator The sum gives a delta function on k and k ′ and the first term is the Hamiltonian, which after the use of a commutator the RHS term is where the 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 [23].
The tetrad E µ ν is the amplitude of the field. This plays a role similar to the minimal electric field E = ω/V ǫ 0 in box normalization [24]. A plausible choice for tetrad is then where α ′ is the string parameter and ω the frequency. For α ′ ≪ 1/ω this is a small term.
The curvature in quantum modes is then which is O(α) in the scale parameter. In fluctuations of the curvature the metric is g ∼ δL/L, for δL > L p . The connection terms are of order Γ ∼ δL/L 2 and curvatures are R ∼ δL/L 3 . The wave vectors are k ∼ 1/L and the scaling parameter is α ′ ω ∼ δ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 The second order term is formed from the total contraction on the Riemann tensor R αßµν R αβ σρ , it is (26) This term is to O(α ′2 ) and contributes a term O(α ′3 ) to the Lagrangian.
Consider the operator matrix operation R µνσρ R µνσρ |m . The first term has the operator matrix elements where n |n n| is a completeness sum and the momentum values assumed in the states |m and |n . This contributes an energy-squared. A similar analysis for m|R µνσρ R µνσρ gives where z * = z = i((1/4)Π µνσρ (k)Π µνσρ ). Hence the action phase due to the action e iS contains a squeeze operator. The final operator term is more complicated. The operator terms are evaluated by commuting operators and this leads to the square of number opera- are then a product of terms which represent a squeezed state operator.The squeeze operator S(z) acts upon the displacement operator which effectively creates a modified displacement operator The action of the squeezed state operator on b is SbS † = b cosh(|z|) + b † sinh(|z|), which is a Bogoliubov transformed operator [26]. For a set of bosons, here linear gravitons, with the same state there is then n α n / √ n|n states with the operator acting on this n α n (a † ) n / √ n acting on the vacuum. This operator may be formed from the S(z)D(α)S † (z) for α small and |z| ≫ |α| with and where we may then define z * α/4 → α, and the Bogoliubov transformation of the operator b−b † constructs a displacement operator [27]. 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 [28] W (φ µ ν ) = exp(´iφ µν e µ dx ν ). In the path integral The infinitesimal shift in the field φ → φ+δφ adjusts Z[φ, W ] → Z[φ+δφ, W ] = W and the expansion is where the invariance of the expectation gives This formula is only well defined for a polynomial function. So we make the following approximation. The loop is considered to be very small and in that way we can approximate the Wilson loop with so that the functional derivative of W (φ) is for ǫ µν a unit area. The solution is then Now, let us consider the second order expansion which gives the result where by continuing the series this leads to 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(α). The Wilson loop is a form of the Stokes' law [29] and In vacuum the canonical h ++ and h ×× polarizations obey h ++ = h ×× = 0 [17]. The longitudinal modes due to R 2 terms obeys [15,16] where the mass is a topologically induced mass. The longitudinal h c = φ 2 then defines the equation [15,16,20] where a Lorenz gauge sets terms with φ = 0 [15,16,20]. This term plays a role similar to the Helmholtz potential in electromagnetism but in the case of f (R) theories it results an effective potential through the identifications [15,20] Φ → f ′ (R) and which give a Klein -Gordon equation for the effective Φ scalar field [15,20] The φ µν which physically contributes to the Wilson integral has a source term which is the topological mass.

Potential imprint in Cosmic Microwave Background
We recall that CMB is thermal radiation filling the observable Universe almost uniformly [21,22]. Precise measurements of 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 ∼ 2.7 K. The peak of the spectrum results in the microwave range frequency of 160.2 GHz [21,22]. At the present time, the best available data on CMB arise from the WMAP Mission [21,22], which used symmetric, rapid-multi-modulated scanning, rapid switching radiometers to minimize non-sky signal noise. WMAP works within the Solar System and to take into account weak potential effects on CMB by relic massive GWs we can use the weak field approximation (the linearized theory). In the linearized theory, the coordinate system in which the space-time is locally flat has to be used and the distance between any two points, say A and B, is given simply by the difference in their coordinates in the sense of Newtonian physics [17]. 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 space-time within the Solar System is locally flat with respect to the global distribution of CMB. Our goal is to understand how relic massive GWs perturb the trajectories of CMB photons between A and B. The global effect results a particular gravitational lensing [30] 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 [30,31] n ≡ 1 + 2|V |.
In the usual Geometrical Optics, the condition n > 1 implies that the light in a medium is slower than in vacuum [32]. Then, the effective speed of light in a gravitational field is expressed by [30,31,32] Thus, one can obtain the Shapiro delay [33] by integrating over the optical path between the source and the observer: The situation is analogous to the prism [32].

Gravitational lensing in the direction of the propagating gravity wave
For the sake of simplicity, 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 the observer A and the spatial axes are centered in 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] The connection between Newtonian theory and linearized gravity is well known [13] g 00 = 1 + 2V, where V is the Newtonian potential. Let us consider the interval for photons propagating along the z -axis The condition for a null trajectory (ds = 0) gives the coordinate velocity of the photons which to first order is well approximated by 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 presence of the graviton, can be defined: 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 space-time, i.e. in absence of GW, and t ′ is the delay time which corresponds to the unperturbed photon trajectory: (i.e. t is the time at which the photon arrives in the position T , so T − z = t − t ′ ).
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 which is the equation for geodesic deviation in this frame [17], and R i 0k0 is the linearized Riemann tensor [17].
On the other hand, with an opportune choice of the Lorentz gauge, the linearization process of f (R) theories which generates the third longitudinal mode h c = h c (t − v G z) enables a conformally flat line element [15,16,20] where v G represents the group velocity of the massive GW. In fact, the velocity of every standard massless tensorial modeh µν is the light speed c, but the dispersion law for the modes of h c is that of a massive field which is a wavepacket [15,16,20]. Also, the group-velocity of a wave-packet of h c centered in − → p is [15,16,20] which is exactly the velocity of a massive particle with mass m (see Eq. (44)) and momentum − → p . This group-velocity is function of both of the mass and frequency of the wave-packet [15,16,20] Even if the coordinates (52) are different from the coordinates (60), we recall that the linearized Riemann tensor is gauge invariant [17]. Hence, we can calculate it directly from Eq. (60). Following [16] it is: that, in the case eq. (60), begins [16] R α 0γ0 = the different elements are (only the non zero ones will be written) [16]: By putting these results in Eq. (64) one gets [16] which shows that the field is not transversal. In fact, Eq. (59) implies [16] 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 → 0 the longitudinal force vanishes.
Equivalently we can say that there is a gravitational potential [16,17]: which generates the tidal forces, and that the motion of the test mass is governed by the Newtonian equation [16,17] − → r = − ▽ V.
Thus, the variation of the proper distance between A and B from its unperturbed value T which is due by the presence of the massive GW h c is eq. (77) can be integrated in the frequency domain by using the Fourier translation and derivation theorems where is the longitudinal response function for relic gravitons.
In order to use eqs. (79) and (80) we recall that relic gravitons represent a stochastic background [18,20]. Hence, one has to use average quantities [18,20]. The well known equation for the characteristic amplitude [18], adapted for the third component of GWs can be used [20]: obtaining, for example at 100 HZ and taking into account the bound (20), Considering a graviton propagating with a speed of v G = 0.999 (ultrarelativistic case), if we insert these values in eqs. 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 by the transverse effect of the GW (eqs. (71) and (72)), are present in the x and y directions. Thus, eqs. (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 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.

Chaos and Relativity in Orbital and Optical Systems
The consequences of GWs form f (R) 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. This indicates that relic gravitons could have a larger effect on the early distribution of matter. One of the early tests of general relativity was the prediction of the perihelion precession in the orbit of Mercury [17]. This is a departure from Newtonian gravity that is largely post-Newtonian to O(1/c 2 ). These general relativistic corrections are completely integrable and with 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 [34]. In effect a Lyapunov exponent λ = log(Λ) will have a relativistic correction Λ = Λ 0 + Λ(O(c −2 )), and this correction then amplifies the chaotic behavior of the system. This is extended to optical systems. Einstein lenses [35] 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 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 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 [36]. 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. It may further lead to an understanding of how a stochastic distribution of massive gravitons can lead to the distribution of galaxies in filaments and walls.

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] 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 ∼ (u/c) 2 . The anomaly angle θ obeys the dynamical equation [17] dθ ds and for E the potential energy per unit mass of the particle "at infinity" = constant.
For GM/c 2 ≪ 1 we may solve this problem by perturbation methods. The solution of interest is O(1) plus O(c −2 ), which would be Newton plus first order GR correction. The expansion is carried out with the variables u, θ according to Here the term ǫ = 1/c 2 gives the order of the expansion. The differential with respect to θ to first order in ǫ is taken as If we input the expansion for u in equation (86) into the differential equation of motion (83) the following two equations are obtained: The term d 2 u 0 /dθ 0 dθ 1 = 0 since u 0 is not a function of θ 1 . Further, the term κ = GM . The O(1) differential equation (89) has the solution which is the standard Newtonian solution for the radial velocity for a particle with orbital eccentricity ǫ ′ 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 where 1−2κǫu is the Schwarzschild transformation between proper and standard time coordinates [17]. This differential equation has the two contributing parts: We are primarily concerned at this point in the solution to order O(ǫ) for the orbit of a test mass in a GR orbit, where the Newtonian solution u 0 is given by equation (91). The square of u 0 in the non-homogenous term is which by elementary trigonometric identities is The solution is elementary at this point. The first non-homogeneous term is going to give a solution and the quadratic trigonometric functions determine the solution: The hard part is the perturbation of the third planet. The Jovian planet obeys a similar dynamical equation, but where c → ∞ and Newtonian dynamics is recovered as 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] dθ The angle θ ′ 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 which is approximately To find the distance | r 1 − 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 r 1 = r 1 e iθ1 and r 2 = r 2 e iθ2 and so the distance between the two masses is given by The potential energy defines the force in equation (102) by F = −∇U . For r 2 ≫ r 1 the denominator in the potential may then be cast in the u, v variables Here ω i = dθ 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 where the first three terms are harmonic oscillator Hamiltonians We now have two order parameters ǫ = 1/c 2 , and another κ ′ = Gm, where the mass m is the mass of the "Jovian" planet. The Hamiltonian term that scales according to ǫκ ′ for θ = θ 0 + ǫθ 1 is To compute the Lypunov 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 ≪ 1, u 1 ≪ u 0 these are then to order (v/u 0 ) 2 These are the forces F = −∇H due to the three configuration variables u 0 , u 1 and v. The last right hand side terms in ∇ u1 H are dependent upon both the general relativistic correction O(1/c 2 ) and the gravitational coupling with the Jovian planet κ ′ . We consider the change in the phase space flow The change in momenta due to the perturbation from the Jovian planet is where the last term is a coupling of general relativistic O(1/c 2 ) effects and planetary perturbation. To O(κ ′ /c 2 ) ∆u ∝ ∆p. Define ∆p(t) to be the deviation in momentum due to planetary perturbation, and let δp(t) be the deviation due to the O(1/c 2 ) coupling term. The Lyapunov exponent is then so that with λ 0 defined for ǫ = 0. The exponential divergence in phase space between nearby trajectories has then contribution in addition to λ 0 with Z(t) ∼ Z(t 0 )e λ0t e ǫ(v/u0) 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.

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 [37] 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 sublensing as well if the density of dark matter 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 [30,31] shows that a light ray which approaches within a radius r ≫ 2GM/c 2 will be deflected approximately by an angle θ = GM/rc 2 . In a more general setting the deflection of light is given by the Einstein angular radius 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 [30]. The distance relationships are determined by θd s = βd s + α ′ d ls . The reduced angle of deflection α(θ) = (d ls /d s )α ′ (θ) gives a relationship between the angles of importance α(θ) + β = θ. Complex distributions of matter can act similar to a compound lens in a weak gravitational limit. However, light rays which pass close to clumps of matter to exhibit O(1/c 2 ) deviations will exhibit deviations from this linear summation. A light ray which passes through a set of random lenses will display caustics which are similar to the caustics seen on the bottom of a swimming pool. 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 [38] 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 x, the propagation of light along the z axis from this source then reduces the visual appearance of the object to ξ = (ξ x , ξ y ) along the axis of optical propagation. The weak gravitational lensing of light [30] then indicates that the deflection of the appearance of this object along the axis of optical propagation is given by for ξ the position of the image with the mass present and Φ(ξ) the gravitational potential. The difference in the vector position of the image ξ i − ξ 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 The integration over the direction of propagation then gives the mass density in the plane of the image, often called the surface mass density Σ( ξ). The angle of deflection α is then determined by the Poisson equation and the potential as for Σ( ξ) a mass/area density distribution in the image. The function Σ( ξ) 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 Σ( ξ) is a constant. The deflection angle is simply where for small angles | ξ| = ξ = d l θ and for the critical mass density Σ c = (c 2 /4πG)(d s /d ls d l ). This is the minimal mass density which might be distributed in the area of an Einstein ring [39].
For a more complex arrangements of gravitational lenses, such as large density nonlinearities, the mass density Σ has a general form The position of the image then plays the role of the vector r in equations (83)-(84) and beyond in the above discussion. Let the reciprocal of the vector norm | ξ| = 1/ν plays the role of u. The analogue of the Newtonian equation of motion in equation 1 with c → ∞ is then for ξ the impact parameter. Now the Newtonian description of gravitational lens deflection has the effective photon angular momentum per mass term j = (cξ) −1 .
The general relativistic extension of this equation is then To order expansions the analogue of equations (89) and (90) are for u 0 the reciprocal of ξ. The term κ = GM for a general distribution in the plane of the Einstein ring is κ/ξ ≃ 4πGΣd ls d l /d s , which reproduces the Einstein ring case in the first approximation. The second order term is ∆ν = ∆ξ/ξ 2 and α ≃ 2∆ξ/ξ, which reproduces the weak field gravity lens result. The O(1/c 2 ) correction gives an effective general relativistic correction term 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 [40].

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 δθ, 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 ≃ 2 n δθ 1 , the next by ≃ 2 n+1 δθ 2 , where for large n and θ i ≃ θ ∀i the total angular error in computing a ray trace will be approximately 2 n+1 δθ. This is analogous to the arena problem of computing the trajectory of a ball.
The vector ξ 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 [30], but nature may provide local clumping of material which introduce some chaos in the ray tracing. Further, the over all 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 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 ξ. 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 ψ( r) = ψ 0 e i k· r−ωt . The occurrence of a gravitational lens perturbs the the wave front according to ψ ′ ( r) = χ( r)e iφ( r) ψ( r).
Here the φ( r) is the change in the wave front phase and χ( r) is the change in the wave front amplitude. The vectors describing the visual appearance of the image are ξ = ∇ r || φ( r) for r || coordinate directions along the wave front. This means that ∆φ(ξ) = Φ(ξ), which is a Poisson equation [13]. The phase deviations are caused by an effective index of refraction [30,31] 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 The vector ρ = ξ + z, where ξ = r − ρ, so D φ ( ρ) is a phase variance between two different direction in the aperture plane. The phase terms obey a Poisson equation, where some distribution of 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 where each ξ i ≪ ξ . Each one of these sources gives a solution for the variable ξ i = r, and the solution converges to a point course in the limit σ → 0. Each of these perturbing changes on the aperture vector is due to a succession of matter clumps. A photon which 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 α( ξ) → α( ξ) + δ α( ξ) is determined by the Gaussian distribution as . For simplicity the angle of deviation δ α( ξ) will now be treated as a scalar and with ξ ≫ ξ i the angle deviation is such that α ≃ ξ −1 . This is a partition function analogous to that in the Ising model [41], 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 δξ j = ξ j − ξ j−1 in the exponent, the product of any two variations vanish δξ j δξ j ≃ 0, so that for i = j the sum of these stochastic variables is 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 for β = 1/σ 2 . β is analogous to the Boltzmann factor, which is determined by the scale at which matter is lumped together. A correlation length scale is which for β ≪ 1, or equivalently for large σ is λ ≃ σ. 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 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 conjugate to action variables J, J ′ , J", and the path is given in a classical setting by the action (H + H ′ + H ′′ )dt = (Jdα(ξ) + J ′ dα ′ (δξ) + J ′′ dα ′′ (ξ 1 ).
For the angular momentum variables J ≃ J ′ ≃ J ′′ in this equation then the action is entirely governed by the angle deviation, which for dα = (dα/dt)dt expresses this as a principle of least time. Just as in the case of planetary motion. The angle of deviation due to clumpiness of matter is approximated as 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 The magnification for η = 2π 2|R|/Σ c is M = 1 + ηΘ + O(η 2 ). 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, θ} and {v, θ ′ }. Further, as the angular momentum vector is given by l = r 1 ∧ p 1 /m (similar for L), then L is rotated relative to l by the Euler angles α, β and γ. The rotation matrix is then [R] = cos(γ) 1 cos(α) .
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.

Conclusion remarks
In this paper the production of massive relic coherent gravitons in in a particular class of f (R) gravity which arises from string theory and their possible imprint in 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 domain wall in the Universe today. In that case, the effect on 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.