text:
# Time resolved excitonic coherent polarization dynamics in GaAs multiple quantum wells

## Abstract

We investigate pump induced modulation in the temporal evolution of the excitonic coherent polarization generated by a probe pulse in GaAs multiple quantum wells. For this, we perform fs time resolved pump-probe reflectivity measurements using frequency upconversion at both positive and negative pump-probe delays. A contribution arising from interference of excitonic linear and nonlinear polarization, which is usually not considered in describing pump-probe experiments, is required to understand the results.

Linear and nonlinear optical response of excitonic resonances in GaAs like semiconductor multiple quantum wells (MQW) to excitation by ultrashort laser pulses and the behavior of the primary and secondary radiation emitted by the excited macroscopic polarization has been of much interest in recent years. Measurement of time resolved emission in the reflection weber ; haacke ; hayes  and transmission kim ; lyngnes  directions as a result of resonant excitation of excitons in quantum wells (QW) by an ultrashort laser pulse showed a time-profile significantly different from that of the incident pulse as a result of coherent emission from the excitonic polarization surviving up to several ps. A modification in the temporal evolution of the coherent emission in the probe reflection and transmission direction occurs in pump-probe differential (PPD) measurements. This is caused by linear and nonlinear interactions among the populations and polarizations, excited in the MQW by the probe and the pump pulses, and the propagating electric fields. An interesting manifestation of the nonlinear optical effects in the pump-probe differential transmission fluegel ; joffre ; sokoloff ; neukirch  (PPDT) and reflection guenther ; bipul  (PPDR) is the transient oscillations observed in spectrally resolved (SR) measurements for negative ($``$ve) pump-probe delay ($`\tau `$). (Here, $``$ve delay corresponds to the situation when the probe pulse precedes the pump pulse.) As the PPD measurements entail detection of intensity either of the spectrally or time resolved signals, and as these intensities are not simply Fourier related, important phase related effects are not the same in the two measurements. Therefore, additional insights into the modification of the emission in the probe reflection by a pump pulse may be obtained from time domain measurements. However, to our knowledge, such measurements have not been reported so far.

The purpose of this paper is two-fold. First, we present our experimental investigations of the time resolved PPDR from excitons in GaAs MQWs performed using the frequency upconversion technique. We find that the modulation in the dynamics of probe induced excitonic coherence by a pump pulse at low intensities occurs with a rise time of a few ps. However, the rise becomes faster as the pump intensity increases. No signal having the pump pulse polarization is emitted along the probe reflection direction for pump-probe cross-linear polarization. A PPDR signal is observed for $``$ve delay in both time resolved (TR) and time integrated (TI) measurements. The other main objective of this paper is to point out that the theoretical expressions usually employed to explain the PPD measurements srink ; shah  are inadequate in explaining the TR- and TI-PPDR measurements. Traditionally, the PPD signal is expressed as an interference of the probe electric field with the excitonic nonlinear polarization. This was able to explain the SR-PPDT results lindberg ; sokoloff ; koch  in the regime of the third order nonlinear optical susceptibility. One prediction of this description is that no signal in TR-PPD measurements can occur for $`\tau <0`$. This theory allows nonzero signal only for $`\tau >0`$ and that too at the time of arrival of the probe pulse in a narrow time domain limited by the probe pulse width. Furthermore, the model also predicts that the TI-PPD (or equivalently, spectrally integrated (SI) PPD) signal at $``$ve delay should vanish. However, all these predictions are contradictory to our experimental results. To understand the origin and nature of the signals observed in TR- and TI-PPDR measurements, we reexamine the theory for calculating the PPD signals. We show that an additional nonlinear optical contribution arising from interference of excitonic linear and nonlinear polarization, which is usually ignored in the literature, has to be invoked to explain our results.

Our experiments are performed at 8 K on 17.5 nm thick GaAs multiple-QWs, 20 in number, separated by 15 nm Al0.33Ga0.67As barriers, with a 330 nm AlGaAs layer on top. Continuous wave photoluminescence spectra measured at 8 K reveal heavy hole (hh) exciton emission at 1.5305 eV, with a spectral full width at half maximum (FWHM) 0.8 meV. Ti-Sapphire laser pulses of width of 180 fs (spectral FWHM of 10 meV) and a repetition rate of 82 MHz are used to resonantly excite the hh excitons. However, as the hh and light hole (lh) exciton energies are separated by 5.7 meV, coherent excitation of lh excitons also occurs. Dependence of TI- and TR-PPDR signal on the pump-probe delay is measured using the experimental arrangement shown schematically in Fig. 1(a). The TR-PPDR signal as well as the time profile of the reflected probe pulse in absence of a pump pulse are measured using type-I frequency upconversion in a LiIO3 crystal on which the reflected probe pulse and a reference pulse are cofocussed. The pump, probe and reference pulses are all derived from the same laser beam. Figure 1(b) schematically shows the pulse propagation directions for the incident and reflected probe beam in the MQWs.

Delay dependence of the magnitude of the time (and spectrally) integrated PPDR signal (which is actually $``$ve) is plotted in Fig. 1(c). Note that the TI-PPDR signal does not vanish at $``$ve delay up to about $`6`$ ps. The signal decay for $`\tau <0`$ is nearly exponential with a time constant of about 2 ps. For large positive ($`+`$ve) delay, the signal decay is slow and is presumably governed by the exciton life time (full data is not shown). At small $`+`$ve delay, an additional fast decaying component is seen. The oscillations seen in the data at $``$ve and small $`+`$ve delay are the signature of quantum beats due to coexcitation of hh and lh excitons.

The TR-PPDR measurements are shown in Figs. 2 and 3 for both $``$ve and $`+`$ve $`\tau `$. Here, we take the probe pulse to arrive at $`t=0`$ and the pump pulse at $`t=\tau `$, so that the pump pulse arrives after (before) the probe pulse when $`\tau <0`$ ($`\tau >0`$). The average pump (probe) intensity $`I_1`$ ($`I_2`$) is 2 (0.2) mW in Fig. 2. For the data in Fig. 3, $`I_1`$ is varied, keeping $`I_2`$ fixed at 0.3 mW. The data shows quantum beat oscillations (Figs. 2 and 3) with a period of about 0.75 ps, corresponding closely with the hh-lh exciton energy separation of 5.7 meV. Also, shown in Fig. 2(a) (dash-dotted line) is the cross-correlation trace of the laser pulse, giving the time resolution of the upconversion measurement of about 250 fs. Note that the TR-PPDR signal in Figs. 2 and 3 is $``$ve at all $`t`$. This means that the coherent polarization induced by the probe pulse in the MQWs is reduced due to the influence of the pump pulse for both $``$ve and $`+`$ve delay. We show later that such a reduction may be explained by invoking effects like phase space filling and excitation induced dephasing within the framework of optical Bloch equations. (However, local field or exciton energy shift have no role in this.)

Consider first the case of $`\tau <0`$. The emission along the probe reflection direction following excitation by the probe pulse may survive so long as the exciton polarization radiates coherently. A pump pulse incident during this process modifies the subsequent emission. Figure 2(a) shows the time evolution of the PPDR signal measured at three $``$ve delays. As expected from causality, the onset of the TR-PPDR signal correlates with the arrival of the pump pulse, delayed with respect to the probe pulse. No signal is observed when $`\tau `$ exceeds about 8 ps as the probe pulse induced polarization dephases within this time. At low intensities, the signal shows a rise time larger than the laser pulse width. After reaching the peak, it decays approximately exponentially with a time constant of about 2.5 ps. Figure 3(a) plots the TR-PPDR signal measured at three values of average pump intensity $`I_1`$ at $`\tau =1`$ ps. While the decay rate of the signal appears to be rather insensitive to $`I_1`$, a reduction in rise time with increased $`I_1`$ is clearly seen.

The observation of a nonzero TR-PPDR signal even when the pump pulse is incident on the sample well after the probe pulse has passed may be compared with a similar situation noted previously in the case of four wave mixing (FWM) measurements leo ; wegener . However, within the regime of the third order nonlinear optical susceptibility, the nonlinear optical polarization relevant to the PPD is not the same as that in FWM. For example, transient oscillations are seen in SR-PPD measurements at $``$ve delays, but not in SR-FWM. The exact origin of the signals at $``$ve $`\tau `$ is different in the two cases. In the FWM case, the polarization induced by the pulse along $`𝐤_2`$ gets diffracted along $`2𝐤_2𝐤_1`$ after the other excitation pulse arrives along $`𝐤_1`$. This is a result of many body Coulomb interactions wegener ; wang ; mayer  (such as local field, biexcitons, excitation induced dephasing, for example). On the other hand, the TR-PPDR signal at $`\tau <0`$ is in principle possible within the third order nonlinear optical regime even without introducing the many body Coulomb interactions, as shown presently.

It is sometime stated that the SR-PPD signal observed for $`\tau <0`$ is a result of pump pulse being scattered in the probe reflection/transmission direction from a grating created by the probe induced exciton polarization and the pump pulse sokoloff ; koch ; shah . If a transfer of energy from the pump beam to the probe reflection direction occurs, this may give rise to a $`+`$ve TR-PPDR signal, contrary to the data of Figs. 2(a) and 3(a). We further investigate this point by performing the TR-PPDR measurements for the pump polarization parallel and perpendicular to the probe (and the reference pulse) polarization. We find that the TR-PPDR signal which is measured using type-I upconversion is identical for the two cases of pump polarization. This confirms that the emission in the probe reflection direction always has the polarization of the probe pulse sieh .

Figure 2(b) displays the modification caused by a pump pulse in the time evolution of radiation from the excitonic polarization induced in the MQWs by a probe pulse incident after a $`+`$ve delay ($`\tau >0`$). Unlike $`\tau <0`$, the TR-PPDR signal in this case has a large contribution at $`t=0`$, and then a tail for $`t>0`$. The signal at $`t=0`$ is essentially caused by the interference of the probe electric field with the excitonic nonlinear polarization induced by the pump and probe pulses. The signal in the tail part shows a rise and then a decay with $`t`$. The rise for $`\tau >0`$ \[Fig. 2(b)\] is slower compared to the $`\tau <0`$ case seen in Fig. 2(a). In fact, the overall shape of the $`t`$ dependence is different in the two cases \[Figs. 2(a) and (b)\]. Also, the TR-PPDR signal shows only a small dependence on the value of $`\tau `$ when $`\tau `$ is $`+`$ve. (We verify this by varying $`\tau `$ from 1 to 30 ps.) The pump pulse excites its own coherent exciton polarization which turns into incoherent excitons after several ps. The rather small sensitivity of the TR-PPDR signal to change in $`\tau `$ ($`>0`$) exhibited in Fig. 2(b) indicates that the modulation of emission in the probe reflection direction by coherent and incoherent excitons induced by the pump pulse is not strongly distinguishable. In Fig. 3(b), the TR-PPDR signal at small $`I_1`$ shows an initial rise before decaying for larger $`t`$. As $`I_1`$ is increased, the rising part is not seen, the signal decaying continuously with $`t`$. If the decay of the coherent emission in the reflected probe pulse direction were exponential in $`t`$, and if the pump modulation only enhances the exponential decay rate, the TR-PPDR signal should be given simply as a difference of two exponentially decaying terms which may show an initial rise and then a decay. We find that the data of Fig. 2(a) ($`\tau <0`$) can be fitted by the two-exponential form. However, it is difficult to fit the $`+`$ve delay data \[Fig. 2(b)\] this way.

To compare experiments with theory, coupled Maxwell-Bloch equations have to be solved stroucken ; jahnke  with appropriate boundary conditions and many body Coulomb effects to obtain the PPDR signal for layered semiconductor heterostructures as used in our experiments. However, the following simplified phenomenological approach may be useful in understanding the origin of the PPDR signal as a starting step. Here, the MQW structure is replaced by an effective single layer of exciton resonances. Neglecting the effects of multiple reflections within the barrier and well regions and within the cap layer and taking the same background refractive index for the barriers and the QWs as a first approximation, the reflectivity $`R`$ for the model sample structure (Fig. 1) can be written as

$$R\left|r_0_2+_{QW}e^{i\varphi }\right|^2,$$
(1)
where only terms linear in the reflection coefficient $`r_0`$ ($`<1`$) of the cap layer on top of the QWs are retained. Here $`_2`$ is the amplitude of the probe electric field incident on the MQW and $`_{QW}`$ is that of the field emitted from the MQW layer. The phase factor $`\varphi =2dn\mathrm{\Omega }/c`$ arises from propagation delay for the probe pulse in the region between the vacuum-sample and top layer-MQW interfaces, $`n`$ is the refractive index of the cap layer, $`c`$ is the speed of light in vacuum and $`\mathrm{\Omega }`$ is the circular frequency of the incident light. The PPDR signal $`\mathrm{\Delta }R`$ is defined as $`R_{\text{on}}R_{\text{off}}`$ where $`R_{\text{on}}`$ ($`R_{\text{off}}`$) is obtained from Eq. (1) when the QWs are excited (not excited) by the pump pulse. In calculating $`R_{\text{on}}`$ we assume that the pump pulse causes modulation $`\mathrm{\Delta }_{QW}`$ in $`_{QW}`$ but $`r_0`$ remains unaffected.

For a thin sample with negligible depletion of incident pulse intensities, and using slowly varying envelope approximation in time, it can be shown bloembergen ; stroucken  that $`_{QW}`$ in absence of pump pulses is proportional to $`i𝒫_2^{(1)}`$ and the modulation $`\mathrm{\Delta }_{QW}`$ caused by the pump pulse is proportional to $`i𝒫^{NL}`$, where $`𝒫_2^{(1)}`$ and $`𝒫^{NL}`$ are the amplitudes of the linear polarization induced by the probe pulse and nonlinear polarization induced by the pump and probe pulses, respectively. The lowest order term in $`𝒫^{NL}`$ is the third order polarization $`𝒫^{(3)}`$ which has products of the three amplitudes $`_2`$, $`_1`$ and $`_1^{}`$ in the present case ($`_1`$ is the amplitude of the pump electric field incident on the MQW).

Keeping only the lowest order term in $`𝒫^{NL}`$, we get the TR-PPDR signal as

$$\begin{array}{cc}\hfill \mathrm{\Delta }R(t,\tau )& \mathrm{Im}[_2^{}(t)𝒫^{(3)}(t,\tau )e^{i\varphi }]\hfill \\ & +C\mathrm{Re}[𝒫_2^{(1)}(t)𝒫^{(3)}(t,\tau )],\hfill \end{array}$$
(2)
where $`t`$ is the real time and $`C=2\pi \mathrm{\Omega }n/cr_0`$. Similarly, the SR-PPDR signal is given by

$$\begin{array}{cc}\hfill \mathrm{\Delta }R(\omega ,\tau )& \mathrm{Im}[_2^{}(\omega )𝒫^{(3)}(\omega ,\tau )e^{i\varphi }]\hfill \\ & +C\mathrm{Re}[𝒫_2^{(1)}(\omega )𝒫^{(3)}(\omega ,\tau )],\hfill \end{array}$$
(3)
where $`_2(\omega )`$, $`𝒫_2^{(1)}(\omega )`$, and $`𝒫^{(3)}(\omega )`$ are the Fourier transforms of $`_2(t)`$, $`𝒫_2^{(1)}(t)`$, and $`𝒫^{(3)}(t)`$ respectively. Eqs. (2) and (3) apply also for the PPDT case when the phase term $`\varphi `$ is set to zero and $`r_0`$ in $`C`$ is replaced by unity. The literature reports which are mostly concerned with SR-PPD measurements usually ignore sokoloff ; lindberg ; koch  the second term in Eq. (3). In our investigation of the SR-PPDR case bipul , it was found that the second term is indeed smaller than the first but is not negligible. We wish to emphasize here that unlike the SR-PPD case, the second term is very crucial for the TR-PPDR (and TR-PPDT) measurements. In fact, within the above model, it is easy to show that the TR-PPDR measurements can not be explained without the second term in Eq. (2).

Before evaluating $`\mathrm{\Delta }R(t,\tau )`$, we may examine Eq. (2) to understand qualitatively how the TR-PPDR signal arises. For this, we refer to the schematic pictures of Fig. 4. Let us assume that the probe and pump pulses are $`\delta `$-like in time. Consider the $``$ve delay case. The probe pulse is taken to arrive on the sample at $`t=0`$ and the pump pulse later at $`t=|\tau |`$. $`𝒫^{(3)}(t,\tau )`$ is nonzero (and decaying with $`t`$) only for $`t|\tau |`$. In the first term of Eq. (2) $`_2(t)`$ is a $`\delta `$-like function at $`t=0`$. There is no overlap between $`_2(t)`$ and $`𝒫^{(3)}(t,\tau )`$ and hence no TR-PPDR signal is expected \[Fig. 4 (upper panel)\] from the first term. However in the second term, $`𝒫_2^{(1)}(t)`$ is nonzero (and decaying with $`t`$) for $`t0`$. If $`|\tau |`$ is less than the exciton dephasing time, then there is an overlap in $`t`$ between $`𝒫_2^{(1)}(t)`$ and $`𝒫^{(3)}(t,\tau )`$ for $`t|\tau |`$ and we get nonzero TR-PPDR signal \[Fig. 4 (upper panel)\] due to the second term in Eq. (2). In the case of $`+`$ve delay, the pump pulse arrives on the sample at $`t=|\tau |`$ and the probe pulse at $`t=0`$. $`𝒫^{(3)}(t,\tau )`$ is nonzero only for $`t0`$ and decays with $`t`$. There is an overlap between $`_2(t)`$ and $`𝒫^{(3)}(t,\tau )`$ at $`t=0`$ and between $`𝒫_2^{(1)}(t)`$ and $`𝒫^{(3)}(t,\tau )`$ for $`t0`$. Hence for $`+`$ve delay, the first term contributes only at $`t=0`$ and the second term contributes for $`t0`$ \[Fig 4 (lower panel)\]. In view of the above, observation of nonzero TR-PPDR signal (and in fact TR-PPDT signal as well) at $`t>0`$ for both $`+`$ve and $``$ve $`\tau `$ in Figs. 2 and 3 can not be explained if only the first term in Eq. (2) is considered.

Furthermore, it is obvious from the $`t`$ ordering of $`_2(t)`$ and $`𝒫^{(3)}(t,\tau )`$ that the time integral of the first term in Eq. (2) vanishes for $`\tau <0`$, giving no TI-PPDR signal. Thus, once again, the observed nonzero TI-PPDR signal for $``$ve $`\tau `$ in Fig. 1(c) can arise only from the $`t`$ integration of the second term in Eq. (2) (which does not vanish). Equivalence of TI- and SI-PPDR signal can be established by using Parseval’s theorem for Eqs. (2) and (3). This immediately shows that the first term in Eq. (3) will lead to zero SI-PPDR signal, though it contributes to SR-PPDR measurements.

We now assess how far a calculation of $`𝒫_2^{(1)}(t)`$ and $`𝒫^{(3)}(t,\tau )`$ in Eq. (2) based on a perturbative solution of the coupled optical Bloch equations for a non-interacting two-level assembly up to third order in electric fields can reproduce the features shown by the experimental data. Many body effects like local fields (LF), excitation induced dephasing (EID) and resonance energy shift (ES) are included phenomenologically bipul ; wegener ; wang  in this simplified approach. Although not done here, the coexcitation of hh and lh excitons can be easily included by modeling the excitonic assembly as a three-level system srink ; bipul . $`\delta `$-like pulses are assumed for computational ease although finite pulse widths can also be used. We calculate the time evolution of the PPDR signal \[$`\mathrm{\Delta }R(t,\tau )`$\] for both $``$ve and $`+`$ve $`\tau `$ (Fig. 5). Since the first term in Eq. (2) contributes only for $`+`$ve delay at $`t=0`$ in a time domain of the size of the pulse width, we display in Fig. 5 only the calculation of the second term in Eq. (2). We find that only EID contributes to $`\mathrm{\Delta }R(t,\tau )`$ apart from the phase space filling related nonlinearity. The peak behavior seen in Fig. 5 is in fact a consequence of EID in this model. The LF and ES related terms lead the EID term by a phase of $`\pi /2`$ and do not contribute to the TR-PPDR. (However, they do contribute to the SR-PPDR signal together with EID bipul .) In this model, the signal peak position in $`t`$ is essentially determined by an exciton density independent dephasing rate. The EID introduced in the Bloch equations does not lead to density dependent decay of $`𝒫^{(3)}(t,\tau )`$. This is consistent with the data of Fig. 3. Also, the model does lead to a delayed rise of the signal as observed in Fig. 2. Within this theory, the $`t`$ integration of the second term in Eq. (2) produces a nonzero signal for $`\tau <0`$, with the behavior $`\mathrm{exp}(2\tau /T_2)`$, $`T_2`$ being the exciton dephasing time. This agrees with the data plotted in Fig. 1(c), giving $`T_24`$ ps, which is close to the value obtained from FWM measurements on this sample.

However, some features in the data are not explained by the simple theory. For example, the origin of the different rising behavior in Figs. 2 and 3 for $``$ve and $`+`$ve $`\tau `$ is not clear. Furthermore, although the TR-PPDR signal decay appears to be not very sensitive to increase in the pump intensity (Fig. 3), the rise is. This is not accounted for in the above theory. Introduction of a phenomenological dependence of the dephasing time for the probe induced exciton polarization on pump generated excitation density may be considered, in analogy with the approach followed in the literature for FWM wang . This will indeed cause a faster rise of the TR-PPDR signal in $`t`$ with increased pump intensity ($`I_1`$). However, this will also lead to a faster decay in $`t`$ of the signal with increasing $`I_1`$, which is not evident in Fig. 3. A quantitative understanding of the $`t`$ and $`I_1`$ dependence exhibited by the data of Figs. 2 and 3 may be possible when an appropriate theory of the nonlinear optical interactions of the excitonic polarization, populations and propagating electric fields, that also takes into account the MQW nature of the sample and higher order Coulomb interactions, is applied. This is beyond the scope of this paper.

Among the important effects that may play some role in the TR-PPDR measurements is the following. It has been suggested stroucken ; baumberg ; prineas  that radiative coupling of excitations in individual QWs in a MQW structure may give rise to a few bright polaritonic modes which reradiate in the reflection direction more efficiently than other modes. One consequence of this theory is that the tail of the emission from the linear polarization may show a finite rise time and beats due to polaritonic interference. Indeed, it was deduced recently baumberg  from measurements of time integrated reflectivity of a pair of phase locked pulses using an interferometric technique that the excitonic reflection in GaAs MQWs shows a finite rise time at low excitation intensities.

To investigate this aspect in our case, the temporal profile of the reflected probe pulse is measured in absence of a pump pulse (Fig. 6). The strong peak at $`t=0`$ is essentially due to reflection of the probe pulse at the vacuum-sample interface. The oscillations seen are quantum beats due to coherent excitation of hh and lh excitons. The exciton polarization emits coherently for several ps. The emission is shown for two values of the average intensity $`I`$ of the incident pulse, namely $`I=200`$ and 5 $`\mu `$W. A simple estimate gives an exciton density of $`8\times 10^6`$ cm-2 per $`\mu `$W per QW. At low intensities, the emission shows a finite rise time, the signal peaking at about 3 ps. This behavior is consistent with the possibility that polaritonic interference has an important effect on the dynamics of excitonic polarization in MQWs and may influence TR-PPDR measurements. This feature is not seen for the high intensity case, possibly due to reduced photon induced coupling between the excitons in the individual QWs caused by exciton dephasing. However, even in this case, the decay of the signal shows a departure from a simple exponential behavior, as indicated in Fig. 6. We should mention that the $`t`$ dependence seen in Fig. 6 is different from the nonexponential decay that may be caused by exciton inhomogeneous broadening when it is comparable to or more than the homogeneous broadening.

To conclude, the TR-PPDR results presented here reveal new features related to the time dynamics of excitonic polarization in MQWs and the dependence of pump modulation effects on pump intensity, pump-probe relative polarization and delay. A nonlinear optical contribution arising from interference of excitonic linear and nonlinear polarization, which was not considered earlier in the literature, was shown to be essential for understanding the TR-PPDR measurements.

We are grateful to L. N. Pfeiffer and J. Shah for the MQW sample.


format: text

text:
# Dimensionality effects in Turing pattern formation

## I Introduction

Alan Turing is best known for his contribution to the foundations of computer science but he has also played a key role in the birth of nonlinear dynamical theory. As he pursued his dream of artificial brain he had to concern himself with the problem of biological growth. In 1952 he published a paper, where he proposed a mechanism by which genes may determine the structure of an organismturing . His theory did not make any new hypothesis, but it merely suggested that the fundamental physical laws can explain many of the emerging features. For example in case of a human body the problem is as follows: Can $`3\times 10^9`$ base-pairs of DNA code for approximately $`10^{11}`$ neurons, $`10^{15}`$ synaptic connections and estimated $`10^{13}`$ cells in total? In contrast to a layman’s belief that everything, i.e, structure, function and behaviour, is determined by genes, it is evident that there have to be some quite general physico-chemical mechanisms involved.

In order to develop a model for this kind of biological growth Turing was aware that he had to simplify things a lot, as becomes clear from the second sentence of his seminal paper: “This model will be a simplification and an idealization, and consequently a falsification.” In order to construct a manageable model he neglected the mechanical and electrical properties of tissue, and instead considered the chemical reactions and diffusion as the crucial factors in biological growth for formulating a descriptive theory as a system of coupled reaction-diffusion equations. Turing showed that this kind of system may have a homogeneous stationary state that is unstable against perturbations. In fact, any random deviation from this stationary state leads to a symmetry break and spatial concentration patterns due to a mechanism called diffusion-driven instability. In his paper Turing discusses the blastula stage of an embryo, which appears almost spherical but shows some deviations from the perfect symmetry to different (random) directions in each embryo of a certain species. Hence he stated that they can not be of great importance to morphogenesis: “From spherical initial state chemical reactions and diffusion certainly cannot result in an organism such as a horse, which is not spherically symmetrical.” Later he proved all this to be false because of the random deviations in the embryo, and stated: “It is important that there are some deviations, for the system may reach a state of instability in which these irregularities tend to grow.”

The fact that all the animals of a certain species do not have exactly the same coating patterns supports Turing’s idea that there is randomness involved in the morphogenesis. For example, all tigers have similar periodic patterns of stripes, but the stripes are not in the same exact positions in different tigers. Turing assumed the function of genes to be purely catalytic, and indeed the genome is known to be the blue-print for the biological structure that is to be formed. The genes of a tiger state that it should have stripes. As a result, the genes code proteins that set the reaction and diffusion rates of the morphogens in such a way that the homogeneous state of the chemical system becomes unstable due to local or temporal inhomogeneities, which eventually lead to complex pattern formation. Here a morphogen stands for a chemical related to biological pattern formation.

These reaction-diffusion, or Turing, systems have mostly been applied in mathematical biology for explaining pattern formation in biological systemsmurray . If one considers the Turing pattern of morphogens as a pre-pattern according to which the melanocytes (pigment producing cells) in the epidermis differentiate, one could explain how some animals get their nonuniform patterns. The problem is that the morphogens have not been identified in vivo and thus even their existence can be questioned. Yet, a lot of research has been done based on Murray’s hypothesis and patterns of butterfliessekimura , fishkondo ; barrio1 , and lady beetlesliaw  have been imitated by numerically solving Turing models. Apart from these numerical studies, also dynamic and stability aspects of the Turing models have been investigatedjudd ; lengyel ; callahan . In these, one applies the linear stability analysis or more complicated nonlinear bifurcation analysis to obtain insight of the system dynamics. Typically the parameters of a Turing model are adjusted on the basis of such an analysis. Despite the simple form of the equations theoretical analysis of the time-dependent evolution is very complicated and can only be done by analyzing the bifurcation from the stationary state. Therefore, the effect of introducing, e.g. inhomogeneous diffusion coefficientsbarrio2 , growing domainsvarea1  or curvature of the domainvarea2 , had to be studied numerically. As other extensions of the Turing systems we have studied the formation of Turing structures in three dimensionsteemu , the effect of noise to Turing structuresteemu2  and the dependence of the structural characteristics of the morphologies on the parameter selectionteemu3 . The first experimental observation of Turing patterns was presented in 1990 in a single-phase open gel reactor with chloride-iodide-malonic acid (CIMA) reactioncastets . For a review of pattern formation in biochemical systems, we refer the reader to the extensive article by Hesshess .

In the next section of this paper we explain the idea of Turing instability in a quite general way. After that we carry out linear analysis in the case of multiple stationary states and briefly discuss the possibility of bistability in the system. The parameter selection is followed by the presentation of numerical results of two- and three-dimensional Turing systems. The differences between 2D and 3D systems are considered with respect to the pattern selection and robustness against noise. In addition, preliminary results on the transition between two and three dimensions by increasing the thickness of the system are introduced and their connection to chemical experiments is discussed.

## II Diffusion-driven instability

Diffusion-driven instability or Turing instability is the mechanism by which the random motion of the molecules, i.e., diffusion may make a stable state of a chemical system unstable. Turing’s idea of diffusion-driven instability was ground-breaking for the following reasons: 1.) It is counter-intuitive: Typically diffusion stabilizes (e.g. a droplet of ink dispersing into water due to diffusion). 2.) From any random initial state diffusion-driven instability will result in the same morphology according to intrinsic characteristics of the system (such as the random deviations in the embryo). 3.) It set the basis for research in mathematical biology, nonlinear dynamical theory and chemical pattern formation; in 1950s non-equilibrium phenomena, symmetry-breaking or complex systems were not fashionable.

The reaction-diffusion mechanism resulting in an instability was illustrated by Turing with the problem of missionaries and cannibals in an island. Missionaries come to the island by boat and want to evangelize the cannibals. The rules are as follows: If two or more missionaries meet one cannibal they can convert him to a missionary. If the relative strength is the other way around, the missionaries get killed and eaten by the cannibals. As the missionaries die, more missionaries are brought to the island. In addition, cannibals reproduce cannibals. The mechanism for instability, i.e. diffusion, is in this example the movement of cannibals and missionaries. The missionaries are assumed to have bicycles and thus they move faster, i.e., they represent the inhibitor of a reaction by slowing down the reproduction of cannibals, which in turn are the activator. Should the missionaries not have bicycles, they would always get killed as they meet cannibals, but by having bicycles they have a chance to escape and return when there are more missionaries around. With these definitions the auto-catalytic nature of the Turing mechanism becomes evident: In areas with a lot of cannibals the number of cannibals will increase due to reproduction, then being more effective in killing missionaries. On the other hand, the predominance of the cannibals means that more missionaries will be brought to the island to convert them. In most cases the cannibals and missionaries will finally find a stationary pattern, which corresponds to a map of the island where the areas with cannibal dominance can be marked by one color and the areas with missionary dominance by another color.

## III Formal Turing models

A Turing model describes the time variation of the concentrations of two chemical substances or morphogens due to reaction-diffusion processes. Generally a Turing system can be written as

$`U_t`$ $`=`$ $`D_U^2U+f(U,V)`$
$`V_t`$ $`=`$ $`D_V^2V+g(U,V),`$ (1)
where $`U=U(\stackrel{}{r},t)`$ and $`V=V(\stackrel{}{r},t)`$ are the spatially varying concentrations, and $`D_U`$ and $`D_V`$ are the diffusion coefficients of the morphogens setting the time scales for the system. The reaction of the chemicals is described by the two nonlinear functions $`f`$ and $`g`$, which in general can be derived from chemical reaction formulas using the law of mass actionmurray .

In this paper we use the generic Turing model introduced by Barrio et al.barrio1  for which the reaction kinetics is obtained by Taylor expanding the nonlinear functions around a stationary solution $`(U_c,V_c)`$ defined by $`f(U_c,V_c)=0`$ and $`g(U_c,V_c)=0`$. From Eq. (1) one can see that this condition, indeed, defines the stationary state (time derivative equal to zero) in the absence of diffusion. If terms above the third order are neglected, the equations of motion read as follows

$`u_t`$ $`=`$ $`\delta D^2u+\alpha u(1r_1v^2)+v(1r_2u)`$
$`v_t`$ $`=`$ $`\delta ^2v+v(\beta +\alpha r_1uv)+u(\gamma +r_2v),`$ (2)
where $`u=UU_c`$ and $`v=VV_c`$. Terms $`r_1`$ and $`r_2`$ set the amplitudes of the nonlinearities, and they describe the reaction, such that quantity $`r_1`$ enhances stripe formation while $`r_2`$ enhances spots in two dimensionsbarrio1  and lamellae and spheres in three dimensionsteemu , respectively. $`D`$ is the ratio of the diffusion coefficients of the two chemicals, $`\delta `$ acts as a scaling factor, and $`\alpha `$, $`\beta `$ and $`\gamma `$ contribute to the mode selection. Here $`D1`$ is a necessary, but not sufficient condition for the diffusion-driven instability in two or more dimensionsvastano .

Previously, this model has been studied by setting $`\alpha =\gamma `$ to keep $`(0,0)`$ the only stationary solutionbarrio1 ; teemu . By relaxing this condition one can find two more stationary states defined by $`f(u_c,v_c)=g(u_c,v_c)=0`$, reading as follows

$$u_c^i=\frac{r_2+(1)^i\sqrt{r_2^2+4\alpha r_1(K\beta \gamma )}}{2K\alpha r_1}$$
(3)
and

$$v_c^i=Ku_c^i,$$
(4)
where $`K=\frac{\alpha +\gamma }{1+\beta }`$ and $`i=1,2`$.

Sufficient conditions for diffusion-driven instability to occur are widely knownmurray . In our case restricting the parameter selection such that $`\alpha (0,1)`$, $`\beta (1,0)`$ and $`\gamma (1,0)`$ we are left with only two conditions

$`1<\beta <\gamma ,`$
$`\alpha 2\sqrt{D|\gamma |}>\beta D.`$
Now the dispersion relation of the system (Eq. (2)) can be found by solving the eigenvalues of the linearized system of equations. This can be done by substituting a trial solution of the form $`w(\stackrel{}{r},t)=_kc_ke^{\lambda t}w_k(\stackrel{}{r},t)`$ into Eq. (2), which results in

$$\left(\begin{array}{c}\lambda u_k\\ \lambda v_k\end{array}\right)=\left(\begin{array}{cc}D_u& 0\\ 0& D_v\end{array}\right)\left(\begin{array}{c}k^2u_k\\ k^2v_k\end{array}\right)+\left(\begin{array}{cc}f_u& f_v\\ g_u& g_v\end{array}\right)_{u_c,v_c}\left(\begin{array}{c}u_k\\ v_k\end{array}\right),$$

where $`f_u`$, $`f_v`$, $`g_u`$ and $`g_v`$ denote the partial derivatives of the reaction kinetics evaluated at the stationary state $`(u_c,v_c)`$, which in the matrix form are as follows

$`\left(\begin{array}{cc}f_u& f_v\\ g_u& g_v\end{array}\right)_{u_c,v_c}=\left(\begin{array}{cc}\alpha (1r_1v_c^2)r_2v_c& 1r_2u_c2\alpha r_1u_cv_c\\ \alpha r_1v_c^2+\gamma +r_2v_c& \beta +2\alpha r_1u_cv_c+r_2u_c\end{array}\right),`$ (9)
where $`u_c`$ and $`v_c`$ are given by Eqs. (3) and (4). Now the dispersion relation $`\lambda (k)`$ can be solved from

$`\lambda ^2+\left[(D_u+D_v)k^2f_ug_v\right]\lambda +D_uD_vk^4k^2(D_uf_u+D_vg_v)f_vg_u=0,`$ (10)
where $`k^2=\stackrel{}{k}\stackrel{}{k}`$ and for the generic model $`D_u=\delta D`$, $`D_v=\delta `$.

If one chooses the parameters such that the two unstable states defined by Eqs. (3) and (4) are close to each other, preliminary numerical simulations show temporal changes in the concentration patterns as the system moves between the two stationary states. This bistability is most probably due to complex nonlinear coupling of the modes growing from the two unstable states. However, quantitative prediction of this interesting behavior requires nonlinear bifurcation analysis. The analysis and numerical results of the temporal solutions will be presented elsewhere. Following Ref. teemu , we fix $`\alpha =\gamma `$ from now on.

By observing the boundaries of the region with positive growth rate, i.e., $`\lambda (k)>0`$ one can analytically derive the modulus of the critical wave vector

$$k_c^2=\frac{1}{\delta }\sqrt{\frac{\alpha (\beta +1)}{D}}.$$

By adjusting the parameters and allowing only a few modes to be unstable one can end up with several different parameter sets. Here we use the parameters $`D=0.516`$, $`\alpha =0.899`$, $`\beta =0.91`$ and $`\delta =2`$, which correspond to critical wave vector $`k_c=0.45`$ and which we have used earlierteemu . These selections fix the characteristic length of the pattern to be $`2\pi /k_c`$. The corresponding dispersion relation can be calculated based on Eq. (10), which in the case of $`u_c=v_c=0`$ reduces to

$$\lambda ^2(\alpha +\beta \delta k^2(1+D))\lambda +(\alpha \delta Dk^2)(\beta \delta k^2)\gamma =0.$$
(11)
The real part of this dispersion relation is plotted in Fig. 1. For a more detailed discussion, see Ref. teemu .

## IV Numerical results

The numerical simulations were carried out by discretizing the system into a spatial mesh or lattice. The Laplacian was calculated by finite difference method, with $`dx=1.0`$, and the equations were iterated in time using Euler’s scheme with time step $`dt=0.05`$. In a discretized three-dimensional system, the wave number values are not continuous but of the form

$$|\stackrel{}{k}|=2\pi \sqrt{\left(\frac{n_x}{L_x}\right)^2+\left(\frac{n_y}{L_y}\right)^2+\left(\frac{n_z}{L_z}\right)^2},$$
(12)
where $`L_{x,y,z}`$ are the system dimensions to the corresponding directions and $`n_{x,y,z}`$ are the wave number indices. In this paper we have used periodic boundary conditions for the concentration fields, which were initially set as random perturbations around the stationary state, i.e., random numbers with zero mean and variance of $`0.05`$.

### 2D patterns and 3D structures

In the simulations the characteristic length of the pattern is fixed by the parameter selection and the morphology of the resulting pattern can be adjusted by the nonlinear parameters $`r_1`$ and $`r_2`$. The coefficient $`r_1`$ of the cubic term enhances stripe or lamellae formation, whereas $`r_2`$ enhances the formation of spherical shapes. Figure 2 shows the results of numerical simulations for two-dimensional $`128\times 128`$ systems, where the nonlinearities favour either spots or stripes.

Extending the pattern formation problem from 2D to 3D is by no means straightforward. Figure 3 shows 3D structures corresponding to the patterns of Fig. 2. The two-dimensional stripes become complex and aligned lamellae in the three-dimensional system instead of pure lamellar planes that one would expect (see Fig. 6). This is due to one more degree of freedom. Planes are formed, but the resulting structure is, as a matter of fact, a combination of aligned planes crossing each other. The system dynamics is unable to organize the three-dimensional structure into a more regular shape. The complex lamellar structure satisfies the symmetry requirements imposed by the nonlinearities though the resulting structure is not optimal. In the case of enhanced quadratic nonlinear interaction one obtains 3D spherical shapes as one would expect (Fig. 3), but the packing is not FCC.

### Dimensionality transition

The transition from 2D patterns to 3D structures is a more challenging problem than the direct comparison of the final results (Figs. 2 and 3). It has been shown experimentally that an open gel reactor with CIMA reaction may have a bistability, i.e., both spot and stripe patterns may appear into the gel on different heights. This has been explained by a concentration gradient that is imposed by the reactorquyang . Turing patterns have also been studied in ramped systems, where the thickness of the gel is increased gradually. In that case one observes qualitatively different patterns corresponding at different thicknessesdulos .

The problem has also been addressed more quantitatively by Dufiet and Boissonadedufiet . They modeled the pattern formation of a three-dimensional experimental reactor by imposing a permanent gradient on one of the bifurcation parameters of their model. This corresponds to the fact that in an experimental reactor the concentrations are kept constant only on the feed surfaces, whereas the concentration inside the gel is governed by reaction and diffusion. They claim that the patterns in quasi-2D reactors can, to a certain extent, be interpreted as two-dimensional patterns when the thickness of the gel $`L_z`$ is less than the characteristic wavelength $`\lambda _c=2\pi /k_c`$ of the pattern.

As for the effect of dimensionality in the generic Turing model of Eq. (2) we have observed that the transition from a 2D to a 3D system is not as simple as it has been thought. There seems to be some stochastic behaviour in the transition: As one increases the thickness of the system, it will lose the correlation between the bottom and top plates gradually, but the transition thickness is not well-defined, i.e., $`L_z>\lambda _c`$ does not guarantee that the structure becomes three-dimensional. Figure 4 shows the resulting structures for two systems of nearly the same thickness $`L_z>\lambda _c`$ as one starts from random initial conditions. One can easily observe that the leftmost structures is three-dimensional, whereas the rightmost is quasi-two-dimensional, although the latter structure grows in a thicker domain. The preliminary results presented in Fig. 4 were obtained by straightforward simulations of the Lengyel-Epstein modellengyel , as follows

$`u_t`$ $`=`$ $`{\displaystyle \frac{1}{\sigma }}(au4{\displaystyle \frac{uv}{1+u^2}}+^2u)`$
$`v_t`$ $`=`$ $`b(u{\displaystyle \frac{uv}{1+u^2}})+d^2v,`$ (13)
where $`a`$, $`b`$, $`d`$ and $`\sigma `$ are adjustable parameters. We used the parameter set $`\sigma =50`$, $`d=1.07`$, $`a=8.8`$ and $`b=0.09`$, which is known to correspond to hexagonal patternsrudovics .

Based on our numerical studies using both the generic Turing model (Eq. (2)) and the Lengyel-Epstein model (Eq. (13)) it seems that the probability that the pattern is three-dimensional is proportional to the number of wave vectors $`\stackrel{}{k}`$ (Eq. (12)) with $`n_z0`$ such that the mode $`k`$ is unstable according to the dispersion relation (Eq. (11)). Our results show that the number of three-dimensional systems does not change continuously as a function of the system thickness, but the dependence is very complicated. It should also be noted that linear analysis is not sufficient to explain the complex dynamics of the dimensionality transition. Thus some sort of non-linear analysis seems necessary. However, from the numerical simulations one can calculate interesting statistics of the morphologies (e.g. correlations, structure factors), although this kind of structural analysis is difficult to carry out for chemical patterns in a gel. In addition visualization of numerical results makes it possible to see into the structure. We have used the Lengyel-Epstein model reactionlengyel ; rudovics  to analyze further the structures that have been reported from experiments in ramped systemsdulos . It seems that the experimental patterns can be easily misinterpreted, because the depth information is lost in the 2D projection made for visualization. What seems to be a non-harmonic modulation may be just an aligned lamellae, and what seems to be a combination of stripes and spots may as well be tubes, which appear to the observer as stripes if seen from the side and as spots if seen from the end.

As discussed above an important feature of experiments, which is absent in the numerical simulations is the gradients in the reactor. In a computer a given model can be solved to the computational precision, from which it follows that if a model captures the phenomenon accurately, the results are accurate and there are no such artifacts as concentration gradients. If one thinks of applying Turing system in biological modeling, it is of great importance that the models imitate the real processes and there is nothing that skews the results. In biological tissue the sources controlling the parameters of the reaction would in part be in the cells and not always in the boundaries. Thus the gradients are not always present in biological morphogenesis, from which it follows that the patterns obtained from numerical simulations may in fact be more realistic than the patterns arising in experimental reactors.

### The effect of noise on lamellar structures

We have recently reported results of the general Turing system concerning the effect of noise on spherical structuresteemu2 , which are more robust against noise than stripes or lamellar structures. Here we present results on the effect of noise specifically in the presence of cubic nonlinear interaction, which favours stripes or lamellae. In order to study this effect we introduce uncorrelated Gaussian noise sources $`\eta (\stackrel{}{x},t)`$ such that the equations of motion read as follows

$`u_t`$ $`=`$ $`D\delta ^2u+f(u,v)+\eta _u`$
$`v_t`$ $`=`$ $`\delta ^2v+g(u,v)+\eta _v,`$ (14)
where the first and the second moment of the noise is defined as $`\eta (\stackrel{}{x},t)=0`$ and $`\eta (\stackrel{}{x},t)\eta (\stackrel{}{x}^{},t^{})=A^2\delta (\stackrel{}{x}\stackrel{}{x}^{})\delta (tt^{})`$, with the angular bracket denoting the average and $`A`$ the noise intensity. It is noted that the noise was added to each lattice site at every time step of the simulation. Due to discretization the noise has to be normalized such that

$$\eta =\frac{A}{(dx)^{d/2}\sqrt{dt}},$$
(15)
where $`d`$ is the dimension of the system, $`dx`$ the lattice constant and $`dt`$ the time step.

In Figure 5 we show the two-dimensional patterns corresponding to different noise intensities for nonlinear parameters $`r_1=3.5`$ and $`r_2=0`$. The noise was applied all the way throughout the evolution, and one can observe that the stripes fade into the noise gradually. It should be noted that there are complex dynamics involved in the evolution, but the patterns are still formed according to the system parameters, and even a substantial noise does not distort the evolution completely, but only makes the resulting patterns look noisy.

We have shown earlier that by using a special initial condition in the mid-plane of the system, one can obtain a three-dimensional structure of pure parallel planesteemu . Figure 6 shows such a planar structure and the same structure in the presence of noise. One can clearly see that the planes get holes on them and align even in the presence of a small amount of noise. The noise intensity in the rightmost structure of Fig. 6 corresponds to that of Fig. 5B, where the stripes are intact. This explains the fact that the three-dimensional pattern selection is very sensitive and the pure planes are not obtained from random initial conditions.

In general, we find that the three-dimensional case is much more robust against noise than the two-dimensional one, as can be observed with the help of a structure factor analysisteemu2 . The interfaces of the chemical domains become distorted, but one can still localize the structure even for noise intensities up to 50% of the amplitude of the chemical concentration wave. Such high level robustness of Turing patterns is very interesting as such but also from the point of view of applying Turing systems to biological growth, in which morphogenesis seem to sustain well against external random distortions.

## V Discussion

In this paper, we have studied the characteristics of two-dimensional, quasi-two-dimensional and three-dimensional pattern formation in Turing systems. Our preliminary results indicate that the system may behave as a three-dimensional system for certain thickness and as one increases the thickness, the system may behave quasi-two-dimensionally with the same set of parameters. Also we find that the pattern selection in a 3D system is much more complex than in 2D. Although spherical structures are analogous to spotty patterns, similar generalization from stripes to lamellae is not as easy. A simulation of a three-dimensional system in the presence of a cubic nonlinear interaction results in aligned lamellae with holes. This may be due to the fact that it is more difficult to maintain long-range order in three dimensions, which is required for lamellae without holes. One can observe the breaking of long-range order also in the case of spherical structures: The two-dimensional spots are organized in an almost perfect triangular lattice, whereas in 3D spherical structures seem to be unable to establish regular ordering.

In the study of the effect of noise on striped and planar structures we found that these structures and the dynamics of the instability are proportionally speaking quite robust against random perturbations. However, spherical structures can sustain more noise in absolute values since the amplitude of the concentration wave of a spherical structure is larger. We have also shown that the three-dimensional planar structure did not sustain noise as well as the corresponding two-dimensional stripe pattern although three-dimensional structures should be more robust than two-dimensional patterns as indicated by our preliminary studies.

In this paper we have also found that the dimensionality can affect the pattern formation process in a very profound way and thus understanding the differences between two-dimensional pattern formation and three-dimensional structure formation is of great importance. The real biological processes are always taking place in a three-dimensional domain, although the symmetry restrictions, instability and dynamics governing the process may actually be quasi-two-dimensional. Therefore, choosing whether a system should be treated as three-dimensional or quasi-two-dimensional is by no means trivial and thus needing further theoretical and numerical studies for deeper understanding of morphogenesis modeling.

## Acknowledgements

One of us (R. A. B.) wishes to thank the Laboratory of Computational Engineering at Helsinki University of Technology for their hospitality. This work has been supported by the Finnish Academy of Science and Letters (T. L.). and the Academy of Finland through its Centre of Excellence Program (T. L. and K. K.).

Animations of Turing systems are available at
http://www.lce.hut.fi/research/polymer/turing.shtml


text:
# Comparing and combining Wilkinson Microwave Anisotropy (WMAP) probe results and Large Scale Structure

## I Introduction

The first year results from the Wilkinson Microwave Anisotropy probe (WMAP) were announced on February 11 2003. On the same day the telescope was renamed in honor of Prof. David Wilkinson, member of the science team and pioneer in the study of cosmic microwave background (CMB) radiation. The primary goal of the WMAP mission is to produce a high-fidelity all-sky polarization-sensitive map of the cosmic microwave background (CMB) radiation to determine the cosmology of our Universe. After a year of observation it has produced a full sky map of the microwave sky in 5 frequencies, with a resolution a factor 30 higher than the previous full sky map as produced by the COBE satellite Bennettetal96 . This is the cleanest picture of the early Universe; the structures on the CMB –the pattern of hot and cold spots– carry information about the composition, geometry, age etc. of our Universe (e.g., Jungmanetal96a ; Jungmanetal96b ). Here I will illustrate how it is possible to constrain inflationary models by combining WMAP’s CMB observations with external data sets: in particular, the smaller scales CMB experiments CBI and ACBAR, and observations of large-scale structure probes such as the power spectrum of Two-degree-field galaxy redshift survey (2dFGRS; Collessetal01 ) and the linear matter power spectrum as recovered by Croftetal02  and Gnedinhamilton02  from Ly$`\alpha `$ forest observations. We also consider the Type 1A supernoave results Riessetal01  and HST Key project Hubble constant determinations Freedmanetal01 . The results of this analysis were reported in Verdeetal03 ; Spergeletal03 ; Peirisetal03 .

## II WMAP

If the CMB is Gaussian (see Komatsuetal03 ) then the information in the mega-pixel maps can be compressed in the power spectrum Hinshawetal03 . We show the WMAP temperature power spectrum in Fig. 1 (left), where the points with error-bars are the (band power) data with the error given by the instrumental noise. The solid line is our best fit model and the gray area shows the cosmic variance uncertainty. At the moment we can reliably measure only the temperature-polarization (TE) cross-correlation power spectrum which is shown in Fig. 1 (right) where the points with error-bars are the band-power data. The temperature power spectrum makes precise predictions for the TE power spectrum on large scales. The solid line in Fig. 1 (right) is this prediction from the best fit to the TT data only. For $`\mathrm{}>10`$ the data are in remarkable agreement with the TT predictions, this is a triumph for the standard cosmological model (see Peirisetal03 ); the WMAP detection of TE anti-correlation at $`\mathrm{}50150`$ rules out a broad class of active models and implies the existence of super-horizon adiabatic fluctuations at decoupling: both inflation and the Ekpyrotic scenario predict the existence of super horizon fluctuations.

For $`\mathrm{}<10`$ we detect a signal that cannot be seen in the TT data alone: this is the signal of early re-ionization Kogutetal03 . Such a strong signal at low multipoles gives a constraint on the integrated optical depth to the last scattering surface $`\tau =0.17\pm 0.04`$, for instantaneous reionization this corresponds to a reionization redshift of $`z_{re}=17\pm 3`$.

Why WMAP has such a good signal to noise? A part from the fact that it is an all sky experiment which observed the sky for a full year, this is the result of an obsession of all those involved in designing and building the satellite for minimizing systematics and maximizing the signal to noise. During the designing and building of the satellite the emphasis was on minimizing systematic errors so that the data could be straightforwardly interpreted. For example there are 10 different channels from 20 to 95 GHz, so that foreground subtraction is not a problem Bennettetal03 ; Hinshawetal03 . Undesirable $`1/f`$ noise is minimized by the design of WMAP radiometers Jarosiketal03 , uncorrelated noise between pixels, accurate determination of the beam patterns Pageetal03 ; Barnesetal03  and well understood properties of the radiometers are invaluable for this analysis.

The rest of the analysis must therefore have the same level of “obsession”. For example, to extract as much information as possible from the maps without introducing systematics, the TT power spectrum is obtained by optimally combining “only” 28 cross-channels correlations. This method has then be checked against another 2 independent methods Hinshawetal03 . Also, the point source subtraction has been done in 3 different ways Hinshawetal03 ; Bennettetal03 ; Komatsuetal03 (they all agree!). Finally in the covariance matrix we propagate all the effects we know about (beam errors, noise, sky cut etc.) To cut a long story short, we spent a long time worrying about 2% error on the error! There are effects that we are aware of, which we neglect; these effects give about 0.5 to 1% error on the error.

## III The Method

The analysis method must match the level of accuracy and precision of the data. If the CMB sky is Gaussian, then the spherical harmonics coefficients $`a_\mathrm{}m`$ are normally distributed, but the angular power spectrum, $`C_{\mathrm{}}`$, is not. For large $`\mathrm{}`$ the Central Limit Theorem will make the Gaussian likelihood an increasingly good approximation, but the error-bars in the WMAP data are so small, that previous likelihood approximations extensively used in the literature are not suitable. We introduce a different form of the likelihood function that is accurate to better than 0.1%. The likelihood function and its covariance matrix has been calibrated on 100,000 Monte Carlo realization of the sky as seen by WMAP Verdeetal03 . This calibration is also important because we can now compute the reduced $`\chi ^2`$ for the “best fit model” and use it to compute the probability that the temperature pattern observed in the sky is a realization of the best fit model.

We perform a likelihood analysis for the cosmological parameters: at each point in parameter space (i.e., for a given set of cosmological parameters) a new model $`C_{\mathrm{}}`$ is computed using the publicly available code CMBfast 4.1 cmbfast . Even with a supercomputer, at 2 seconds per evaluation, with a naive grid-based approach it would take a prohibitive amount of CPU time to perform the analysis for 6 or more parameters. We resort to Markov Chain Monte Carlo (MCMC): this approach has become the standard tool for CMB analyses CM ; Cetal ; LB ; Kosowskyetal02 . In fact, for a flat reionized universe we can evaluate the likelihood $`120,000`$ times in $`<2`$ days: this is adequate for finding the best fit model and reconstruction of the 1- and 2- $`\sigma `$ confidence levels for the cosmological parameters. While more details are in Verdeetal03  and references therein, here we just remind the reader that a MCMC is similar to a random walk in parameter space: at each step a new point is sampled in parameter space but the density of steps (points in parameter space) is proportional to the a posteriori likelihood. After an initial “burn-in” period, the chain has no dependence on the starting location. The important point to keep in mind about MCMC is that a MCMC cannot tell us about regions of parameter space where it has not been: it is very important that the chain achieves good mixing (i.e., moves rapidly in parameter space) and fully explore the likelihood surface. Using an MCMC that has not fully explored the likelihood surface for determining cosmological parameters will yield wrong results, thus it is important to have a convergence criterion and mixing diagnostic. For any analysis of the WMAP data we strongly encourage the use of a convergence criterion, for example we use the method proposed by GR92 .

We speed up convergence and improve mixing by using the physical parameters of BE  and Kosowskyetal02  (e.g., physical cold dark matter density $`\mathrm{\Omega }_ch^2`$ instead of $`\mathrm{\Omega }_c`$ and the characteristic angular scale of the acoustic peaks $`\mathrm{\Theta }_A`$ instead of the Hubble constant).

### III.1 WMAPonly results

We first present the cosmological parameter constraints from WMAP data only. Then we consider external data sets (smaller scale CMB experiments and large scale structure power spectrum measurements) and compare these data to the predictions of the best fit model for WMAP data. We find that this extrapolation of WMAP predictions to low redshift and smaller scales works remarkably well. We thus proceed to perform a joint likelihood analysis of WMAP data with these external data sets. Apart from the fact that the CMB seems to be Gaussian Komatsuetal03 , consistent with the generic prediction of inflationary models, and that the Universe seems to have been reionized much earlier that previously though Kogutetal03 , the most striking result, probably, is that the “vanilla”, flat LCDM model still fits the data amazingly well: 6 parameters fit 1348 data-points (these are the power spectrum multipoles, one could be bolder and say that 6 parameters fit $`10^6`$ data, these being the pixels!). The parameters that best fit WMAP data are summarized in table 1.

Remarkably, this same model fits not only WMAP data but also other small scale CMB experiments and a host of other cosmological observations: the Hubble constant constraint ($`H_0=72\pm 5`$) is in agreement with the Hubble key project Freedmanetal01  determination of $`H_0=72\pm 3\pm 7`$ at the 68% confidence level(C. L.); the age determination of $`13.4\pm 0.3Gyr`$ is in good agreement with other age-determination based on stellar ages; the baryon abundance $`\mathrm{\Omega }_bh^2=0.024\pm 0.001`$ implies a primordial deuterium abundance in good agreement with D/H determination from quasar absorption lines; the determination of amplitude of fluctuations parameterized by $`\sigma _8`$ (the r.m.s. fluctuations smoothed on a scale of $`8Mpc/h`$) is in good agreement with clusters and weak lensing estimates.

## IV External data sets

We also consider smaller scales CMB experiments ACBAR Kuoetal02  and CBI Masonetal02 ; Pearsonetal02  on scales where they do not overlap with WMAP measurements and they are not affected by possible SZ signal. We combine CMB data with measurement of the low redshift universe. Galaxy redshift surveys allow us to measure the galaxy power spectrum at $`z0`$ and observations of the Lyman $`\alpha `$ (Ly$`\alpha `$) forest allows us to probe the dark matter power spectrum at $`z3`$. We use the Anglo-Australian Telescope Two Degree Field Galaxy redshift Survey (2dFGRS;Collessetal01 ), and the linear matter power spectrum as recovered by Croftetal02  and Gnedinhamilton02  from Ly$`\alpha `$ forest observations. When analyzing the 2dFGRS survey we find that we need to model the non-linear evolution of the power spectrum and non-linear redshift space distortions, even when limiting the analysis to a regime traditionally considered to be still well described by linear theory. We tested our modeling extensively on mock catalogs of the 2dFGRS. The goal of our modeling is to relate the shape and the amplitude of the observed power spectrum to that of the underling linear dark matter power spectrum as constrained by CMB data. Thus we use the information about the redshift space-distortion parameters and galaxy bias as in Peacoketal01  and Verdeetal02  (for more details see Verdeetal03 ).

For the matter power spectrum derived from Ly$`\alpha `$ forest observations we follow the approach on Gnedinhamilton02  and increase the error-bars where their analysis is not in agreement with that of Croftetal02 ). We stress here that the addition of Ly$`\alpha `$ forest data appear to confirm trends seen in the other data sets and tighten cosmological constraints however our results do not rely on Ly$`\alpha `$ data. More observational and theoretical work is still needed to confirm the emerging consensus that the Ly$`\alpha `$ forest data traces the large scale structure (see L. Hui and U. Seljack talks)

The key cosmological parameters that are relevant for single field inflation are: the constraint on the total matter-energy density of the universe relative to the critical density $`\mathrm{\Omega }_{tot}`$, the spectral slope of the primordial power spectrum $`n_s`$ the running of the spectral index $`dn/d\mathrm{ln}k`$ and the tensor to scalar ratio $`r`$. These are reported in table 1 and in figure .

## V Results and conclusions

I will present here 2 sets of conclusions:

* 1) for the pessimists I have shown how to pollute the clean CMB data with the messy gastrophysics of the low redshift probes of large-scale structure.
* 2) for the optimists, I have shown that all these different and independent data sets are in remarkable agreement. Within the context of the inflation + LCDM paradigm cosmological parameters constraints obtained from these different datasets are consistent. The results obtained at $`z1088`$ can be extrapolated forward to $`z0`$ and describe a host of astronomical observations.

This gives us the confidence to extrapolate our results backwards, to well before $`z1088`$ and place constraints on Inflation, as presented in Peirisetal03  (see also D. Spergel talk and fig. 3).

For both optimists and pessimists I hope I have illustrated how it is possible to enhance the scientific value of CMB data from $`z1088`$ by combining it with measurements of the low redshift universe. When including these external data sets one should keep in mind that the underlying physics for these data sets is much more complicated and less well understood than for CMB data and especially WMAP data. For both pessimists and optimists I hope I have highlighted the importance of improving the rigor of large-scale structure data analysis and of the interpretation of these observations, and how crucial it is to quantify an d account for the level of systematic contamination.

The data set, the WMAP temperature and temperature-polarization angular power spectra are available at http://lambda.gsfc.nasa.gov. If using these data set please refer to Hinshaw et al. (2003) for the TT power spectrum and Kogut et al. (2003) for the TE power spectrum. We also make available a subroutine that computes the likelihood for WMAP data given a set of $`C_{\mathrm{}}`$, please refer to Verde et al. (2003) if using the routine.

## VI Acknowledgments

I thank my collaborators in this work: C. Barnes, C. Bennett (PI), M. Halpern, R. Hill, G. Hinshaw, N. Jarosik, A. Kogut, E. Komatsu, M. Limon, S. Meyer, N. Odegard, L. Page, H. Peiris, D. Spergel, G. Tucker, L. Verde, J. Weiland, E. Wollack, E. Wright. (the WMAP science team). The WMAP mission is made possible by the support of the Office of Space Science at NASA Headquarters. LV is supported by NASA through Chandra Fellowship PF2-30022 issued by the Chandra X-ray Observatory center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073. The WMAP mission is made possible by the support of the Office of Space Science at NASA Headquarters. I thank the 2dFGRS consortium for help and discussion.


text:
# 1 Introduction

## 1 Introduction

The explanation of the anomalous perihelion shift of Mecury’s orbit represents one of the most famous classical tests of General Relativity (GR)- the relativistic gravitation theory in curved space-time. This is the only one involving relativistic effects on massive bodies - all other classical tests of GR being based on light propagation effects only. The perihelion advance of Mercury had been an unsolved problem in celestial mechanics for over half a century, since the announcement by Leverrier in $`1859`$  that, after the perturbing effects of the precession of the equinoxes on the astronomical co-ordinate system had been subtracted, there remained in the data an unexplained advance in the perihelion of Mercury. The modern value of this discrepancy is $`43`$ arc-seconds per century. A number of ad hoc proposals were made in an attempt to account for this excess, including , among others, to postulate  a magnetic-type component in gravity ( the gravitomagnetic field) for the gravitational influence of the Sun on the motion of planets ( the magnitude of this ad hoc component could be adjusted so as to account for the excess perihelion motion of Mercury), the existence of a new planet Vulcan near the Sun, a ring of planetoids, a solar quadrupole moment and a deviation from the inverse-square law of gravitation, but none was successful. The great interest in this paradoxical result vanished in 1915 when Einstein in his GR showed that the excess advance of Mercury’s perihelion could be explained as a relativistic effect in curved space-time. In the GR the non-Newtonian “excess” advance of perihelion of Mercury’s orbit is explained by using space curvature and Schwartzschild metric. Since the special relativistic approach to the problem of perihelion advance in the Kepler motion could not yield the observed effect (the existing approaches yield only one sixth of Einstein’s predicted value in GR, see for example  ), the observation of this effect played a role of one of the seemingly successful tests of the GR. In this work we, by the way reporting the establishment of a compatibility of Newtonian gravity with special relativity in a non-general-relativistic way, offer a flat space-time relativistic explanation for the observed anomalous perihelion advance of Mercury from a new angle. In this new approach to the old and general relativistically solved problem we make use of the notions of gravitational Lienard-Wiechert potentials (GLWP) and the gravitational Thomas precession(GTP) which naturally follow from some extended special relativistic considerations. To this end we devote Sec.$`2`$ to the notion of gravitational Lienard-Wiechert potential (GLWP). In Sec.$`3`$ we show that the GLWP together with the GTP invoked in  is sufficient to explain the observed perihelion advance of Mercury within the framework of a flat space-time relativistic toy model vector theory of gravity named as “Maxwellian Gravity” by the authors of .In Sec.$`4`$ we have some concluding remarks of interest for relativistic gravity.

## 2 The notion of GLWP

The notion of gravitational Lienard-Wiechert potential (GLWP) is the gravitational analogue of Lienard-Wiechert potential in classical electromagnetism. So to invoke the GLWP one has to establish a correspondence between the gravitation theory and electromagnetic theory. The close formal analogy between Newton’s law of gravitation and Coulomb’s law of electricity led many authors, in the past and also more recently , to investigate further similarities , such as the possibility that the motion of gravitational mass/charge could generate the analogs of a magnetic field.The magnetic field is produced by the motion of electric-charge,i.e. the electric current: the motion of gravitational mass/charge would produce what is called “gravitomagnetic” field. It is to be noted that such a correspondence had been proposed by J.C. Maxwell in 1865 and later further studied by Oliver Heaviside and O. Jefimenko.In the linearized versions of Einstein’s gravitational field equations, Maxwell-Lorentz-like equations for gravity have been obtained by several authors, see for examples. However from a closer look at these equations the reader can find that the linearized theory of GR is not perfectly isomorphic with electromagnetism which is commonly understood as a limitation of Linearized GR. On the other hand there is no such limitation in the linear relativistic gravity developed in  under the name of “Maxwellian Gravity”(MG).In its formulation the relativistic nature of gravity and the source of gravity have been re-investigated by re-designing an often cited  thought experiment  involving the motion of two point-like charged particles in Minkowski space-time.In the thought experiment we have considered a system of two point-like charged particles having such amount of rest masses and charges that the force of electrostatic repulsion balances that of gravitational attraction between the particles in an inertial frame of reference $`K^{}`$ in which both the particles are at rest under equilibrium condition. Then we investigated in the condition of equilibrium of the said particle system in $`K^{}`$-frame as well as in another inertial frame $`K`$ moving with uniform velocity with respect to the $`K^{}`$-frame of reference. From the requirement of the frame-independence of the equilibrium conditions, we not only obtained a Lorentz-force law for gravitational interaction between the moving masses as expected  but also unexpectedly found the Lorentz-invariant rest masses of the interacting particles as representing their gravitational masses (or gravitational charges) in complete analogy with the Lorentz force law and the invariant electric charges of the classical electromagnetic theory. These findings are in conformity with Poincaré’s remark that if equilibrium is to be a frame-independent condition, it is necessary for all forces of non-electromagnetic origin to have precisely the same transformation law as that of the Lorentz-force. Having recognized these findings we have obtained four Faraday-Maxwell-type toy-model linear equations of gravity describing what we call “Maxwellian Gravity” following the known procedures of the electromagnetic theory. The equations have a surprisingly rich and detailed correspondence with Faraday-Maxwell’s field equations of the electro-magnetic theory. The field equations can be written in the following Faraday-Maxwellian-form:

$$\stackrel{}{}\stackrel{}{E}_g=4\pi G\rho _0=\rho _0/ϵ_{0g}whereϵ_{0g}=1/4\pi G$$
(1)
$$\stackrel{}{}\times \stackrel{}{B}_g=\mu _{0g}\stackrel{}{j}_0+(1/c^2)(\stackrel{}{E}_g/t),where\mu _{0g}=4\pi G/c^2$$
(2)
$$\stackrel{}{}\stackrel{}{B}_g=0$$
(3)
$$\stackrel{}{}\times \stackrel{}{E}_g=\stackrel{}{B}_g/t$$
(4)
Where $`\rho _0`$ = rest mass (or proper mass) density; $`\stackrel{}{j_0}`$ = rest mass current density; $`G`$ is Newton’s universal gravitational constant; $`c`$ is the speed of light in empty space; the gravito-electric and gravito-magnetic fields $`\stackrel{}{E_g}`$ and $`\stackrel{}{B_g}`$ respectively are defined by the gravitational Lorentz force on a test particle of rest mass $`m_0`$ moving with uniform velocity $`\stackrel{}{u}`$ as

$$\frac{d}{dt}[m_0\stackrel{}{u}/(1u^2/c^2)^{1/2}]=m_0[\stackrel{}{E_g}+\stackrel{}{u}\times \stackrel{}{B_g}]$$
(5)
where the symbols have their respective meanings in correspondence with the Lorentz force law in its relativistic form.It is to be noted that these set of gravitational Maxwell-Lorentz equations coincide with those speculated by Maxwell, Heaviside and discussed also by Peng in the weak field and slow motion limit of Einstein’s field equations. However the treatment made in  suggests the validity of these equations for relativistic speeds as well. In covariant formulation, introducing the space-time four vector $`x_\mu =(x,y,z,ict)`$, proper mass current density four vector $`j_\mu =(j_{0x},j_{0y},j_{0z},ic\rho _0)`$ and the second-rank antisymmetric gravitational field strength tensor

$`F_{\mu \nu }=\left(\begin{array}{cccc}0& B_{gz}& B_{gy}& iE_{gx}/c\\ B_{gz}& 0& B_{gx}& iE_{gy}/c\\ B_{gy}& B_{gx}& 0& E_{gz}/c\\ iE_{gx}/c& iE_{gy}/c& iE_{gz}/c& 0\end{array}\right)`$ (6)
The field equations (1-4) can now be represented by the following two equations:

$$\underset{\nu }{}F_{\mu \nu }/x_\nu =\mu _{0g}\stackrel{}{j}_\mu ,\mathrm{where}\mu _{0g}=4\pi G/c^2$$
(7)
$$F_{\mu \nu }/x_\lambda +F_{\nu \lambda }/x_\mu +F_{\lambda \mu }/x_\nu =0$$
(8)
while the gravitational Lorentz force law (5) assumes the form :

$`c^2(d^2x_\mu /ds^2)=F_{\mu \nu }(dx^\nu /ds)`$ (9)
The absence of the rest mass of the test particle in its co-variant equation of motion (9) in the external gravitational field $`F_{\mu \nu }`$ describes clearly the universal nature of gravitational interaction in conformity with Galileo’s empirical law of universality of free fall (UFF) in a uniform gravitational field in the relativistic case as well. The law of UFF is one of the most fascinating character of gravity. For the Newtonian theory, the universality of free fall is of course a consequence of the equality of inertial mass and gravitational mass. Being fascinated by the mystery of the UFF observed in Galileo-Newtonian physics, Einstein in his development of the GR postulated that it holds generally, in particular also for large velocities (the relativistic case) and strong fields and for any form of mass-energy (i.e. he postulated the equality of inertial mass and gravitational mass for the UFF to hold generally). It is to be carefully noted that the general validity of the UFF is now a natural consequence of Maxwellian Gravity without the requirement of an Einsteinian postulate on the equality of inertial and gravitational masses. This is one of the interesting theoretical revelations of MG exacting the real UFF. The origins of the UFF may be traced to the Lorentz-invariant nature of the gravitational mass-charge (i.e. the rest mass or any form of Lorentz-invariant mass-energy) and the laws of physics as dictated by special relativity.

The fields $`\stackrel{}{E}_g`$ and $`\stackrel{}{B}_g`$ of MG are derivable from potential functions

$$\stackrel{}{B}_g=\times \stackrel{}{A}_g,\stackrel{}{E}_g=\mathrm{\Phi }_g\stackrel{}{A}_g/t$$
(10)
where $`\mathrm{\Phi }_g`$ and $`\stackrel{}{A}_g`$ represents respectively the gravitational scalar and vector potential of MG.These potentials satisfy the inhomogeneous wave equations :

$$^2\mathrm{\Phi }_g\frac{1}{c^2}\frac{^2\mathrm{\Phi }_g}{t^2}=\mathrm{\hspace{0.17em}4}\pi G\rho _0=\rho _0/ϵ_{0g}$$
(11)
$$^2\stackrel{}{A}_g\frac{1}{c^2}\frac{^2\stackrel{}{A}_g}{t^2}=\frac{4\pi G}{c^2}\stackrel{}{j}_0=\mu _{0g}\stackrel{}{j}_0$$
(12)
if the gravitational Lorenz gauge condition

$$\stackrel{}{}\stackrel{}{A}_g+\frac{1}{c^2}\frac{\mathrm{\Phi }_g}{t}=\mathrm{\hspace{0.17em}0}$$
(13)
is imposed. These will determine the generation of gravitational waves by prescribed gravitational charge and current distributions. Particular solutions (in vacuum) are

$$\mathrm{\Phi }_g(\stackrel{}{r},t)=G\frac{\rho _0(\stackrel{}{r}^{},t^{})}{|\stackrel{}{r}\stackrel{}{r}^{}|}𝑑v^{}$$
(14)
$$\stackrel{}{A}_g(\stackrel{}{r},t)=\frac{G}{c^2}\frac{\stackrel{}{j}_0(\stackrel{}{r}^{},t^{})}{|\stackrel{}{r}\stackrel{}{r}^{}|}𝑑v^{}$$
(15)
where $`t^{}=t|\stackrel{}{r}\stackrel{}{r}^{}|/c`$ is the retarded time.These are called the retarded potentials. Thus we saw that retardation in gravity is possible in Minkowski space-time in the same procedure as we adopt in electrodynamics.This result seems to conflict with the view  that Newtonian gravity is entirely static, retardation is not possible until the correction due to deviations from Minkowski space is considered. Now in analogy with the electromagnetic case we have here the Gravitational Lienard-Wiechert potentials for a point particle of rest mass $`m_0`$ moving with velocity $`\stackrel{}{v}`$ as

$$\mathrm{\Phi }_g=\frac{Gm_0}{r(\mathrm{\hspace{0.17em}1}\frac{v^2\mathrm{sin}^2\theta }{c^2})^{1/2}}$$
(16)
$$\stackrel{}{A}_g=\frac{Gm_0\stackrel{}{v}}{c^2r(\mathrm{\hspace{0.17em}1}\frac{v^2\mathrm{sin}^2\theta }{c^2})^{1/2}}$$
(17)
where $`r`$ is the magnitude of the instantaneous position vector $`\stackrel{}{r}`$ of the field point from the position of the particle and $`\theta `$ is the angle between $`\stackrel{}{r}`$ and $`\stackrel{}{v}`$ at the instant of time. Thus for a planet with rest mass $`m_0`$ having relative velocity $`\stackrel{}{v}`$ with respect to the the Sun (with rest mass $`M_{}`$) the instantaneous gravitational Lienard-Wiechert (LW) potential energy is given by

$$U_{gLW}=\frac{GM_{}m_0}{r(\mathrm{\hspace{0.17em}1}\frac{v^2\mathrm{sin}^2\theta }{c^2})^{1/2}}$$
(18)
Considering the angular momentum of the planet as $`L=m_0rv\mathrm{sin}\theta `$ Eq.(18) can be re-written as

$$U_{gLW}=\frac{GM_{}m_0}{r(\mathrm{\hspace{0.17em}1}\frac{L^2}{m_0^2c^2r^2})^{1/2}}$$
(19)
which in the first approximation can be reduced to

$$U_{gLW}=\frac{GM_{}m_0}{r}\frac{GM_{}L^2}{2m_0c^2r^3}=\frac{k}{r}\frac{h_0}{r^3}$$
(20)
where

$$k=GM_{}m_0\text{ and}h_0=\frac{GM_{}L^2}{2m_0c^2}$$
(21)
Thus we saw that the consideration of the GLWP introduced a $`\mathrm{\hspace{0.17em}1}/r^3`$ potential into the Kepler problem. Its effect on the planetary motion will be considered in the following section together with the GTP effect.

## 3 The GTP, GLWP and The perihelion advance

The Thomas precession is purely kinematical in origin . If a component of acceleration $`(\stackrel{}{a})`$ exists perpendicular to the velocity $`\stackrel{}{v}`$, for whatever reason, then there is a Thomas Precession, independent of other effects . When the acceleration is caused by a gravitational force field, the corresponding Thomas Precession is reasonably referred to as the Gravitational Thomas Precession (GTP). Given the physics involved in the Thomas Precession, the possibility of the existence of the GTP in planetary motion can not be ruled out in principle. The Thomas Precession frequency $`\stackrel{}{\omega }_T`$ in the non-relativistic limit (i.e., when $`v<<c`$) is given by

$$\stackrel{}{\omega }_T=\frac{1}{2c^2}(\stackrel{}{a}\times \stackrel{}{v}),$$
(22)
where the symbols have there usual meanings. For a planet (say Mercury) moving around the Sun, the acceleration $`\stackrel{}{a}`$ is predominately caused by the Newtonian gravitational field of the Sun,viz.,

$$\stackrel{}{a}=\frac{GM_{}}{r^3}\stackrel{}{r},$$
(23)
where the symbols have their usual meanings. Thus, from Eqs.$`(22)`$ and $`(23)`$ we get the GTP frequency of the planet in question as

$$\stackrel{}{\omega }_{gT}=\frac{GM_{}}{2c^2r^3}(\stackrel{}{r}\times \stackrel{}{v}),$$
(24)
where $`\stackrel{}{v}`$ is the velocity of the planet. In terms of the angular momentum of the planet $`\stackrel{}{L}=m_0(\stackrel{}{r}\times \stackrel{}{v})`$, Eq.$`(24)`$ can be re-written as

$$\stackrel{}{\omega }_{gT}=\frac{GM_{}}{2m_0c^2r^3}\stackrel{}{L},$$
(25)
If, as Thomas first pointed out, that coordinate system rotates, then the total time rate of change of the angular momentum $`\stackrel{}{J}`$For Thomas $`\stackrel{}{J}=\stackrel{}{S}`$, the spin angular momentum; but here we consider a more general term $`\stackrel{}{J}=\stackrel{}{L}+\stackrel{}{S}`$, $`\stackrel{}{J}`$ representing the total ( orbital + spin ) angular momentum of the particle under consideration.  or more generally, any vector $`\stackrel{}{A}`$ is given by the well known result ,

$$\left(\frac{d\stackrel{}{A}}{dt}\right)_{\text{nonrot}}=\left(\frac{d\stackrel{}{A}}{dt}\right)_{\text{rest\hspace{0.17em}\hspace{0.17em}frame}}+\stackrel{}{\omega }_T\times \stackrel{}{A}$$
(26)
where $`\stackrel{}{\omega }_T`$ is the angular velocity of rotation found by Thomas. When applied to the total angular momentum $`\stackrel{}{J}`$, Eq.$`(26)`$ gives an equation of motion:

$$\left(\frac{d\stackrel{}{J}}{dt}\right)_{\text{non-rot}}=\left(\frac{d\stackrel{}{J}}{dt}\right)_{\text{rest\hspace{0.17em}\hspace{0.17em}frame}}+\stackrel{}{\omega }_T\times \stackrel{}{J}$$
(27)
The corresponding energy of interaction is

$$U=U_0+\stackrel{}{J}\stackrel{}{\omega }_T=U_0+\stackrel{}{L}\stackrel{}{\omega }_T+\stackrel{}{S}\stackrel{}{\omega }_T$$
(28)
where $`U_0`$ is the energy corresponding to the coupling of $`\stackrel{}{J}`$ to the external fields - say the Coulomb field in atomic case ,nuclear field in nuclear case and the Newtonian gravitational field in the planetary case. The origin of the Thomas precessional frequency $`\stackrel{}{\omega }_T`$ is the acceleration experienced by the particle as it moves under the action of external forces. Since the nature of the external forces is not specified, the result obtained in Eq.(28) is valid for all type of force fields which cause accelerations of whatever nature. When applied to the gravitodynamic problems in solar system where the acceleration of a planet with respect to the Sun is predominately caused by a force arising out of the Newtonian scalar potential, Eq.(28) takes the form

$$U_g=U_{0g}+\stackrel{}{J}\stackrel{}{\omega }_T=U_{0g}+\stackrel{}{L}\stackrel{}{\omega }_{gT}+\stackrel{}{S}\stackrel{}{\omega }_{gT}$$
(29)
where $`U_{0g}`$ is the Newtonian potential energy of the planet under consideration and $`\stackrel{}{\omega }_{gT}`$ is given by Eq.(25).We then have

$$U_g=\frac{k}{r}\frac{h_1}{r^3}\frac{h_2}{r^3}$$
(30)
where $`k=GM_{}m_0`$ and

$$h_1=\frac{GM_{}L^2}{2m_0c^2}$$
(31)
$$h_2=\frac{GM_{}}{2m_0c^2}(\stackrel{}{L}\stackrel{}{S}).$$
(32)
Thus we see the gravitational Thomas precession in the non-relativistic limit introduced two potentials of the form $`\mathrm{\hspace{0.17em}1}/r^3`$ into the classical Kepler problem. If in place of the $`U_{0g}`$ in Eq.(29) we take the gravitational Lienard-Wiechert (LW) potential energy $`U_{gLW}`$ given by Eq.(20), we would then have the following equation in place of Eq.(30) :

$$U_g=\frac{k}{r}\frac{h_0}{r^3}\frac{h_1}{r^3}\frac{h_2}{r^3}=\frac{k}{r}\frac{h_{MG}}{r^3}$$
(33)
where

$$h_{MG}=h_0+h_1+h_2=\frac{GM_{}L^2}{m_0c^2}\left[\mathrm{\hspace{0.17em}1}+\frac{(\stackrel{}{L}\stackrel{}{S})}{2L^2}\right].$$
(34)
What effect will result from the introduction of the potential(s) of the form $`\mathrm{\hspace{0.17em}1}/r^3`$ into the Kepler problem ? It is shown in  that if a potential with $`\mathrm{\hspace{0.17em}1}/r^3`$ form is added to a central force perturbation of the bound Kepler problem, the orbit in the bound problem is an ellipse in a rotating coordinate system. In effect the ellipse rotates, and the periapsis appears to precess. If the perturbation Hamiltonian is

$$H=\frac{h}{r^3},(\text{ h\hspace{0.17em}=\hspace{0.17em}some constant })$$
(35)
then it predicts  a precession of the perihelion of a planet arising out of the perturbation Hamiltonian (of the form as in Eq.(35)) at an average rate of

$$\dot{\stackrel{~}{\omega }}=\frac{6\pi km_0^2h}{\tau L^4}$$
(36)
where $`k=GM_{}m_0`$ and $`\tau `$ is the classical period of revolution of the planet around the sun. It is worth-noting from  that the so-called Schwarzschild spherically symmetric solution of the Einstein field equations corresponds to an additional Hamiltonian in the Kepler problem of the form of Eq.(35) with

$$h=h_E=\frac{GM_{}L^2}{m_0c^2}$$
(37)
so that Eq.(36) becomes

$$\dot{\stackrel{~}{\omega }}_E=\frac{6\pi k^2}{\tau L^2c^2}=\frac{6\pi GM_{}}{\tau c^2a(\mathrm{\hspace{0.17em}1}e^2)}$$
(38)
where we have used the relation $`L^2=GM_{}m_0^2a(\mathrm{\hspace{0.17em}1}e^2)`$. Eq.(38) represents Einstein’s expression for the anomalous perihelion advance of a planet’s orbit. Likewise the contributions to the perihelion advance arising out of the GLWP and GTP in the framework of Maxwellian Gravity can be estimated by taking the $`h`$ in Eq.(36) as

$$h=h_{MG}=\frac{GM_{}L^2}{m_0c^2}\left[\mathrm{\hspace{0.17em}1}+\frac{(\stackrel{}{L}\stackrel{}{S})}{2L^2}\right]=h_E\left[\mathrm{\hspace{0.17em}1}+\frac{(\stackrel{}{L}\stackrel{}{S})}{2L^2}\right].$$
(39)
Then Maxwellian Gravity can predict the relativistic perihelion advance of a planet at

$$\dot{\stackrel{~}{\omega }}_{MG}=\dot{\stackrel{~}{\omega }}_E\left[\mathrm{\hspace{0.17em}1}+\frac{(\stackrel{}{L}\stackrel{}{S})}{2L^2}\right].$$
(40)
For Mercury, the value of $`\dot{\stackrel{~}{\omega }}_E=\mathrm{\hspace{0.17em}42}98`$ arc-seconds/century - a well known data . Hence the relativistic perihelion advance of Mercury’s orbit in the flat space-time Maxwellian Gravity could be predicted at

$$\dot{\stackrel{~}{\omega }}_{MG}=\mathrm{\hspace{0.17em}42}98\left[\mathrm{\hspace{0.17em}1}+\frac{(\stackrel{}{L}\stackrel{}{S})}{2L^2}\right]\text{arc-seconds/century}$$
(41)
The additional $`LS`$ term viz., $`\frac{(\stackrel{}{L}\stackrel{}{S})}{2L^2}`$, that appears in Eq.(41) has a numerical value of the order of $`\mathrm{\hspace{0.17em}10}^{10}`$ when the physical and orbital parameters of Mercury are used and is therefore utterly negligible. So by neglecting this $`LS`$ term we get

$$\dot{\stackrel{~}{\omega }}_{MG}=\dot{\stackrel{~}{\omega }}_E=\mathrm{\hspace{0.17em}\hspace{0.17em}42}98\text{arc-seconds/century \hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}\hspace{0.17em}Mercury}.$$
(42)
## 4 Concluding Remarks

In this work we saw the possibility of explaining the observed anomalous advance of the perihelion of Mercury’s orbit in flat space-time relativistic gravity. This new approach to the old gravitodynamic problem of perihelion advance may serve as a test of the validity of special relativity in the domain of gravitation. Again this also implies a test of the physics of Thomas Precession in gravitational phenomena. It is to be noted that forces and accelerations (of whatever origin) are well within the scope of special relativity (SR) because the SR in its entirity no where forbids one to study the force of gravity within its versatile scope. In this connection we would like to quote an important observation made by Denisov and Logunov:

> “ …,it must be noted that the literature not infrequently contains statements claiming that the special theory of relativity deals with the description of phenomena in the inertial reference frames, while the description of phenomena in non-inertial reference frames is the prerogative of the GTR.
> These statements are wrong. …… . Because of this, it is quite conceivable to describe the physical phenomena either by the special theory of relativity or within non-inertial reference frames. This point was trasparently clear to Fock . ”

By the way we remark that we are not proposing a new theory of gravity as MG. In our work on MG we only investigated some unexplored aspects of relativistic gravity in flat space-time and eleveted the status of Maxwell-Heaviside’s gravity to that of a test theory for testing the foundations of both special and general relativity. It is to be carefully noted that MG is now a natural outcome of some well established principles,theories and methods of study in physics. Therefore the predictions of the MG may not be totally false. The theory might be working somewhere in some domain of physics yet unexplored. So we have to explore the situations where and when the MG was/is/might be operating in the evolution of the physical world. The authors make an appeal to the redears not to consider the MG as an alternative theory of gravity to the GR, because MG has to be made compatible with many other experinental data or observational results for its elevation to that status. So we now prefer MG to be treated as a toy model vector theory of gravity in flat space-time.

Acknowledgments The authors wish to thank Prof. Lewis Ryder, Phys. Laboratory, University of Kent, Canterbury, UK; Prof. B. Mashhoon, Dept. of Phys. and Astronomy, Univ. of Missouri-Columbia, Columbia, Missouri, USA and Prof. M. Gasperini, Dept. of Phys., Univ. of Bari, Bari, Italy for their careful reading of the basic paper  and for their kindly e-mailed valuable critical remarks and constructive criticisms which we found encouraging and of some use for this work. The first author would like to thank Prof. N. Barik and Prof. L. P. Singh, both of Dept. of Phys., Utkal University, Vani Vihar, Bhubaneswar, Orissa, for their kind co-operation in some thought provoking and insightful discussions on the subject with them. The authors also acknowledge the help received from the Institute of Physics, Bhubaneswar for using its library and computer Centre for this work.