text:
# Ion–beam induced current in high–resistance materials

## I Introduction

The irradiation of semiconductors and insulators by ion beams is a tool that is frequently used for studying or modifying the electric properties of these high-resistance materials (see e.g.  and references therein). As s result of ion irradiation, heavily doped charged layers can be formed in a sample , with the density of ions in the layer reaching $`10^{20}`$cm<sup>-3</sup>. The formation of strongly nonuniform charge distributions is the standard situation in the ion-beam irradiation, while for achieving a uniform volume concentration of implanted ions special tricks are required .

Under the influence of ion irradiation, electric properties of high-resistance materials can be essentially changed , with additional peculiarities resulting from a nonuniform distribution of charge carriers \[6–9\]. High-resistance materials mean, as usual, the materials with a very poor concentration of carriers, such as semiconductors, semimetals, or insulators that, being irradiated, acquire conducting properties . The electron transport in such materials can be well described in the semiclassical drift-diffusion approach . The existence of the implantation induced damages can be effectively taken into account by means of characteristic parameters entering the drift-diffusion equations , for instance, by specifying the resistivity or mobility .

When an energetic ion strikes a solid surface, there is a probability of electron capture resulting in the neutralization of a part of implanted ions that become electrically inactive. However, the neutralized ions can easily be activated again by irradiating the material with laser beams \[12–15\].

In the present paper, we study the electric transport in high-resistance materials with a strongly nonuniform initial distribution of charge carriers, which can be formed by ion-beam irradiation. As follows from the above discussion, it is always possible to prepare such conditions, e.g. employing laser activation \[12–15\], that the injected ions could be the principal charge carriers. For concreteness, we consider below positive ions, though this assumption is not principal and negative ions could be treated in the same footing.

## II Peculiarities of Electric Transport

The transport properties of high-resistance materials, such as extrinsic semiconductors, are usually described by the drift-diffusion approach  consisting of the continuity and Poisson equations, respectively,

$$\frac{\rho }{t}+\stackrel{}{}\stackrel{}{j}+\frac{\rho }{\tau }=0,\epsilon \stackrel{}{}\stackrel{}{E}=4\pi \rho ,$$
(1)
where $`\rho =\rho (\stackrel{}{r},t)`$ is the charge density of carriers; $`\stackrel{}{j}=j(\stackrel{}{r},t)`$ is the electric current density; $`\tau `$ is the relaxation time; and $`\epsilon `$ is the dielectric permittivity. The total current density is

$$\stackrel{}{j}_{tot}=\stackrel{}{j}+\frac{\epsilon }{4\pi }\frac{\stackrel{}{E}}{t},\stackrel{}{j}=\mu \rho \stackrel{}{E}D\stackrel{}{}\rho ,$$
(2)
where $`\mu `$ is the mobility of carriers and $`D=(\mu /\epsilon )k_BT`$ is the diffusion coefficient. The average current through the considered sample is given by the integral

$$\stackrel{}{J}(t)=\frac{1}{V}\stackrel{}{j}_{tot}(\stackrel{}{r},t)𝑑\stackrel{}{r}$$
(3)
over the sample volume $`V`$.

Let us consider a plane sample of the thickness $`L`$ and area $`A`$, which is biased with an external constant voltage $`V_0>0`$. Because of the plane symmetry, the three-dimensional picture will be reduced everywhere below to the one-dimensional description. For what follows, it is convenient to simplify the notation passing to dimensionless quantities. The return to dimensional quantities can be easily done as follows. We keep in mind that the coordinate $`x`$ is measured in units of the thickness $`L`$ and the time $`t`$, in units of the transit time $`\tau _0`$, so that for returning to the corresponding dimensional variables one has to make the substitution

$$x\frac{x}{L},t\frac{t}{\tau _0},\tau _0\frac{L^2}{\mu V_0}.$$

For other physical quantities the return to dimensional units is done by accomplishing the following transitions: for the diffusion coefficient,

$$D\frac{D}{D_0},D_0\mu V_0,$$

for the electric field,

$$E\frac{E}{E_0},E_0\frac{V_0}{L},$$

for the total accumulated charge,

$$Q\frac{Q}{Q_0},Q_0\epsilon AE_0=\frac{\epsilon AV_0}{L},$$

for the charge density,

$$\rho \frac{\rho }{\rho _0},\rho _0\frac{Q_0}{AL}=\frac{\epsilon V_0}{L^2},$$

and for the electric current,

$$j\frac{j}{j_0},j_0\frac{Q_0}{A\tau _0}=\frac{\epsilon V_0}{L\tau _0},$$

with the same transformation for the average current (3), $`JJ/j_0`$.

Employing the dimensionless notation, for the plane case considered, we have from equations (1)

$$\frac{\rho }{t}+\frac{}{x}(\rho E)D\frac{^2\rho }{x^2}+\frac{\rho }{\tau }=0,\frac{E}{x}=4\pi \rho ,$$
(4)
where the space and time variables are such that

$$0<x<1,t>\mathrm{\hspace{0.33em}0}.$$

The total current density (2) reads

$$j_{tot}=\rho ED\frac{\rho }{x}+\frac{1}{4\pi }\frac{E}{t}.$$
(5)
The condition that the sample is biased with a constant external voltage now writes

$$_0^1E(x,t)𝑑x=1.$$
(6)
Note that the parameters $`\epsilon `$ and $`\mu `$ do not appear in Eqs. (4) because of the usage of the dimensionless units.

An initial condition to the continuity equation is defined by the distribution of ions after the irradiation process. The distribution of implanted species can be modelled by the Gaussian form

$$\rho (x,0)=\frac{Q}{Z}\mathrm{exp}\left\{\frac{(xa)^2}{2b}\right\},$$
(7)
in which $`0<a<1`$ and

$$Q=_0^1\rho (x,0)𝑑x,Z=_0^1\mathrm{exp}\left\{\frac{(xa)^2}{2b}\right\}𝑑x.$$

The distribution centre, $`a`$, is close to the mean projected range of ions, although may be not exactly coinciding with it.

Our aim is to study the behaviour of the electric current through the sample,

$$J(t)=_0^1j_{tot}(x,t)𝑑x,$$
(8)
as a function of time, when the initial distribution of carriers is given by the form (7). In particular, we shall show that the nonuniformity of the initial distribution may lead to quite unexpected behaviour of the current (8) when it turns against the applied voltage becoming negative, This will be shown to be the result of the solution of the transport equations (4).

First of all, let us emphasize that the occurrence of the negative current, if any, can be due only to a nonuniform distribution of carriers. Really, the current (8), with the current density (5), can be written as

$$J(t)=_0^1\rho (x,t)E(x,t)𝑑x+D\left[\rho (0,t)\rho (1,t)\right].$$
(9)
If $`\rho (x,t)`$ does not depend on $`x`$, then Eq. (9), with the help of the voltage condition (6), immediately shows that the current is positive since we are dealing with a positive charge density $`\rho `$. As far as the diffusion and relaxation terms in the continuity equation tend to smooth the initial nonuniform distribution of carriers, this tells us that the negative current can happen only at the initial stage of the process while the charge density $`\rho (x,t)`$ is yet essentially nonuniform. The initial stage means that $`t1`$, when the influence of the diffusion and relaxation terms is yet negligible. In other words, the favouring conditions for the occurrence of the negative electric current are $`D1`$ and $`\tau 1`$. It is worth recalling here that, according to the discussion given in the Introduction, the influence of the implantation-induced damages is assumed to be taken into account by the corresponding values of parameters in the transport equations and that the implanted ions are assumed to be prepared as electrically active, which can be done by means of laser irradiation \[12–15\].

To demonstrate explicitly that the negative current really can occur, let us analyse the case of a very narrow layer of injected ions, such that the standard deviation $`b`$, also called the straggling , is small, $`b1`$. Then distribution (7) is close to

$$\rho (x,0)=Q\delta (xa).$$
(10)
In the beginning of the process, when $`t0`$, we have from Eqs. (9) and (10)

$$J(0)=QE(a,0).$$

Integrating the second equation in Eq. (4) yields

$$E(a,0)=1+4\pi Q\left(a\frac{1}{2}\right).$$

With these conditions, the current (8) becomes negative provided that

$$4\pi Q\left(\frac{1}{2}a\right)>\mathrm{\hspace{0.33em}1}.$$
(11)
From here, the inequality

$$a<\frac{1}{2}\frac{1}{4\pi Q}$$
(12)
follows for the location of the ion layer. Since this location is inside the sample, we also have the inequality

$$Q>\frac{1}{2\pi }(0<a<1)$$
(13)
for the initial charge.

The above analysis demonstrates that the negative current really can occur provided some special conditions, as (12) and (13), hold true. To our knowledge, there have been no experimental measurements demonstrating the appearance of such a negative current. Therefore the picture we describe here is a theoretical prediction of a novel effect. The occurrence of this negative current is very sensitive to the characteristics of the irradiated sample as well as to the initial nonuniform distribution of carriers, which suggests the possibility of using this effect for studying the mentioned properties. For instance, the mean projected range of irradiating ions could be measured in this way. This can be realized as follows. Assume that for an ion-irradiated material we observed the occurrence of the negative current $`J(0)`$ at initial time. Let us compare the observed value $`J(0)`$ with formula (9) where we have to substitute the distribution (7). The peak of the latter is close to the mean projected range of the irradiating ions. With the given $`\rho (x,0)`$, the right-hand side of expression (9) can be explicitly calculated. This is because the electric field satisfying Eq. (4), with the voltage condition (6), can be presented as the functional

$$E(x,t)=1+4\pi \left[Q(x,t)_0^1Q(x,t)𝑑x\right]$$

of the density $`\rho (x,t)`$ entering through

$$Q(x,t)=_0^x\rho (x^{},t)𝑑x^{}.$$

Hence, for a given $`\rho (x,0)`$, the electric field $`E(x,0)`$ is uniquely defined by the above functional. Thus, equating the calculated $`J(0)`$ from formula (9) with the corresponding measured value, we obtain an equation for the dimensionless mean projected range $`a`$. For example, in the case of a narrow ion layer, we find

$$a=\frac{1}{2}\frac{1}{4\pi Q}+\frac{J(0)}{4\pi Q^2}.$$

Returning to dimensional units, for the mean projected range $`\lambda =aL`$ we obtain

$$\lambda =\frac{L}{2}\frac{\epsilon AV_0}{4\pi Q}+\frac{\epsilon A^2L^2}{4\pi \mu Q^2}J(0).$$
(14)
This formula directly relates the mean projected range of ions, $`\lambda `$, with the known parameters of the system and the measured initial current $`J(0)`$. In deriving Eq. (14), we have only used the fact that the initial distribution of ions is narrow, $`b1`$, while the current $`J(0)`$ can, in general, be of any sign. The advantage of using the effect of the negative electric current is that the latter always requires a relatively narrow initial distribution. Therefore, as soon as we observe a negative current $`J(0)`$, we can employ formula (14) as a realistic approximation for the ion mean projected range.

When the initial ion layer is not narrow and the current $`J(0)`$ is not necessary negative, so that the dependence of the formulas on the standard deviation $`b`$ becomes important, then two situations can happen. One is when $`b`$ is known from other experiments. Then, in order to find the mean projected range, one should proceed exactly as is described above equating the measured and calculated currents $`J(0)`$. If the straggling $`b`$ is not known, it can be found in the following way. Assume that at the initial time $`t=0`$ there appears the negative current. Since this is a transient effect, existing only in a finite time interval $`0tt_0`$, there is some time $`t=t_0`$ when the current changes its sign becoming normal, i.e. positive, for $`t>t_0`$. At the same time $`t=t_0`$, the current then is zero, $`J(t_0)=0`$. Therefore, the latter equation, with the given experimentally measured $`t_0`$, may serve as an additional condition defining the straggling $`b`$.

It may also happen that the injected ion layer is narrow but located close to the surface, so that $`ab`$. Then formula (14) can be corrected by taking into account the second term in Eq. (9), which gives for the mean projected range

$$\lambda =\frac{L}{2}\frac{\epsilon A}{4\pi Q}\left[V_0+\frac{AL}{\mu Q}\left(D\mathrm{\Delta }\rho LJ\right)\right],$$
(15)
where

$$JJ(0),\mathrm{\Delta }\rho \rho (0,0)\rho (1,0).$$

To find the electric current (8) for arbitrary times, we need to solve equations (4), with the voltage condition (6) and the initial condition (7). Those should also be complimented by boundary conditions which can be taken in the Neumann form . We have accomplished such numerical calculations whose results are presented in the figures, where all values are given in dimensionless units, and the conditions $`D1`$ and $`\tau 1`$, favouring the occurrence of the negative current are assumed. The figures correspond to the setup explained in detail in the Introduction and specified in the present section.

Fig. 1 shows the dependence of the electric current $`J(t)`$ on time for the varying location $`a`$ of the initial distribution peak of ions. Formulas (14) and, respectively, (15) for the mean projected range are valid with a very good accuracy, with an error being less than $`1\%`$. In order to distinguish the lines, they are numbered.

Fig. 2 presents the electric current as a function of time for different widths of the initial distribution of ions. Formulas (14) and (15) are again very accurate for the negative current $`J(0)`$. But when the current $`J(0)`$ is positive and at the same time, the initial distribution is not narrow, then Eqs. (14) and (15) are not valid, as it should be.

Fig. 3 demonstrates the time dependence of the electric current for different values of the charge $`Q`$. Since $`b1`$, formulas (14) and (15) work perfectly for both cases, for the negative current, as in Fig. 3, as well as for positive current, as in Fig. 4; the errors being not more than around $`1\%`$. The current in Fig. 4 becomes positive because condition (13) does not hold. Increasing the accumulated charge $`Q`$, so that condition (13) becomes true, the negative current appears provided that condition (12) is satisfied, as is clearly illustrated in Fig. 5.

## III Discussion

We have shown that in high resistance materials irradiated by ion beams an unusual transient effect of the negative electric current can occur. This effect can be used for studying the characteristics of irradiated materials as well as those of irradiating ions. For instance, the mean projected range of ions can be measured whose value is well approximated by formula (15). This new way of measuring the ion mean projected range is, certainly, not the sole possible, but it can provide additional information being complimented by other experimental methods. Throughout the paper, we have been talking about ions. But, as is clear, the same consideration concerns any kind of charged particles, say electrons, positrons, or muons.

It is worth stressing that here we have advanced a theoretical proposal for observing and employing a novel effect. The general prerequisites for realizing the latter look, as follows from the above discussions, achievable. Moreover, there exist so many types of various high-resistance materials, such as semiconductors, semimetals, and insulators, there are so many ways of fabricating such materials with specially designed characteristics, and also there are so many methods of varying the properties of these materials by involving additional external sources, as electromagnetic fields or laser beams, that the feasibility of preparing the desired conditions for realizing the suggested effect looks quite realistic. Even if this realization is difficult today, it might be accomplished in future. The history shows that practically any effect that is possible in principle sooner or later becomes realizable in experiment.

Although our main aim in this paper has been to advance a theoretical proposal for the principal possibility of observing a novel effect, we have also paid attention to the feasibility of realizing it experimentally. In addition to the above discussions, we would like to touch several other important points that could substantiate the feasibility of observing the considered effect.

First, the space charge due to the accumulated ions will act against further implantation, and the ion doses will certainly be limited. At the same time, to observe the effect of negative electric current, one needs to reach some threshold value for the accumulated charge, as in the inequality (13). In order to understand how the latter condition can be satisfied, it is necessary to rewrite it in the corresponding dimensional units introduced in the beginning of Sec. 2, in this case, in the units of $`Q_0`$. This gives us the condition

$$2\pi \epsilon A\frac{V_0}{L}Q>1.$$

From here, it is evident that this condition can always be satisfied for any given charge $`Q`$, which can be achieved by increasing either the sample area $`A`$ or the applied voltage $`V_0`$. Of course, for a given sample, one has to increase the voltage $`V_0`$.

The implanted ions do not form a homogeneous layer and their distribution may be not exactly Gaussian. In this paper we consider the Gaussian distribution (7) which usually describes well the profile of implanted ions. However this form is not compulsory for realizing the negative-current effect. The latter persists as well for other distributions, for which one needs only that the straggling $`b`$ be less than the projected range $`a`$.

For the realization of the effect, it is also not necessary that the ion distribution be centered at the mean projected range of ions. For simplicity, we called the distribution center $`a`$ the projected range. However, if these two do not coincide, all consideration remains valid with merely slightly changing the terminology, so that $`a`$ is to be called the distribution center.

To measure currents through the sample, one has to take into account that near contacts there often occurs the space charge build-up caused by electron or hole injection from contacts. The current due to the carriers injected from contacts can mask the current produced by ions. How to ascertain that the considered effect is caused by the implanted ions? This problem can be easily resolved by measuring the current through the sample before ion irradiation. Then one may study the influence of contacts on creating electric currents. In this way, one can always separate the influence of contacts from physical effects resulting from ion irradiation.

In conclusion, we do not see principal difficulties for realizing the negative-current effect. And, as follows from the above discussions, different technical problems seem to be resolvable. Being realized, this effect will make it possible to have an additional tool for studying the transport properties of high-resistance materials as well as to measure characteristics related to irradiated ions.

Acknowledgement

We appreciate financial support from the São Paulo State Research Foundation.

Figure Captions

Fig. 1.

Electric current through the sample as a function of time for the parameters $`Q=1,b=0.05`$, and different initial locations of the ion layer: $`a=0.1`$ (curve 1, solid line), $`a=0.25`$ (curve 2, long–dashed line), $`a=0.5`$ (curve 3, short–dashed line), and $`a=0.75`$ (curve 4, dotted line).

Fig. 2.

The dependence of the electric current on time for $`Q=0.5,a=0.1`$, and varying widths of the initial ion distribution: $`b=0.05`$ (curve 1, solid line), $`b=0.1`$ (curve 2, long–dashed line), and $`b=0.5`$ (curve 3, short–dashed line).

Fig. 3.

Electric current as a function of time for $`a=0.25,b=0.05`$ and different charges: $`Q=3`$ (curve 1, solid line), $`Q=1`$ (curve 2, long–dashed line), and $`Q=0.5`$ (curve 3, short–dashed line).

Fig. 4.

The time dependence of the electric current for $`Q=0.1,b=0.05`$, and two different initial locations of the ion distribution: $`a=0.1`$ (curve 1, solid line) and $`a=0.25`$ (curve 2, long–dashed line).

Fig. 5.

Electric current vs. time for $`Q=0.5,b=0.1`$, and varying initial locations of the peak of the ion distribution: $`a=0.1`$ (curve 1, solid line), $`a=0.25`$ (curve 2, long–dashed line), $`a=0.5`$ (curve 3, short–dashed line), and $`a=0.75`$ (curve 4, dotted line).


format: text

text:
# Can we predict the fate of the Universe?

## I Introduction

The issue of the present state, future dynamics and final fate of the Universe, or at least our patch of it, has been recently pushed to the front line of research in cosmology. This is mostly due to observations of high redshift type Ia supernovae, performed by two independent groups (the “Supernova Cosmology Project” and the “High-Z Supernova Team”), which allowed accurate measurements of the luminosity-redshift relation out to redshifts up to about $`z1`$ . It should be kept in mind that these measurements are done on the assumption that these supernovae are standard candles, which is by no means demonstrated and could conceivably be wrong. There are concerns about the evolution of these objects and the possible dimming caused by intergalactic dust , but we will ignore these for the purposes of this paper, and assume that the quoted results are correct.

The supernovae data, when combined with the ever growing set of CMBR anisotropy observations, strongly suggest an accelerated expansion of the Universe at the present epoch, with cosmological parameters $`\mathrm{\Omega }_\mathrm{\Lambda }0.7`$ and $`\mathrm{\Omega }_\mathrm{m}0.3`$. A further cause of concern here is the model dependence of the CMBR analysis, but we shall again accept the above results for the purpose of this paper.

Taken at face value, these results would seem to show that the universe will necessarily enter an inflationary stage in the near future. However, as pointed out by Starkman, Trodden and Vachaspati  this is not necessarily so. We could be living in a small, sub-horizon bubble, for example. And even if we were indeed inflating, it would not be trivial to demonstrate it. In the above work, these authors looked at the crucial question of ‘How far out must we look to infer that the patch of the universe in which we are living is inflating?’ Their analysis is based on previous work by Vachaspati and Trodden  which shows that the onset of inflation can in some sense be identified with the comoving contraction of our minimal anti-trapped surface (MAS)***The MAS of each comoving observer is a sphere centered on him/her, on which the velocity of comoving objects is $`c`$. For the particular case of an homogeneous universe, the MAS has a physical radius $`cH^1`$.. They then argue that if one can confirm cosmic acceleration up to a redshift $`z_{MAS}`$ and detect the contraction of our MAS, then our universe must be inflating. Unfortunately, even if we can do the former (for $`\mathrm{\Omega }_\mathrm{\Lambda }=0.8`$ the required redshift is $`z_{MAS}1.8`$), it turns out that there is no way to presently confirm the latter, because the accelerated expansion hasn’t been going for long enough for the MAS to contract. Only if we had $`\mathrm{\Omega }_\mathrm{\Lambda }0.96`$ would we be able to demonstrate inflation today.

As in the proverbial mathematicians joke, the method outlined by Starkman et al. provides an answer that is completely accurate but will take a long time to find, and hence is of no immediate use to us. In this paper, however, we will explore a different possibility. Our main aim is to provide what could be called a physicists version of the “mathematical” question of Starkman et al. . In other words, we are asking, ‘If we can’t know for sure the fate of the universe at present, what is our best guess today?’. As we will discuss, we can answer this question, although it will involve making some crucial additional assumptions.

In order to answer the above question we compare the particle and event horizons. We show that for a flat universe with $`\mathrm{\Omega }_\mathrm{\Lambda }\text{ }>0.14`$ the particle horizon is greater than the distance to the event horizon meaning that today we may be able to observe a larger portion of the Universe than that which will ever be able to influence us. We argue that if we find evidence for a constant vacuum density up to a distance from us equal to the event horizon then our Universe will necessarily enter an inflationary phase in the not too distante future, assuming that the potential of the scalar field which drives inflation is time-independent and that the content of the observable universe will remain ‘frozen’ in comoving coordinates.

Note that Starkman et al. argue that inflation can only take place if the vacuum energy dominates the energy density on a region with physical radius not smaller than that of the MAS at that time. However, they did not assume that the content of the observable universe would remain frozen in comoving coordinates and so they found that the larger is the contribution of a cosmological constant for the total density of the Universe, the larger is the redshift out to which one has to look in order to infer that our portion of the universe is inflating. This result seems paradoxical, until one realizes that the size of the MAS at a given time does not, by itself, say anything about inflation. The main reason why $`z_{MAS}`$ grows with $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is simply because the scale factor has grown more.

The plan of this paper is as follows. In the next section we introduce the various lengthscales that are relevant to our discussion, and provide a qualitative discussion of our test for inflation. We also discuss the assumptions involved and compare our ‘physical’ test with the ‘mathematical’ one recently proposed by Starkman et al. . In section III we provide a more quantitave analysis of our criterion. We also discuss in more detail our crucial assumption of an energy-momentum distribution which remains frozen in comoving coordinates. Finally, in section IV we summarize our results and discuss some other outstanding issues.

## II A ‘physical’ test for inflation

The dynamical equation which describes the evolution of the scale factor $`a`$ in a Friedmann-Robertson-Walker (FRW) universe containing matter, radiation and a cosmological constant can be written as

$$H^2=H_0^2(\mathrm{\Omega }_{m0}a^3+\mathrm{\Omega }_{r0}a^4+\mathrm{\Omega }_{\mathrm{\Lambda }0}+\mathrm{\Omega }_{k0}a^2).$$
(1)
where $`H=\dot{a}/a`$ and the density parameters $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_r`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ express respectively the densities in matter, radiation and cosmological constant as fractions of the critical densityA dot represents a derivative with respect to the cosmic time $`t`$. The subscript ‘$`0`$’ means that the quantities are to be evaluated at present epoch, and we have also taken $`a_0=1`$.. Naturally one has $`\mathrm{\Omega }_k=1\mathrm{\Omega }_m\mathrm{\Omega }_r\mathrm{\Omega }_\mathrm{\Lambda }`$.

The distance $`d`$, to a comoving observer at a redshift $`z`$ is given by

$$d(z)=c_{t(z)}^{t_0}\frac{dt^{}}{a(t^{})}=cH_0^1_0^z\frac{dz^{}}{\left[\mathrm{\Omega }_{m0}(1+z^{})^3+\mathrm{\Omega }_{r0}(1+z^{})^4+\mathrm{\Omega }_{\mathrm{\Lambda }0}+\mathrm{\Omega }_{k0}(1+z^{})^2\right]^{1/2}},$$
(2)
and is related to the ‘radius’ of the local universe which we can in principle observe today. The distance to the event horizon can be defined as

$$d_e=c_{t_0}^{\mathrm{}}\frac{dt^{}}{a(t^{})}=cH_0^1_1^{\mathrm{}}\frac{da}{(\mathrm{\Omega }_{m0}a+\mathrm{\Omega }_{r0}+\mathrm{\Omega }_{\mathrm{\Lambda }0}a^4+\mathrm{\Omega }_{k0}a^2)^{1/2}}.$$
(3)
and represents the portion of the Universe which will ever be able to influence usIn writing the upper integration limit as infinity we are of course assuming that the universe will keep expanding forever; an analogous formal definition could be given for an universe ending in a ‘big crunch’.. On the other hand, the particle horizon, $`d_p`$, is defined by (from eqn. (2))

$$d_p\underset{z\mathrm{}}{lim}d(z),$$
(4)
and it represents the maximum distance which we can observe today.

If today the distance to the event horizon is smaller than the particle horizon ($`d_e<d_p`$) this means that today we are able to observe a larger portion of the Universe than that which will ever be able to influence us. We can do this if we look at a redshift greater than $`z_{}`$ defined by (see also )

$$d(z_{})=d_e.$$
(5)
In a flat universe solutions to this equation are only possible for $`z_{}1`$ and for $`\mathrm{\Omega }_{\mathrm{\Lambda }0}\text{ }>0.14`$. Hence, assuming that the energy-momentum distribution within the patch of the Universe which we are able to see remains unchanged in comoving coordinates, our Universe will necessarily enter an inflationary phase in the future if there is a uniform vacuum density permeating the Universe up to a redshift $`z_{}`$.

This assumption obviously requires some further discussion. One can certainly think of a universe made up of different ‘domains’, each with its own values of the matter and vacuum energy density. Furthermore, by cleverly choosing the field dynamics, one can always get patches with time-varying vacuum energy densities, or patches where the vacuum energy density is non-zero for only short periods. In all such cases, the domain walls separating these patches can certainly have a very complicated dynamics, and in particular it is always possible that a domain wall will suddenly get inside our horizon sometime between the epoch corresponding to our observations and the present day. On the other hand, it should also be pointed out that a certain amount of fine-tuning would be required to have a bubble coming inside our horizon right after we have last observed it. In these circumstances, the best that can be done is to impose constraints on the characteristics of any bubble wall that could plausibly have entered the patch of the universe we are currently able to observe, given that we have so far seen none. We shall analyse this point in a more quantitative manner in the following section.

We think that the results obtained in this way, even if less robust from a formal point of view, are intuitively more meaningful than those obtained in  in the sense that, among other things, in this case the minimum redshift $`z_{}`$ out to which one must observe in order to be able to predict an inflationary phase (subject to the conditions mentioned earlier) decreases as $`\mathrm{\Omega }_\mathrm{\Lambda }`$ increases—see Fig. 1. In other words, the larger the present value of the cosmological constant, the easier it should be to notice it.

It is perhaps instructive to compare our test with that of  in more detail. Starkman et al. require the contraction of the MAS. Now, in order to see the MAS one has to look at a redshift defined by:

$$a(z_{MAS})d(z_{MAS})cH^1(z_{MAS})=0.$$
(6)
This finds the redshift, $`z_{MAS}`$, for which the physical distance to a comoving observer at that redshift, evaluated at the corresponding time $`t_{MAS}`$, is equal to the Hubble radius at that time. However, if the vacuum density already dominates the dynamics of the Universe at the redshift $`z_{}`$ then eqn. (5) reduces to:

$$d(z_{})cH^1(z_{})=0.$$
(7)
(recall that $`a_0=1`$) because during inflation the physical size of the event horizon is simply equal§§§This is only exactly true when the vacuum energy density is the only contributor to the energy density, in which case exponential inflation occurs. to $`cH^1`$ (this is ultimately the reason for the choice of criterium for inflation by Vachaspati and Trodden ). As has been discussed above, eqns. (6) and (7) have totally different solutions ($`z_{}=1`$ while $`z_{MAS}\mathrm{}`$).

To put it in another way, the main difference between our approach and that of ref. lies on the fact that we assume that the energy-momentum content of the observable Universe does not change significatively in comoving coordinates. This allows us to use the equation of state of the local universe observed for a redshift $`z`$ (looking back at a physical time $`t(z)`$) to infer the equation of state of the local universe at the present time. In the following section, we shall discuss these points in somewhat more detail.

We should also point out that if we were to relax the assumption of a co-movingly frozen content of the observable universe, then the equation—analogous to (5)—specifying the redshift out to which one should look in order to be able to predict the future of the Universe would be

$$d(z_+)=d_e(z_+).$$
(8)
This equation has no solution, so a stronger test of this kind is not feasible in practice.

## III Discussion

Here we go through some specific aspects of our test in more quantitative detail. To begin with, we have solved numerically eqn. (5); the numerical results were obtained for choices of cosmological parameters such that $`\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }=0.7,1.0,1.3`$, with an additional $`\mathrm{\Omega }_m=0.3`$ for illustration. We are interested only in a matter–dominated or $`\mathrm{\Lambda }`$dominated epoch of the evolution of the universe, and therefore we have dropped the radiation density parameter $`\mathrm{\Omega }_0^\mathrm{r}`$ of eqn. (1), in the calculations.

These results are displayed in Fig. 1 as a function of $`\mathrm{\Omega }_{\mathrm{\Lambda }0}`$. The cases with constant total density are shown in solid curves (with the top curve corresponding to the higher value of the density), while the case of a fixed $`\mathrm{\Omega }_\mathrm{m}`$, is shown, for comparison purposes, by a dotted curve.

As expected, as the universe becomes more $`\mathrm{\Lambda }`$dominated and/or less matter–dominated, the comoving distance to the event horizon decreases, which is reflected in the decrease of the redshift $`z_{}`$ of a comoving source located at that distanceNote that pushing $`\mathrm{\Omega }_{\mathrm{\Lambda }0}`$ down to zero, the value of $`z_{}`$ tends to infinity, since in such universes an event horizon does not exist.. For the observationally preferred values of $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, the required redshift will be $`z_{}1.8`$.

For comparison, the redshift, $`z_{MAS}`$, defined by eqn. (6), which is the analogous relevant quantity for the criterion of Starkman et al. is shown, for the same choices of cosmological parameters, in Fig. 2. Note that in this case, as the universe becomes more $`\mathrm{\Lambda }`$dominated and/or less matter–dominated, the redshift of the MAS will increase. As we already pointed out our test will not be applicable for very low values of the vacuum energy density, and for intermediate values, it requires a higher redshift than $`z_{MAS}`$. However, for high values of the cosmological constant and/or low matter contents, the fact that the universe will be expanding much faster makes the redshift of the MAS increase significantly, and even become larger than $`z_{}`$ for some combinations of cosmological parameters. For the same observationally preferred values of $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ quoted above, the required redshift will be $`z_{MAS}1.6`$. A comparison of the values of $`z_{}`$ and $`z_{MAS}`$ for the spatially flat model is shown in Fig. 3.

We now return to our assumption about the energy-momentum content of the universe, considering the possibility that different regions of space may have different values for the vacuum energy density, which are separated by domain walls. This means that we are assuming the existence of a scalar field, say $`\varphi `$, which within each region sits in one of a number of possible minima of a time-independent potential. It is obvious that if the potential depends on time or if the scalar field did not have time to roll to the minimum of the potential, then it is not possible to predict the fate of the universe without knowing more about the particle physics model which determines its dynamics. For simplicity, we shall assume that we live in a spherical domain with constant vacuum energy density (effectively a cosmological constant) that is surrounded by a much larger region in which the the vacuum density has a different value—for the present purposes we will assume it to be zero. Note that this is the case where the dynamics of the wall will be faster (more on this below).

Is it possible that a region with a radius $`d(z)`$, say centred on a nearby observer, can be inside a given domain at the conformal time $`\eta (z)`$, but outside that domain at the present time, $`\eta _0`$? This problem can provide some measure of how good the assumption of a frozen energy-momentum distribution in comoving coordinates is. In other words, is it likely that a domain wall may have entered this region at a redshift smaller than $`z`$? In order to provide a more quantitative answer to this question, we have performed numerical simulations of domain wall evolution using the PRS  algorithm, in which the thickness of the domain walls remains fixed in comoving coordinates for numerical convenience. See also  for a description of the simulations.

We assume that the domain wall has spherical symmetry, thereby reducing a three-dimensional problem to a one-dimensional one. We perform simulations of this wall in a flat universe on a one-dimensional $`8192`$ grid. The comoving grid spacing is $`\mathrm{\Delta }x=c\eta _i`$ where $`\eta _i=1`$ is the conformal time at the beginning of the simulation. The initial comoving radius of the spherical domain was chosen to be $`R=2048\mathrm{\Delta }x`$, and the comoving thickness of the domain wall was set to be $`10\mathrm{\Delta }x`$. In these simulations we neglect the gravitational effect induced by the different domains and domain walls on the dynamics of the universe, and we also do not consider the possibility of an open universe. We must emphasise, however, that both these effects would slow down the defects, thereby helping to justify our assumption of a constant equation of state in comoving coordinates even more.

We have obtained the following fit for the radius of the domain wall as a function of the conformal time $`\eta `$

$$R(\eta )=R_{\mathrm{}}\left(1\left(\frac{c\eta }{\alpha R_{\mathrm{}}}\right)^n\right)^{2/n},$$
(9)
where

$$\alpha =2.5,n=2.1$$
(10)
and $`R_{\mathrm{}}`$ is the initial comoving radius of the domain wall (with $`R_{\mathrm{}}c\eta _i`$). This fit is accurate to better than $`5\%`$, except for the final stages of collapse. In a flat universe with no cosmological constant the comoving distance to a comoving object at a redshift $`z`$ is given by

$$d(z)=c\eta _0\left(11/\sqrt{1+z}\right),$$
(11)
whereas the radius of the spherical domain wall can be written as a function of the redshift $`z`$, given its initial radius $`R_{\mathrm{}}`$, as follows

$$R(z)=R_{\mathrm{}}\left(1\left(\frac{c\eta _0}{\alpha R_{\mathrm{}}\sqrt{1+z}}\right)^n\right)^{2/n}.$$
(12)
Now, by solving the equation

$$R(z=0)=d(z)$$
(13)
we can find the initial comoving radius of our domain (in units of the present conformal time, that is $`R_i(z)/\eta _0`$) which it would be required to have so that its comoving size today is equal to the comoving distance to an object at a redshift $`z`$. Finally, we can calculate the radius of this domain at the present time $`\eta _0`$ and at the redshift $`z`$ (call it $`R_{\mathrm{max}}(z)`$) with the value of $`R_i(z)/\eta _0`$ obtained from the previous equation. In Fig. 4 we plot the value of $`R_{\mathrm{max}}(z)/d(z)`$ as a function of the redshift, $`z`$.

If the radius of our domain at a redshift $`z`$ was smaller than $`d(z)`$ the domain wall would be in causal contact with us at the present time and we could in principle detect the gravitational effect both of the domain wall and of the different vacuum density outside our bubble. On the other hand, if the radius of our domain at a redshift $`z`$ was greater than $`R_{max}`$ then it would not have time to enter the sphere of radius $`d(z)`$ before today.

When the redshift of the cosmological object we are looking at is small, that is ($`z0`$), its comoving distance from us, $`d(z)`$, is much smaller than the comoving horizon, $`\eta (z)`$, at the time at which the light was emitted. Consequently, a domain wall with a comoving size equal to $`d(z)`$ at the present time would already have a velocity very close to the speed of light by the redshift $`z`$. It is easy to calculate the maximum comoving size, $`R_{max}`$, which our domain would need to have at the redshift $`z`$, in order for the domain wall to enter a sphere of comoving radius $`d(z)`$ centred on a nearby observer sometime between today and redshift $`z`$. This is simply given by

$$\frac{R_{max}(z)}{d(z)}2$$
(14)
when $`z0`$, because $`d(z)`$ is the distance travelled by light from a redshift $`z`$ until today (see fig. 4).

If we assume that the comoving radius of our bubble at a redshift $`z`$ is larger than $`\eta (z)`$, then it will remain frozen in comoving coordinates until its size gets smaller than the horizon. This means that in this case the value of $`R_{max}/d(z)`$ is even smaller, approaching

$$\frac{R_{max}(z)}{d(z)}\frac{\alpha ^2\times (\sqrt{1+4\alpha ^n}1)^{2/n}}{2^{2/n}},$$
(15)
when $`z\mathrm{}`$ (see fig. 4). For a spherical domain we have

$$\frac{R_{max}(z)}{d_{\mathrm{}}}1.12.$$
(16)
We thus see that for the purposes of predicting the fate of the universe it may be a plausible assumption to assume a fixed content in comoving coordinates. The above discussion also suggests, in particular, that one may find a posteriori that it is indeed a reasonable assumption if we can observe the dynamical effects of a uniform vacuum density up to a redshift $`z1`$.

## IV Conclusions

We have provided a simple analysis of the use of cosmological observations to infer the state and fate of our patch of the universe. In particular, in the same spirit of Starkman et al.  , we have discussed possible criteria for inferring the present or future existence of an inflationary epoch in our patch of the universe.

We have presented a ‘physical’ criterion for the existence of inflation, and contrasted it with the ‘mathematical’ one that has been introduced in . Ours has the advantage of being able to provide (in principle) a definite answer at the present epoch, but the disadvantage of ultimately relying on assumptions on the content of the local universe and on field dynamics. We consider our assumptions to be plausible, but we can certainly conceive of (arguably contrived or fine-tuned) mechanisms that would be capable of violating it.

###### Acknowledgements.

C. M. is funded by FCT (Portugal) under ‘Programa PRAXIS XXI’ (grant no. PRAXIS XXI/BPD/11769/97). We thank Centro de Astrofísica da Universidade do Porto (CAUP) for the facilities provided.


text:
# Scarred Patterns in Surface Waves

## I Introduction

Parametrically forced surface waves arising as a result of the Faraday instability have provided an excellent opportunity to study nonlinear pattern formation. One of the special features of this system is that the system size relative to the basic correlation length can be varied so that both the large aspect ratio and small aspect ratio limits can be explored. At large aspect ratio, all of the classic ordered patterns have been found, including stripes, hexagons, and squares; additional exotic structures such as quasicrystalline and superlattice patterns have also been found, as well as secondary instabilities giving rise to spatiotemporal chaos. Extensive references can be found in .

The case of small aspect ratio has also been studied in rectangular and circular containers. Typically the wave patterns found near onset are either normal modes of the container or symmetrical combinations of these modes . For example, in the circular case the normal modes are Bessel functions of the radius multiplied by sinusoidal functions of the azimuthal angle. The effects of container shape can be either a nuisance or a benefit depending on one’s point of view. One example of the usefulness of considering container geometry is the study by Lane et al.  in which the conceptual differences between square symmetry and square geometry were elucidated. On the other hand, finite size effects impede efforts to utilize amplitude equations to describe the wave dynamics.

The influence of the container shape is also a fundamental issue in the field of quantum chaos. There is known to be a close correspondence between certain finite quantum systems (or analogous systems supporting classical waves) and their particle (or ray-optic) counterparts . Of particular interest are non-integrable quantum systems with classical counterparts that are chaotic, such as the billiard formed from two semicircles separated by two straight edges. For almost all initial conditions, a particle launched inside such a billiard will exhibit sensitive dependence on initial conditions - the hallmark of chaos. Experimental, numerical and theoretical studies have shown that the statistical behavior of the wavefunctions of the quantum or wave version of this system is distinctly different from the behavior for “integrable” geometries such as a square or circle . Most notably, regions of high amplitude in the wavefunctions - called scars - are found along some of the deleted word periodic orbits of the classical counterpart .

Effects of this type were explored to a limited extent using parametrically forced surface waves by Blümel et al. . Their experiment utilized water as a working fluid and high frequency excitation. They reported observations of “scarlets”, that are ridge like structures consistent with a random superposition of plane waves  (but are not located along periodic orbits.) On the other hand, no clear evidence for the simpler “scarred” wavefunctions was given. The possible effects of hydrodynamic nonlinearity on the utility of the ray optics approach has also not been discussed. Nonlinearity is in principle important, since even infinitesimally above the onset of instability, saturation of the wave amplitude is produced by nonlinear effects.

In this paper we first present observations of the spatial modes of Faraday waves in a finite non-integrable geometry, close to onset where nonlinearity is as weak as possible, and the waves might be usefully described by quasi-classical (ray optics) methods. A low viscosity fluid in a stadium shaped container is used for this purpose. Scarred patterns that resemble the computed eigenfunctions of the stadium geometry are clearly evident, but some of the linear eigenfunctions are apparently suppressed. The relative suppression of certain modes has been plausibly explained by Agam and Altshuler  in terms of higher dissipation rates for those modes near the boundaries in comparison with the modes that are observed.

We then consider the evolution of the wave patterns as the degree of nonlinearity is increased. Transitions between modes are found at some driving frequencies, along with a general increase in spatial complexity. The scars that are characteristic of the linear eigenmodes are often evident substantially above onset. Finally, we consider the development of spatiotemporal chaos (STC) in the stadium geometry. Though the onset of STC is strongly dependent on the excitation frequency, the boundaries continue to play a large role, leading for example to coherent oscillations between symmetric and asymmetric states, a phenomenon that we study using the concept of pattern entropy.

## II Theoretical background

A fluid layer with a free surface is subjected to an oscillatory vertical acceleration of amplitude $`a`$. The surface is flat until a critical acceleration $`a_c`$ is reached, at which point the surface becomes unstable and standing wave patterns are observed that oscillate at half the driving frequency. The threshold acceleration depends on the frequency and the viscosity of the fluid. It is convenient to define a dimensionless driving parameter $`ϵ=(aa_c)/a_c`$ that measures the departure from onset and hence the degree of nonlinearity. The patterns are time-independent for a range of positive $`ϵ`$ but eventually a secondary instability gives rise to spatiotemporal chaos. In this section, we briefly discuss the linear inviscid theory, the effect of viscosity, and the role of nonlinearity, as they pertain to the present investigation.

The linear stability theory for Faraday waves was first developed by Benjamin and Ursell . We summarize it here because the quantum/classical correspondence occurs for linear waves. They started from the Euler (inviscid) equation of motion and the continuity equation for an ideal fluid with a free surface in an oscillating gravitational field, and simplified the equations by retaining only the linear terms appropriate for small amplitude waves. The surface deformation $`h(𝐱,t)`$ as a function of spatial coordinate $`𝐱`$ and time $`t`$ may be written as a superposition of normal modes $`\psi _i(𝐱)`$ with coefficients $`A_i(t)`$:

$$h(𝐱,t)=\underset{i}{}A_i(t)\psi _i(𝐱),$$
(1)
where $`\psi _i(𝐱)`$ is a complete orthogonal set of eigenfunctions of the Helmholtz equation,

$$(^2+k_i^2)\psi _i(𝐱)=0.$$
(2)
The sidewall boundary condition is imposed by setting the normal component of the velocity of the fluid at the wall to zero. This leads to a quantization condition on $`k_i`$ (the wavenumber). In addition, $`k_i`$ satisfies the dispersion relation which relates the frequency of oscillation $`\omega `$ of the fluid to the wavenumber:

$$\omega ^2=\mathrm{tanh}(k_id)(\frac{\mathrm{\Gamma }}{\rho }k_i^3+gk_i),$$
(3)
where, $`\rho `$ is the fluid density, $`\mathrm{\Gamma }`$ is the surface tension, $`d`$ is the mean fluid depth, and $`g`$ is the gravitational acceleration. In our experiment, the wave number is sufficiently large so that the surface tension term is much greater than the gravity term. The hyperbolic tangent factor is close to unity since $`k_id`$ is large.

The time dependent amplitudes $`A_i`$ of these normal modes satisfy the Mathieu equation:

$$\frac{d^2A_i}{dt^2}+k_i\mathrm{tanh}k_ih(k_i^2\frac{\mathrm{\Gamma }}{\rho }+ga\mathrm{cos}(\omega t))A_i=0.$$
(4)
An instability occurs and the amplitude $`A_i`$ grows exponentially in time when the eigenvalue is in a band (known as the stability tongue) such that the frequency of oscillation of the fluid is half the driving frequency. The instability occurs at arbitrarily small driving amplitude in the absence of viscosity.

Damping, which is provided by a number of distinct mechanisms in addition to bulk viscosity, can be included by means of a phenomenological linear damping term as reviewed in Ref.. Though treating damping in this way may not be fully adequate, the main effect is to reduce the width of the stability tongue in parameter space and raise the critical threshold to a finite amplitude. A proper theoretical treatment of instability in the viscous case has been given in Ref. , where the shapes of the computed stability boundaries were presented.

If the acceleration $`a`$ is slightly higher than $`a_c`$, all modes in a band ($`𝐤\mathrm{\Delta }𝐤/2,𝐤+\mathrm{\Delta }𝐤/2`$) are accessible and can be excited. The wavenumber width $`\mathrm{\Delta }k`$ of the stability band for small $`ϵ`$ has been estimated  to be:

$$\mathrm{\Delta }k=8\sqrt{2}\rho \nu \omega \sqrt{ϵ}/3\mathrm{\Gamma },$$
(5)
where $`\nu `$ is the kinematic viscosity of the fluid. For a suitable choice of $`ϵ`$, $`\rho `$, and $`\omega `$, and assuming no interaction between modes, one then expects to find either single mode patterns or superpositions of a few modes whose thresholds lie in the window $`(𝐤\mathrm{\Delta }𝐤/2,𝐤+\mathrm{\Delta }𝐤/2)`$.

The cumulative number of eigenvalues of the Helmholtz equation $`N(k)`$ is related to the geometry and is given by:

$$N(k)\frac{S}{4\pi }k^2\frac{L}{4\pi }k,$$
(6)
where $`S`$ is the area, $`L`$ is the perimeter of the stadium, and the negative or positive sign corresponds to Dirichlet or Neumann boundary condition respectively . At high $`k`$, the perimeter term is negligible compared to the area term. Taking $`ϵ=0.01`$ and using Eqs. 5, 6 with an area $`S`$ that is appropriate to the experiments reported here, one can estimate that the typical number of accessible modes is about 8 for a driving frequency of 70 Hz.

What are the effects of nonlinearity? A nonlinear theory that describes regular Faraday wave patterns in large containers rather well has been given by Zhang and Viñals  and Chen and Viñals . In this theory, an evolution equation is determined for the time derivative of the amplitude of a typical Fourier mode $`B_1`$ of the interfacial deformation. It may be expressed in the form

$$\frac{dB_1}{dT}=\alpha B_1g_0B_1^3\underset{m1}{}g(\theta _{m1})B_m^2B_1,$$
(7)
where $`T`$ is a slow time variable, the linear term is due to the basic instability, the cubic self-interaction term produces saturation of the wave pattern, and the coupling terms to other modes (which depend on their relative angle $`\theta `$) are also of cubic order. The constants have been computed, and the ratio $`g(\theta )/g_0`$ is of order unity and independent of $`ϵ`$. This implies that coupling effects between the accessible modes may be substantial. The theory was able to explain the striking cascade of 2n-fold patterns discovered by Binks and van de Water . It also explains semi-quantitatively the appearance of striped, square and hexagonal patterns observed in experiments using viscous fluids in large containers . However, the amplitude equation is variational, and is only appropriate near onset. It cannot describe nonuniform patterns, secondary instabilities, or spatiotemporal chaos. An earlier approach that allowed spatially varying patterns was given by Milner .

The amplitude equations also ignore the effects of the boundary. For slightly viscous fluids in small containers, a large fraction of the dissipation occurs in the boundary layer and can in fact be the leading cause of dissipation . In work stimulated by the experiments reported here, Agam and Altshuler  show that the dissipation near the boundary depends strongly on the nature of the mode.

## III Experimental Apparatus

The apparatus is similar to that used by Gluckman et al. in Ref. . Fig. 1 shows a schematic diagram of the experimental setup. The stadium shaped container made of Delrin has the following dimensions: depth $`d=1.25`$ cm, radius of semicircles $`r=3`$ cm, and length of straight edge $`l=4.5`$ cm. The top and bottom plates of the container are made of glass to allow the transmission of light. The fluid is silicone oil of kinematic viscosity $`0.02\mathrm{cm}^2\mathrm{s}^1`$, chosen for its stable surface tension and good wetting characteristics. To minimize meniscus waves, a brim full boundary condition was prepared by machining a ledge in the boundary at the same height as the fluid. Therefore the fluid meets the ledge at $`90^0`$. By maintaining the fluid under brim full conditions, the fluid surface is pinned to the ledge and boundary dissipation is reduced. This situation has been modeled as a Dirichlet boundary condition ($`\psi =0`$. The container is rigidly attached to an electromagnetic shaker (Vibration Test Systems Model 40C) and the acceleration is measured with an accelerometer. The apparatus is placed within a temperature-controlled environment. The driving frequency is selected to be greater than 55 Hz to be in the capillary wave limit, but less than 75 Hz to prevent the density of modes from becoming too high.

The patterns are imaged with shadowgraph techniques. The specific implementation is discussed in depth in Ref.. Light from an expanded and collimated incident beam is collected and imaged onto a CCD (charge-coupled device) video camera via a large collecting lens and the camera lens. The resulting images can be interpreted by considering which rays of light reach the CCD plane after passing through the fluid. The relatively small aperture of the camera lens restricts the rays that reach the CCD. All the rays that leave the fluid surface at an angle measured from the normal that is greater than a critical angle (typically about $`10^2`$ radians) are blocked. Since the critical angle is so small, light is collected only from the nearly horizontal regions of the wave surface. Therefore, the bright regions in an image corresponds to local extrema or antinodes of the wave pattern. Images are averaged over one video frame, 1/30 s, which is more than a full cycle of the standing waves. The imaging process is nonlinear in the wave height, but a quantitative model for the measured intensity was presented and tested in Ref. .

## IV Patterns Near Onset

We made a survey of the time-independent wave patterns near onset over the range 55 to 65 Hz by changing the frequency in 0.1 Hz steps. In order to obtain useful statistics for the surface wave patterns, a systematic procedure was followed: for each selected frequency, $`a_c`$ was first measured to within 0.1% and then $`ϵ`$ was raised to 0.01, the smallest value that could be maintained accurately. The threshold $`a_c`$ is $`3.1\mathrm{m}\mathrm{s}^2`$ at $`f=60`$ Hz and increases weakly with frequency .

A sequence of images from the survey for driving frequency $`f`$ between 60.1 Hz and 62.8 Hz, and with an approximate spacing of 0.4 Hz (i.e. every fourth image), is shown in Fig. 2. This spacing is comparable to the experimentally observed increment (0.3 Hz) typically required to obtain a distinctly different pattern in this frequency range; it is greater then the computed mean level spacing (about 0.1 Hz in this frequency range as estimated from Eqns. 5 and 6.) Most of the observed patterns show the reflection symmetries of the stadium. Regions of large amplitude are often located along lines that would form periodic orbits of the classical analog. Since these regions are similar to those found in other numerical and experimental investigations , we refer to the patterns containing such enhancements as scarred patterns.

We compare the observed patterns with numerically computed eigenstates of the Laplace operator for the stadium geometry, obtained for comparable mean wavenumber using an algorithm due to Heller . Sample computed eigenstates (selected from a large number of distinct patterns) are shown in Fig. 3. We find that some of the computed states resemble observed patterns. On the other hand, a one-to-one correspondence for sequences of eigenstates was definitely not observed. Furthermore, certain computed eigenstates that occur frequently such as the “whispering gallery” mode (Fig. 3d) were not observed in the full experimental frequency range.

Interestingly, most of the observed symmetric patterns resemble one of three basic classes of eigenstates shown in Fig. 3(a-c). For instance, Figs. 2(a,d) resemble the bouncing ball eigenstate Fig. 3(a); Figs. 2(g,h) are close to the longitudinal eigenstate of Fig. 3(b); and Figs. 2(b,e) are a combination of the longitudinal and bowtie eigenstates of Figs. 3(b,c). It is noteworthy that among the observed patterns are states such as Fig. 2(c) that do not have the reflection symmetries of the stadium. We refer to these as disordered patterns. Table I summarizes the percentages of the onset patterns that were visually judged to approximate particular computed scarred eigenstates in the frequency range 55 to 65 Hz at $`ϵ=0.01`$. (Visual comparison was used because automated pattern recognition, which we attempted, was not sufficiently reliable.)

Since the discovery of scars, there have been a number of theoretical attempts to obtain a quantitative measure for scarring  based on eigenstate overlap, Wigner function overlap, and inverse participation ratios for the amplitudes in the vicinity of the scars. To utilize such measures experimentally, the local wave amplitude is required with high accuracy. The shadowgraph technique used here is quantitative but nonlinear  and does not provide this information. Development of a quantitative experimental measure of scarring has proven to be difficult even for linear probes.

We use the concept of “pattern entropy” as a tool to classify the patterns. Egolf, Melnikov, and Bodenschatz  have applied this concept successfully to measure the complexity of patterns observed in Rayleigh-Bénard convection. The pattern entropy is calculated from the power spectrum of the pattern. If $`P(𝐤)`$ is the normalized two dimensional power spectrum of the pattern at time t, then the pattern entropy $`E(t)`$ is defined as:

$$E(t)=\underset{𝐤}{}P(𝐤)\mathrm{ln}(P(𝐤)).$$
(8)
Here $`E(t)`$ measures the spectral complexity of a pattern. If the image consists of just one Fourier mode of amplitude unity, then $`E=0`$; otherwise $`E>0`$. To minimize the effects of experimental noise, we sum contributions only in a band of wavenumbers centered at the mean wavenumber of the pattern with a rage of $`\pm 25\%`$. In Table II, the approximate entropy ranges for the various types of patterns observed in the range 55-65 Hz are given. Note that the patterns are not distinguishable solely by their entropy, since some of the ranges overlap. However, the pattern entropy can be useful in studies of time dependence farther above onset, as we show in Sec. VI.

## V Patterns beyond onset

Here we examine the evolution of the wave patterns farther from onset, where the interactions between different Fourier components of the waves become increasingly nonlinear and the approximation of Eq. (2) becomes inapplicable. The patterns were observed to be time-independent while changing adiabatically with $`ϵ`$ for $`ϵ<0.3`$. On the other hand, they become weakly time-dependent for $`ϵ0.3`$ at most frequencies.

The evolution toward time-dependence with increasing $`ϵ`$ depends on the excitation frequency. Three examples of this evolution are shown in Figs. 4-6. For some driving frequencies, the patterns remain reflection-symmetric as $`ϵ`$ is increased, but exhibit transitions from one spatial mode to another prior to the onset of time dependence, as in Fig. 5. In these cases, the transition to spatial disorder (asymmetry) tends to coincide with the onset of spatiotemporal chaos (STC).

It is important to note that as $`ϵ`$ increases, the width of the stability tongue grows (see Eq.(5)): for example at $`f=74.1`$ Hz and $`ϵ=0.2`$ the number of accessible modes of the container is approximately 35. Therefore, the observed mode switching might be a combined effect of the increase in the number of accessible modes and an increase in the degree of nonlinearity that couples them. It is remarkable that the container boundary continues to influence the patterns even at $`ϵ=0.252`$ (Fig. 5d), where nonlinearity clearly plays a major role.

The variability of the nonlinear development is evident from examining the examples in Figs. 4-6. In Fig. 4, there is a general increase in complexity with $`ϵ`$, but the dominant mode does not change. In Fig. 6, the near onset pattern is nearly obliterated even at $`ϵ=0.025`$ by the growing complexity, and the pattern is also distinctly asymmetric, while remaining time-independent.

In one instance, a complex but symmetric time-independent pattern was observed at an unusually high driving amplitude of $`ϵ=0.8`$ at a frequency of 65.0 Hz, in a regime where spatiotemporal chaos is usually fully developed. The image shown in Fig. 7 was averaged over 3000 images taken over a period of 5 minutes to test for time-dependence. The lack of blurring demonstrates its time-independence.

## VI Role of ordered states in the regime of spatiotemporal chaos

For driving amplitudes just beyond the frequency-dependent onset of spatiotemporal chaos, the time dependence of the pattern is often intermittent; the patterns appear to oscillate between states that are relatively ordered and states that are relatively disordered (see the images in Fig. 8 and a corresponding web-based movie ). The power spectra of these typical patterns are also shown in Fig. 8 and indicate the greater complexity of the disordered case, where the power is distributed more uniformly on the ring corresponding to the dominant wavenumber. The time dependence and complexity of the patterns are monitored using two quantitative measures: (i) the rate of change of the pattern $`R(t)`$, and (ii) the entropy $`E(t)`$ as defined in Eq. (8). The rate of change $`R(t)`$ was calculated by subtracting two consecutive images I($`𝐱,t+\mathrm{\Delta }t`$) and I($`𝐱,t`$) separated by a time interval of $`\mathrm{\Delta }t=0.36`$ seconds, and calculating $`R(t)`$ according to the following formula:

$$R(t)=c(I(𝐱,t+\mathrm{\Delta }t)I(𝐱,t))^2,$$
(9)
where, $`𝐱`$ is the position, and $`c`$ is a constant scaling factor.

In Fig. 9(a), a section of the resulting quantity $`R(t)`$ is shown as function of time for $`f=71.9`$ Hz and $`ϵ=0.55`$. The data has been smoothed by averaging over four adjacent points. At every pronounced minimum of this smoothed data we find that the corresponding pattern is symmetric and appears to have long range order. At all other times the pattern is asymmetric and disordered. A graph of the corresponding entropy $`E(t)`$ is shown in Fig. 9(b). Examples of the ordered (X) and disordered (Y) states are shown in Fig. 8, where the location in time is indicated by symbols X ($`t=124`$ s) and Y ($`t=136`$ s) respectively in Fig. 9. At X, where the pattern entropy is low, $`R(t)`$ is small, while at Y, where the pattern entropy is high, $`R(t)`$ is large. This correlation between $`R(t)`$, and $`E(t)`$ holds true for most of the other strong peaks and valleys.

For higher $`ϵ`$ ($`>0.8`$) the oscillations diminish in strength and uniform STC is observed. In this regime, following Gluckman et al. , we obtained the time averaged pattern after adding 3000 instantaneous images obtained over five minutes. An example of an instantaneous pattern is shown in Fig. 10(a). The resultant average pattern is shown in Fig. 10(b). The average reveals considerable ordered structure including remnants of the bouncing ball states. This phenomenon is seen at most frequencies and persists up to $`ϵ1.6`$. Beyond this point, scars were visually absent and the average patterns are locally parallel to the boundary, as observed in circular or square patterns by Gluckman et al. .

## VII Discussion

In this paper we have discussed the parametrically forced wave patterns formed in a stadium-shaped container containing a low viscosity fluid, as a function of driving frequency and amplitude. The patterns near onset ($`ϵ=0.01`$) were compared to a simple model consisting of linearized equations that reduce to the Helmholtz equation with Dirichlet boundary conditions (see Section II). While a large proportion of observed patterns resemble the numerically computed eigenstates of the stadium, many of the computed eigenstates (for instance the whispering gallery modes) were not observed in a scan with sufficient frequency resolution to detect them. The observed patterns may be broadly classified into three categories: (a) bouncing ball patterns, (b) longitudinal patterns, and (c) bowtie patterns, which have high amplitudes near corresponding periodic orbits. In addition, a significant number of disordered patterns (lacking in symmetry but time-independent) were observed near onset. Furthermore, the observed mode spacing ( $`0.3`$ Hz) is somewhat greater than the mean eigenvalue separation implied by Eq. (6) ($`0.1`$ Hz). These observations imply that the simplest model is inadequate even close to onset.

Recently, Agam and Altshuler have offered an explanation for the selection of modes at onset  by considering the stability of the periodic orbits corresponding to the scars. They argue that a threshold for excitation of a particular scarred pattern is given by:

$$h>\gamma _b+\gamma _p+\lambda /2,$$
(10)
where $`h`$ is proportional to the rate that energy is pumped into the system (i.e. the driving amplitude), $`\gamma _b`$ is the dissipation rate in the bulk of the fluid, $`\gamma _p`$ is the dissipation near the perimeter, and $`\lambda `$ is the Lyapunov exponent of the ray orbit that predominantly scars the pattern. The bulk dissipation $`\gamma _b=\nu k^2`$ is the same for all the patterns and corresponds to approximately 2 sec<sup>-1</sup> in the frequency window used in the experiments. Therefore the appearance of a scarred pattern depends on a combination of the two remaining factors which are orbit-dependent. Stability of a pattern is favored both by a small Lyapunov exponent of the associated scarring orbit, and by small perimeter dissipation $`\gamma _p`$.

In the limit of high wavenumber $`k`$ relevant to our experiments, Agam and Altshuler derive an expression for the damping rate of scars due to boundary effects:

$$\gamma _p=\frac{\sqrt{\omega \nu /2}}{L}\underset{i}{}\frac{(1\mathrm{cos}^2(\varphi _i))}{\mathrm{cos}(\varphi _i)}$$
(11)
where $`\omega `$ is the angular frequency, $`\nu `$ is the viscosity, $`L`$ is the length of the periodic orbit, and the sum is over all the collision points of the orbit with the boundary, $`\varphi _i`$ being the angle between the orbit and a line perpendicular to the boundary at the collision point. The parameter $`\gamma _p`$ and the Lyapunov exponent have been calculated for most of the shortest periodic orbits and a long ergodic orbit in Ref. . The perimeter dissipation $`\gamma _p`$, which can be either smaller or larger than $`\gamma _b`$, varies between 0.2 sec<sup>-1</sup> and 3.0 sec<sup>-1</sup>. The extreme values correspond to patterns scarred by the horizontal orbit and the ergodic orbit respectively.

Most of the features observed in the experiments near onset appear to be captured by Eq. (10). The bouncing ball orbit is prominent because the Lyapunov exponent is zero in that case, and the longitudinal orbit occurs because of relatively low perimeter dissipation. The whispering gallery orbits and others with angles that come close to $`\pi /2`$ have particularly large boundary dissipation and are suppressed. Eq. 10 implies that if one increases the dissipation on the boundary so that it dominates, the longitudinal orbit will be the last to survive. This is precisely what we observe in the experiments when the level of the fluid is lowered, a change that results in higher perimeter dissipation because of motion of the contact line.

The theory just cited  is also able to account for the observed tendency of one scarred pattern to suppress other nearby eigenmodes through nonlinear interactions, as well as the existence of some asymmetric patterns. At some driving frequencies the patterns are observed to switch modes as $`ϵ`$ is increased (Fig. 5). This occurs especially when the bowtie state is observed at onset, an observation that may be related to the larger perimeter dissipation and Lyapunov exponent of the bowtie mode.

At higher driving amplitude, additional nonlinear effects occur as indicated by the growth in spatial complexity, and no adequate theoretical treatment exists. However, the onset patterns are often robust, persisting in the presence of increasing spatial complexity. The boundaries remain influential even beyond the onset of time dependence. At sufficiently high $`ϵ`$ ($`0.2`$), the onset patterns are no longer visible, though they persist in the time-average.

Strong intermittency in the degree of order of the patterns is observed in the regime of spatiotemporal chaos. Furthermore, the rate of change of the pattern $`R(t)`$ just above the STC onset is strongly correlated with the order as characterized by the entropy $`E(t)`$. The more ordered patterns evolve more slowly in time, a striking observation that remains to be explained. This tendency for ordered patterns to be more stable may be related to the complex but highly symmetric and time-independent pattern observed at atypically large acceleration (Fig. 7).

## VIII Acknowledgments

We thank Oded Agam and Boris Altshuler for many useful discussions, and E.J. Heller for software to determine eigenstates of the stadium. Bruce Boyes provided technical help. This work was supported by the National Science Foundation under Grant DMR-9704301. A. K. acknowledges a grant from the Alfred P. Sloan Foundation.


text:
# Small-eccentricity binary pulsars and relativistic gravity

## 1. Introduction

The majority of binary pulsars is found to be in orbit with a white-dwarf companion. Due to the mass transfer in the past, these systems have very small orbital eccentricities and, therefore, neither the relativistic advance of periastron nor the Einstein delay were measured for any of these binary pulsars. In fact, the only post-Keplerian parameters measured with reasonable accuracy for a small-eccentricity binary pulsar are the two Shapiro parameters in case of PSR B1855+09 (Kaspi et al. 1994). However, since the orbital period of this system is 12.3 days, the expected gravitational wave damping of the orbital motion is by far too small to be of any importance for timing observations and, consequently, there is no third post-Keplerian parameter which would allow the kind of test conducted in double-neutron-star systems (Damour & Taylor 1992). On the other hand, many alternative theories of gravity, tensor-scalar theories for instance, predict effects that depend strongly on the difference between the gravitational self energy per unit mass ($`ϵE^{\mathrm{grav}}/mc^2`$) of the two masses of a binary systems (Will 1993; Damour & Esposito-Farése 1996ab). While this difference in binding energies is comparably small for double-neutron-star systems, it is large in neutron star-white dwarf systems since for a white dwarf $`ϵ10^4`$ while for a 1.4$`M_{}`$ neutron star $`ϵ0.15`$.

## 2. Gravitational dipole radiation

Unlike general relativity, many alternative theories of gravity predict the presence of all radiative multipoles — monopole and dipole, as well as quadrupole and higher multipoles (Will 1993). For binary systems scalar-tensor theories, for instance, predict a loss of orbital energy which at highest order is dominated by scalar dipole radiation. As a result, the orbital period, $`P_b`$, of a circular binary system should change according to

$$\dot{P}_b^{(\mathrm{dipole})}\frac{4\pi ^2G_{}}{c^3P_b}\frac{m_pm_c}{m_p+m_c}(\alpha _p\alpha _c)^2,$$
(1)
where $`m_p`$ and $`m_c`$ denote the mass of the pulsar and its companion, respectively, $`G_{}`$ is the ‘bare’ gravitational constant and $`c`$ the speed of light. The total scalar charge of each star is proportional to its mass and its ‘effective coupling strength’ $`\alpha (ϵ)`$ (Damour & Esposito-Farése 1996b). For a white dwarf companion $`|\alpha _c|1`$ and thus the expression $`(\alpha _p\alpha _c)^2`$ in equation (1) can be of the order one if the pulsar develops a significant amount of scalar charge. In this case the gravitational wave damping of the orbit is completely dominated by the emission of gravitational dipole radiation.

PSR J1012+5307 is a 5.3 ms pulsar in a 14.5 h circular orbit with a low mass white-dwarf companion. Since its discovery in 1993 (Nicastro et al. 1995) this pulsar has been timed on a regular basis using the Jodrell Bank 76-m and the Effelsberg 100-m radiotelescope, sometimes achieving a timing accuracy of 500 ns after just 10 min of integration (Lange et al. , this conference). In addition, the white-dwarf companion appears to be relatively bright ($`V=19.6`$) and shows strong Balmer absorption lines. Based on white dwarf model calculations, a companion mass of $`m_c=0.16\pm 0.02`$ and a distance of $`840\pm 90`$ pc was derived (van Kerkwijk et al. 1996, Callanan et al. 1998). Further, a reliable radial velocity curve for the white dwarf has been extracted, which then, in combination with the pulsar timing information, gave a mass for the pulsar of $`m_p=1.64\pm 0.22`$. Since $`\dot{P}_b=(0.1\pm 1.8)\times 10^{13}`$ for this binary system, we find from equation (1)

$$|\alpha _p|<0.02\text{(95\% C.L.)}$$
(2)
Simulations show, that this value should improve by a factor of five within the next three years (Lange et al., in prep.).

## 3. Gravitational Stark effects

### 3.1. Violation of the strong equivalence principle

The strong equivalence principle (SEP) requires the universality of free fall of all objects in an external gravitational field regardless of their mass, composition and fraction of gravitational self-energy. While all metric theories of gravity share the property of universality of free fall of test particles (weak equivalence principle), many of them, which are considered as realistic alternatives to general relativity, predict a violation of the SEP. A violation of the SEP can be understood as an inequality between the gravitational mass, $`m_g`$, and the inertial mass, $`m_i`$, which can be written as function of $`ϵ`$:

$$m_g/m_i1+\delta (ϵ)=1+\eta ϵ+𝒪(ϵ^2).$$
(3)
While the analysis of lunar-laser-ranging data tightly constrains the ‘Nordtvedt parameter’ $`\eta `$ (Müller et al. 1997) it indicates nothing about a violation of the SEP in strong-field regimes, i.e. terms of higher order in $`ϵ`$, due to the smallness of $`ϵ`$ for solar-system bodies. For neutron stars, however, $`ϵ0.15`$ and thus binary-pulsars with white-dwarf companions ($`ϵ10^4`$) provide ideal laboratories for testing a violation of the SEP due to nonlinear properties of the gravitational interaction (Damour & Schäfer 1991).

In case of a violation of the SEP the eccentricity vector of a small-eccentricity binary-pulsar system exposed to the external gravitational field of the Galaxy, $`𝐠`$, is a superposition of a constant vector $`𝐞_F`$ and a vector $`𝐞_R`$ which is turning in the orbital plane with the rate of the relativistic advance of periastron. The ‘induced’ eccentricity $`𝐞_F`$ points into the direction of the projection of the Galactic acceleration onto the orbital plane, $`𝐠_{}`$, and $`e_F(\delta _p\delta _c)P_b^2g_{}`$. However, neither the length of $`𝐞_R`$ nor its rotational phase $`\theta `$ are known quantities. We therefore have to proceed as follows. Given a certain $`(\delta _p\delta _c)\delta _p`$, i.e. a certain $`e_F`$ for a given binary pulsar, the observed eccentricity $`e`$ sets an upper limit to $`|\theta |`$ which is independent of $`e_R`$: $`\mathrm{sin}|\theta |<e/e_F`$ for $`e<e_F`$ and $`|\theta |\pi `$ for $`ee_F`$ (Wex 1997). We now have to calculate an upper limit for $`\theta `$ for every observed small-eccentricity binary pulsar and compare the result with Monte-Carlo simulations of a large number of (cumulative) distributions for the (uniformly distributed) angle $`\theta `$. This way, by counting the number of simulated distributions which are in agreement with the ‘observed’ limits, one obtains the confidence level with which a certain $`\delta _p`$ is excluded. As a safe upper limit for $`|\delta _p|`$ we find

$$|\delta _p|<0.009\text{(95\% C.L.)}$$
(4)
Note, in order to calculate $`e_F`$ for a given binary system, we need also the masses of pulsar and companion and the location and orientation of the binary system in the Galaxy. If there are no restrictions from timing and optical observations, the pulsar masses were assumed to be uniformly distributed in the range $`1.2M_{}<m_p<2M_{}`$, the companion masses were taken from evolutionary scenarios (Tauris & Savonije 1999), and the pulsar distances were estimated using the Taylor-Cordes model assuming a typical error of 25% (Taylor & Cordes 1993). Finally, the orientation of the ascending node in the sky, which is an unobservable parameter for all our binary system, was treated as variable which is uniformly distributed between 0 and $`2\pi `$.

### 3.2. Violation of local Lorentz invariance and conservation laws

If gravity is mediated in part by a long-range vector field or by a second tensor field one expects the global matter distribution in the Universe to select a preferred frame for the gravitational interaction (Will & Nordtvedt 1972). At the post-Newtonian level, gravitational effects associated with such a violation of the local Lorentz invariance of gravity are characterized by two theory dependent parameters $`\alpha _1`$ and $`\alpha _2`$. If $`\alpha _1`$ were different from zero, the eccentricity of a binary system which moves with respect to the global matter distribution in the Universe would suffer a secular change similar to a violation of the SEP. This time, $`|e_F|\alpha _1|m_pm_c|P_b^{1/3}w_{}`$ where $`𝐰`$ denotes the velocity of the binary system with respect to the preferred frame, i.e. the cosmic microwave background. Again, we can perform a Mote-Carlo analysis as outlined in the previous section to derive the upper limit

$$|\alpha _1|<1.2\times 10^4\text{(95\% C.L.)}$$
(5)
This limit is slightly better than the limit obtained from lunar-laser-ranging data (Müller et al. 1996) and, more importantly, also holds for strong gravitational-field effects which could occur in the strong-field regions of neutron stars. Due to its small eccentricity, $`e<1.7\times 10^6`$ (95% C.L.), and high velocity with respect to the cosmic microwave background ($`w500`$ km/s), PSR J1012+5307 turns out to be the most important binary system for this kind of analysis (Lange et al. , in prep). While for PSR J1012+5307 also the radial velocity of the system is known from spectroscopic observations of the white dwarf companion, for all the other binary pulsars no radial velocity information is available and we have to assume an isotropic probability distribution for the 3-d velocity.

In theories of gravity which violate the local Lorentz invariance and the momentum conservation law, a rotating self-gravitating body will suffer a self-acceleration which is given by $`𝐚_{\mathrm{self}}=\frac{1}{3}\alpha _3ϵ𝐰\times 𝛀`$ (Nordtvedt & Will 1972), where $`\alpha _3`$ is a theory dependent parameter and $`𝛀`$ denotes the rotational velocity of the body. Again, binary pulsars are ideal probes for this kind of self-acceleration effects (Bell & Damour 1996). A careful analysis analogous to the previous analyses (SEP, local Lorentz invariance) gives

$$|\alpha _3|<1.5\times 10^{19}\text{(95\% C.L.)}$$
(6)
Note, the statistical tests for gravitational Stark effects presented here for the first time appropriately take care of selection effects by simulating the whole population, therefore, giving the first reliable limits for $`\delta _p`$, $`\alpha _1`$, and $`\alpha _3`$.

#### Acknowledgments.

I am grateful to Kenneth Nordtvedt for pointing out to me the problem of selection effects related with binary-pulsar limits to gravitational Stark effects. I thank Christoph Lange for numerous valuable discussions.

## References

J. F. Bell & T. Damour: Class. Quantum Grav., 13, 3121 (1996)

P. J. Callanan, P. M. Garnavich, D. Koester: MNRAS, 298, 207 (1998)

T. Damour & G. Esposito-Farése: Phys. Rev. D, 53, 5541 (1996a)

T. Damour & G. Esposito-Farése: Phys. Rev. D, 54, 1474 (1996b)

T. Damour & G. Schäfer: Phys. Rev. Lett., 66, 2550 (1991)

T. Damour & J. H. Taylor: Phys. Rev. D, 45, 1840 (1992)

V. M. Kaspi, J. H. Taylor & M. F. Ryba: ApJ, 428, 713 (1994)

J. Müller, K. Nordtvedt & D. Vokrouhlický: Phys. Rev. D, 54, R5927 (1996)

J. Müller, M. Schneider, K. Nordtvedt & D. Vokrouhlicky: In: Proceedings of the 8th Marcel Grossman Meeting, Jerusalem 1997

L. Nicastro, A. G. Lyne, D. R. Lorimer, P. A. Harrison, M. Bailes, B. D. Skidmore: MNRAS, 273, L68 (1995)

K. Nordtvedt & C. M. Will: ApJ, 177, 775 (1972)

T. M. Tauris & G. J. Savonije: A&A, 350, 928 (1999)

J. H. Taylor & J. M. Cordes: ApJ, 411, 674 (1993)

M. H. van Kerkwijk, P. Bergeron & S. R. Kulkarni: ApJ, 467, L89 (1996)

N. Wex: A&A, 317, 976 (1997)

C. M. Will: Theory & experiment in gravitational physics, (Cambridge University Press, Cambridge 1993)

C. M. Will & K. Nordtvedt: ApJ, 177, 757 (1972)