text:
# On disc driven inward migration of resonantly coupled planets with application to the system around GJ876

## 1 Introduction

The recent discovery of extrasolar giant planets orbiting around nearby solar–type stars (Marcy & Butler 1998, 2000) has stimulated renewed interest in the theory of planet formation. The objects observed so far have masses, $`M_p`$, that are characteristic of giant planets

($`0.4\mathrm{M}_J\stackrel{<}{}M_p\stackrel{<}{}\mathrm{\hspace{0.33em}11}\mathrm{M}_J`$), $`M_J`$ denoting a Jupiter mass. The orbital semi-major axes are in the range

$`0.04\mathrm{AU}\stackrel{<}{}a\stackrel{<}{}\mathrm{\hspace{0.33em}2.5}\mathrm{AU}`$, and orbital eccentricities in the range $`0.0\stackrel{<}{}e\stackrel{<}{}\mathrm{\hspace{0.33em}0.67}`$ (Marcy & Butler 2000).

Disc-protoplanet interactions have been invoked to explain the presence of giant planets orbiting close to their host stars through inward orbital migration induced through disc-protoplanet tidal interaction ( eg. Papaloizou, Terquem & Nelson 1999, Lin et al 2000). Up to now, for the most part extrasolar planets appear to be isolated. However, a few multiple systems are known. The configuration of these may contain important information about their origin and possible migration history. Of special interest is the recently discovered system around GJ876. This is found to be close to a $`2:1`$ commensurability. Such a configuration is indeed suggestive of orbital migration. Commensurable satellite systems such as the Galilean satellites are thought to owe their origin to migration induced by tidal interaction with the central planet ( eg. Goldreich 1965).

Recent simulations of single protoplanets in the observed mass range (Kley 1999, Bryden et al. 1999, Lubow, Seibert & Artymowicz 1999) interacting with a disc with parameters thought to be typical of protoplanetary discs, but constrained to be in circular orbit, indicate gap formation and upper mass limit consistent with the observations. Nelson et al. (2000) relaxed the assumption of fixed circular orbits, found inward migration and that the disc-protoplanet interaction leads to strong eccentricity damping. Due to accretion onto the central star, an inner cavity was formed in the disc interior to which the protoplanet orbited.

Simulations of two planets interacting with a disc have been performed by Kley (2000), Bryden et al (2000), and Masset & Snellgrove (2001). So far inward migration of two planets locked into a $`2:1`$ commensurability has not been simulated. However, Bryden et al (2000) found a tendency for gap material between the two planets in fixed circular orbits to be cleared, ending up interior to the inner planet orbit or exterior to the outer planet orbit.

Taken together the above results suggest a natural outcome of two protoplanets interacting with a disc is that they orbit interior to an inner disc cavity while the external disc causes inward migration of the outer orbit. This catches up the inner orbit leading to the possibilty of resonant interaction.

It is the purpose of this paper to investigate such resonant interaction and whether, for reasonable protoplanetary disc models, it leads to a locking such that the planets subsequently migrate maintaining the commensurability. Similar behaviour occurs as a result of the tidally induced migration of the Gallilean satellites (eg. Goldreich 1965, Lin & Papaloizou 1979).

Here we shall assume the prior evolution of the system leads to the orbital separation of the planets being slightly larger than that required for a strict $`2:1`$ commensurability without considering the history in detail as it is beyond the scope of this paper. However, we comment that this might have been complicated with the planet masses varying with time through mass accretion from the disc.

Analytic methods, N body integrations and direct two dimensional numerical simulations are used to investigate the evolution and found to give consistent results.

We consider two protoplanets gravitationally interacting with each other and a protoplanetary disc. The two planets orbit interior to a tidally maintained disc cavity while the disc interaction induces inward migration. This migration is of type II and so is regulated by the magnitude of the disc viscosity. When the migration is slow enough, the more rapidly migrating outer protoplanet approaches and becomes locked in a $`2:1`$ mean motion commensurability with the inner one. The commensurability persists in subsequent evolution. The eccentricities of the protoplanets are increased by the resonant perturbations and damped by circularization which occurs through interaction with the disc (eg. Goldrech & Tremaine 1981). A balance is achieved in which the protoplanet eccentricities are determined by the migration rate and circularization rate induced in the outer planet orbit by the external disc.

We apply our results to the recently discovered resonant planets around GJ876. Simulation shows that migration induced by a disc with parameters expected for protoplanetary discs results in trapping in the $`2:1`$ commensurability when the planets orbit in an inner cavity. Eccentricities in the observed range may be obtained. Further studies using $`N`$ body integrations indicate that the planetary system will remain stable for at least $`2\times 10^7`$ orbits when the external disc is removed.

In section (2) we describe a simple analytic model of two migrating protoplanets in a $`2:1`$ commensurability. The eccentricities of the protoplanet orbits are related to the migration rate, and circularization rate induced by the disc. In section (3) we describe a simulation of two planets orbiting in an inner disc cavity. Parameters appropriate to GJ867 are adopted. This demonstrates resonant trapping and that eccentricities of the observed magnitude may be produced.

In section (4) we describe N body calculations confirming the above conclusions and indicating the long term stability of the system. Finally in section (5) we discuss our results.

## 2  A simple model

We consider a system consisting of $`2`$ planets and a primary star moving under their gravitational attraction. When there are no disc interactions and the motion is conservative the system is conveniently expressed in Hamiltonian form using Jacobi coordinates (eg. Sinclair 1975). The coordinates, $`𝐫_2,`$ of inner planet of reduced mass $`m_2`$ are referred to the central star of mass $`M_{}`$ and the coordinates of the outer planet, $`𝐫_1,`$ of reduced mass $`m_1`$ are referred to the centre of mass of the central star and inner planet. The Hamiltonian can be written correct to second order in the planetary masses as

$`H`$ $`=`$ $`{\displaystyle \frac{1}{2}}(m_1|\dot{𝐫}_1|^2+m_2|\dot{𝐫}_2|^2){\displaystyle \frac{GM_1m_1}{|𝐫_1|}}{\displaystyle \frac{GM_2m_2}{|𝐫_2|}}`$ (1)
$``$ $`{\displaystyle \frac{Gm_1m_2}{|𝐫_{12}|}}+{\displaystyle \frac{Gm_1m_2𝐫_1𝐫_2}{|𝐫_1|^3}}.`$
Here $`M_1=M_{}+m_1,M_2=M_{}+m_2`$ and $`𝐫_{12}=𝐫_2𝐫_1.`$ The Hamiltonian can be expressed in terms of the osculating semi-major axes, eccentricities and longtitudes of periastron $`a_i,e_i,\varpi _i,i=1,2`$ respectively as well as the longtitudes $`\lambda _i,`$ and the time $`t.`$ We recall that $`\lambda _i=n_i(tt_{0i})+\varpi _i,`$ with $`n_i`$ being the mean motion and $`t_{0i}`$ giving the time of periastron passage. The energy is given by $`E_i=Gm_iM_i/(2a_i),`$ and the angular momentum $`h_i=m_i\sqrt{GM_ia_i(1e_i^2)}`$ which may be used to describe the motion instead of $`a_i`$ and $`e_i.`$

Only terms first order in the eccentricities and involving the resonant angles $`\varphi =2\lambda _1\lambda _2\varpi _1,`$ and $`\psi =2\lambda _1\lambda _2\varpi _2,`$ are retained in the expansion (see Sinclair 1975). Then the perturbing part of the Hamiltonian $`(m_1m_2)`$ can be written

$$H_{12}=\frac{Gm_1m_2}{a_1}\left(Be_1\mathrm{cos}\varphi Ce_2\mathrm{cos}\psi \right),$$
(2)
where $`C=2b_{1/2}^{(2)}(\alpha )+(1/2)\alpha db_{1/2}^{(2)}(\alpha )/d\alpha `$ and $`B=1.5b_{1/2}^{(1)}(\alpha )+(1/2)\alpha db_{1/2}^{(1)}(\alpha )/d\alpha 2\alpha .`$ Here $`b_{1/2}^{(n)}`$ denotes the usual Laplace coefficient and $`\alpha =a_2/a_1.`$ From now on we replace $`M_i`$ by $`M_{}.`$

### 2.1 Orbital precession

We may also take into account additional gravitational forces that may produce precession of the planetary orbits. These could result from the mass distribution provided by the disc exterior to the outer planet or the time averaged mass distribution of the planets themselves. These effects are not included in the resonant Hamiltonian given above. Accordingly we add to it another Hamiltonian

$$H_{int}=\frac{1}{2}m_1n_1a_1^2\omega _{pr1}e_1^2\frac{1}{2}m_2n_2a_2^2\omega _{pr2}e_2^2.$$
(3)
Here, as can be verified from the equations of motion, we have adopted a parameterization such that the precession frequency induced in the orbit of $`m_i`$ is $`\omega _{pri}(a_i).`$ This prescription does not allow the precession frequecy to depend on other quantities. However, appropriate matching can be carried out for a particular case under consideration.

### 2.2 Basic Equations

The equations of motion are derived from:
$`dE_i/dt=n_iH^{}/\lambda _i(n_1T+D)\delta _{i1},`$
$`dh_i/dt=H^{}/\lambda _iH^{}/\varpi _iT\delta _{i1},`$
$`d\lambda _i/dt=n_i+n_iH^{}/E_i+H^{}/h_i,`$
$`d\varpi _i/dt=H^{}/h_i,`$
with $`H^{}=H_{12}+H_{int}.`$ These can be obtained from Hamilton’s equations (eg. Brouwer & Clemence 1961) to which we have added, for the outer planet $`m_1,`$ an additional external torque $`T`$ with an associated orbital energy loss rate $`n_1T`$ together with additional orbital energy dissipation rate $`D.`$ The torque and dissipation rate could be produced by tidal interaction with the disc leading to inward migration and orbital circularization.

We thus obtain to lowest order in the planetary eccentricities and perturbing masses.

$`{\displaystyle \frac{dn_1}{dt}}`$ $`=`$ $`{\displaystyle \frac{6n_1^2m_2}{M_{}}}\left(Be_1\mathrm{sin}\varphi Ce_2\mathrm{sin}\psi \right)`$ (4)
$`+`$ $`{\displaystyle \frac{3n_1a_1}{GM_{}m_1}}\left(n_1T+D\right)`$
$$\frac{dn_2}{dt}=\frac{3n_2^2m_1a_2}{M_{}a_1}\left(Be_1\mathrm{sin}\varphi Ce_2\mathrm{sin}\psi \right)$$
(5)
$$\frac{de_1}{dt}=\frac{m_2n_1}{M_{}}B\mathrm{sin}\varphi \frac{Da_1}{GM_{}m_1e_1}$$
(6)
$$\frac{de_2}{dt}=\frac{m_1n_2a_2}{M_{}a_1}C\mathrm{sin}\psi $$
(7)
$$\frac{d\varphi }{dt}=2n_1n_2\frac{m_2}{e_1M_{}}n_1B\mathrm{cos}\varphi \omega _{pr1}$$
(8)
$$\frac{d\psi }{dt}=2n_1n_2+\frac{m_1a_2}{e_2a_1M_{}}n_2C\mathrm{cos}\psi \omega _{pr2}$$
(9)
### 2.3 Stationary solutions

When no migration or circularization occurs ($`T=D=0`$) equilibrium solutions may exist such that $`\psi `$ and $`\varphi `$ are either zero or $`\pi .`$ Each of $`n_1,n_2,e_1,e_2`$ are then constant.

A relation between the eccentricities then follows from (8) and (9) in the form

$`e_2a_1(e_1M_{}\omega _{pr1}+m_2n_1B\mathrm{cos}\varphi )=`$
$`e_1(m_1a_2n_2C\mathrm{cos}\psi +e_2a_1M_{}\omega _{pr2}).`$ (10)
This condition matches the precession rates of the orbits of the two planets. Also $`2n_1=n_2.`$ Noting that the eccentricities are positive, when they are of very small magnitude, the precessional terms become negligible and there is a solution with $`\psi =0,\varphi =\pi `$ or $`\psi =0,\varphi =\pi .`$ In either case we have

$$m_2a_1e_2n_1B=m_1a_2e_1n_2C.$$
(11)
For larger eccentricities the precessional terms may become important in (10) and then solutions with $`\psi =0,\varphi =0,`$ may occur. Then (10) gives

$$m_2a_1e_2n_1B+m_1a_2e_1n_2C=e_1e_2a_1M_{}(\omega _{pr2}\omega _{pr1}).$$
(12)
For stable solutions, when perturbed, the angles may undergo librations about their equilibrium points (eg. Sinclair 1975). There are two frequencies of oscillation $`\nu _1,\nu _2`$ being given for any $`\psi ,\varphi `$ in the limit of small eccentricities by $`\nu _1^2=(m_2n_1B)^2/(M_{}e_1)^2,`$ and $`\nu _2^2=(m_1n_2a_2C)^2/(a_1M_{}e_2)^2.`$ We look for a solutions with migration which are close to stable solutions of this type.

### 2.4 Resonant Migration

We look for solutions of (4-9) corresponding to the situation where the two planets migrate inwards locked in resonance with $`n_1/n_2`$ maintained nearly equal to $`1/2`$ while the eccentricities remain nearly constant. The tendency for the resonant coupling to excite the eccentricities is counterbalanced by circularization through the action of $`D(GM_{}m_1e_1^2)/(a_1t_c)`$ which defines the circularization time for $`e_1.`$ Similarly (4) defines an inward migration timescale $`t_{mig}=GM_{}m_1/(3Ta_1n_1).`$

We begin by supposing that the angle $`\psi `$ executes a libration about zero such that the mean rate of change of $`e_2`$ is zero. Similarly the mean rates of change of $`n_1`$ and $`n_2`$ induced by $`\psi `$ are zero. Such a libration is seen in simulations. We also suppose the angle $`\varphi `$ either librates or circulates but in such a way that the correspondingly induced mean rates of change are not zero. The simplest example is when the angle executes a very small or even zero amplitude libration about a value slightly offset from zero or $`\pi `$ (eg. Lin & Papaloizou 1979).

We suppose the circulation/libration periods to be short compared to the timescale of migration so that averaging is possible. We denote the average of $`n_1e_1\mathrm{sin}\phi `$ by $`\delta .`$ Then averaging (4-9) gives

$$\overline{\frac{1}{n_1}\frac{dn_1}{dt}}=\frac{6m_2B\delta }{M_{}}+\frac{3a_1}{GM_{}m_1}\left(n_1T+D\right)$$
(13)
$$\overline{\frac{1}{n_2}\frac{dn_2}{dt}}=\frac{6m_1Ba_2\delta }{M_{}a_1}$$
(14)
$$\overline{\frac{de_1^2}{dt}}=\frac{2m_2B\delta }{M_{}}\frac{2Da_1}{GM_{}m_1}0$$
(15)
$$\overline{\frac{d\varphi }{dt}}=\frac{m_2}{e_1M_{}}n_1B\mathrm{cos}\varphi \frac{m_1a_2}{e_2M_{}a_1}n_2C\omega _{pr1}+\omega _{pr2}.$$
(16)
Here on the left hand sides the overline denotes the time average and for simplicity of notation this has been dropped from the right hand sides. The resonance condition also implies $`\overline{\frac{1}{n_2}\frac{dn_2}{dt}}=\overline{\frac{1}{n_1}\frac{dn_1}{dt}}.`$ Using this and eliminating $`\delta `$ from the above, we obtain for the rate of increase of $`n_1`$ through migration

$$\overline{\frac{1}{n_1}\frac{dn_1}{dt}}=\frac{3a_1}{GM_{}m_1}\left(n_1T+D\right)\left(\frac{m_1a_2}{m_1a_2+m_2a_1}\right).$$
(17)
Also (15) gives for the eccentricity balance

$$e_1^2=\frac{t_cm_2a_1}{3t_{mig}(2m_1a_2+m_2a_1)}.$$
(18)
The above determines the eccentricity of the outer planet $`e_1`$ as a function of $`t_c`$ and $`t_{mig}.`$ For a system with $`m_1/m_2=3,`$ we get $`e_1\sqrt{0.07t_c/t_{mig}}.`$

For $`e_1`$ in the $`0.01`$ range we need $`t_c10`$ orbits if $`t_{mig}10^4`$ orbits.

The eccentricity of the inner planet is determined by equation(16). For small amplitude librations this is given by (10). That would still apply when $`\varphi `$ is circulating provided the cosines are time averaged and the mean circulation rate is small.

## 3 A simulation of two migrating resonantly coupled planets

The protoplanetary disc is numerically simulated using an Eulerian 2D hydrodynamic code. The code used is a modified version of NIRVANA, which has been described, tested and used successfully elsewhere on a similar problem involving interacting planets (Masset & Snellgrove Masset1 (2001)). Incorparated with the hydrodynamic code is a 4th-order Runge-Kutta integrator which is used to evolve the orbits of the two planets. The gravitational forces calculated from the disk model are used in the equations of motion of the planets, and the disc itself responds to the planetary potential. Hence the system evolves in a self-consistent fashion. In order to obtain the long integration times needed for simulations of this type the FARGO algorithm (Masset Masset (2000)) is applied. However, tests have shown that the results are not affected by this.

We use a 2D cylindrical $`(r,\phi )`$ grid with 200 radial zones distributed uniformly between $`r=0.4`$ and $`r=3.47`$ in dimensionless units and 300 azimuthal zones. We apply outflow conditions at the inner boundary to simulate the accretion of disc material onto the central star.

### 3.1 Physical setup

We attempt to simulate the resonance locking of the system GJ876 via tidally induced migration of the planets, using plausible values of the disc parameters. The disc is assumed to be thin and isothermal, with constant aspect ratio $`h/r=0.07`$, and a constant Shakura & Sunyaev (Shakura (1973)) $`\alpha `$-viscosity prescription with $`\alpha =2\times 10^3`$ is adopted. The two planets are initially in circular orbits coplanar with the disc, at radial locations $`r_1=1.0`$ and $`r_2=0.6.`$ Hence the outer planet is located outside the exact $`2:1`$ commensurability ($`r_2=0.63`$). The planet masses are chosen to correspond to the minimum mass ratios obtained from observations (Marcy et al. Marcy4 (2001)). With masses normalised so that stellar mass is $`M_{}=1`$, this corresponds to $`m_1=6\times 10^3`$ and $`m_2=1.8\times 10^3`$. The planet masses are fixed as the planets are assumed to be no longer accreting material from the disc.

The disc is prescribed an initial surface density $`\mathrm{\Sigma }_0`$ corresponding to what would give a disc mass of $`2\times 10^3`$ within the orbit of the outer planet. However we assume that both planets are located inside a tidally truncated cavity located at $`r<1.3`$, with low surface density $`\mathrm{\Sigma }_{\mathrm{cavity}}=0.01\mathrm{\Sigma }_0`$. This cavity is supposed to have already been cleared by the tidal action of the two planets. Between $`1.3<r<1.5`$ the surface density is prescribed such that $`\mathrm{ln}\mathrm{\Sigma }`$ linearly joins to $`\mathrm{\Sigma }_0.`$.

### 3.2 Results

The tidal interaction of the planets with the disc material causes the planets to migrate inwards (see Figure 1). The inner planet is deep within the cavity and only interacts with low surface density material and thus migrates slowly. The outer planet has its outer $`2:1`$ lindblad resonance located outside the cavity and in the body of the disc. Hence there is more material exerting a negative torque on the planet, and therefore it migrates faster, despite its larger mass. The ratio of semi-major axes $`a_1/a_2`$ of the planets decreases until the planets ’lock’ into a $`2:1`$ commensurability with $`n_22n_1`$ at a time $`t400`$ orbits. Both planets then subsequently migrate inwards a further 10% maintaining this ratio, showing the resonance to be robust.

Figure 1 also shows the calculated values of the resonance angles $`\varphi `$ and $`\psi `$. Once the commensurability lock has occured, these are both librating about zero. The resonant interaction causes a eccentricity growth of both planets, the growth halting at around average values of $`e_1=0.06`$ and $`e_2=0.34`$ although both eccentricities exhibit variations around these average values. The perihelion angles of the two planets oscillate around the alignment position, being the natural stable state.

The disc cavity remains at low density, and the cavity edge diffuses slowly in on the viscous diffusion timescale, following the outer planet. Figure 2 shows a surface density plot of the inner regions of the disc.

### 3.3 Relation to analytic model and observations

The simulation confirms that tidally induced migration of two planets with the given mass ratios leads to trapping into the $`2:1`$ commensurability. The resonance angles both undergo stable periodic librations as per section (2) and the resonance is robust for the length of the simulation. We can conclude that stable resonance trapping would be the probable outcome of the evolution of such a system. The finding of the system GJ876 near to or in such a resonance is probably simply a consequence of past tidally induced migration of the two planets into such a state. The magnitudes of the eccentricities are in broad agreement with fitted model parameters to the observations of GJ876 with $`e_10.1`$ and $`e_20.27`$ (Marcy et al. Marcy4 (2001)). The eccentricity ratio $`e_2/e_15.7.`$ We comment that equation (11) which applies in the small eccentricity limit when one of $`\varphi ,\psi `$ is zero and the other $`\pi `$ gives $`e_1/e_21/11.`$ However, the fact that both angles librate about zero indicates, within the context of the simple model, that non resonant orbital precession needs to be incorporated and we should use equation (10) to make a comparison with the simulation. This is

$`e_1e_2a_1M_{}\omega _{pr1}+m_2a_1e_2n_1B\mathrm{cos}\varphi =`$
$`m_1a_2e_1n_2C\mathrm{cos}\psi +e_1e_2a_1M_{}\omega _{pr2}.`$ (19)
Since the librations of the angles are about a value close to zero, we replace them by zero to make simple estimates. The largest non central mass in the system is $`m_1`$ so we include its effect in causing non resonant orbital precession for $`m_2`$ by including a non zero $`\omega _{pr2}.`$ Also we include the corresponding effect on $`m_1`$ due to $`m_2`$ which produces a non zero value of $`\omega _{pr1}.`$ Then we have

$$\omega _{pr1}+\frac{m_2n_1B}{e_1M_{}}=\omega _{pr2}\frac{n_2m_1a_2C}{e_2a_1M_{}}.$$
(20)
This condition equates the precession rate of the orbit of $`m_2`$ on the right with that of $`m_1`$ on the left. In the former case the first contribution comes from the time averaged orbit and the second is the resonant contribution. In fact with $`e_20.34,`$ the orbits approach each other quite closely to within $`22`$ percent of the outer semi major axis. But the resonant configuration avoids close encounters keeping the planets apart. Thus we expect the resonant and non resonant contributions to the precession rate to show significant cancelation. The large value of $`e_2`$ probably makes strict comparison innacurate. Nonetheless the simulation gives a precession period of about $`70`$ orbits in a retrograde sense. We can estimate that this to be of the same magnitude as would be predicted from the second term on the left hand side which however gives prograde precession. Thus the non resonant precession rate due to $`m_2`$ needs to be be comparable to the resonant effect but of opposite sign.

Considering the precession rate of $`m_2,`$ we estimate that $`\omega _{pr2}`$ corresponds to $`70`$ orbits in the prograde sense while the resonant term corresponds to $`70/2`$ orbits in the retrograde sense. Thus the combined effect produces a precession period of $`70`$ orbits in the retrograde sense as seen in the simulation.

In summary our simulation gives plasuble eccentricity values for the two planets, that can be understood in outline by use of a simplified analytic theory and are consistent with the current observations.

## 4 Orbit integrations

In addition to the analytic model presented in section (2) and the hydrodynamic simulations presented in section (3), we have also performed three-body orbit integrations using a fifth-order Runge-Kutta scheme (e.g. Press et al. 1993).

The basic assumptions of the model are that the two planets exist within the inner cavity of a tidally truncated disc that lies exterior to the outer planet. Tidal interaction with this disc causes inwards migration of the outer planet, and also leads to eccentricity damping of the outer planet. It is further assumed that as the planets migrate inwards and approach their final semi-major axes, the disc disperses on some prescribed time scale $`t_{disp}.`$ In our numerical calculations, a torque was applied to the outermost planet such that it migrated inwards on a time scale of $`t_{mig}`$ local orbital periods as defined in section(2), and a damping force was applied in the radial direction to damp the eccentricity on a time scale of $`t_c`$ local orbital periods also as defined in section (2).

These integrations used initial conditions corresponding to the more massive, outermost planet being located initially at 5 AU, with the lighter inner-most planet located initially at 2.5 AU. The planet masses adopted for the orbit integrations are the same as the minimum masses reported for the planets in the system around GJ876 by Marcy et al. (Marcy4 (2001)) (i.e. 1.87 M<sub>J</sub> and 0.56 M<sub>J</sub>). The stellar mass is taken to be 0.32 M. Whilst these calculations provide only a crude approximation to the detailed physics of disc–companion interactions, their simplicity allows us to perform many calculations, covering a wide area of parameter space, and also to run for much longer time scales than is possible for simulations of the type described in section(3).

A number of calculations have been performed to examine the relationship between the final values of $`e_1`$, $`e_2`$, and their ratio $`e_1/e_2`$, to the various input parameters $`t_{mig}`$, $`t_c`$, and $`t_{disp}`$. The results of some of these calculations are presented in table 1, and are discussed below. The unit of time used in the abscissa of the figures 3 to 5 is the orbital period for an object at 1 AU in orbit around a star with mass $`0.32`$ M, and is denoted as P(1 AU).

### 4.1 Dependence on migration and circularization times

Equation 18 shows that the eccentricity of the outer planet, $`e_1`$, depends on the ratio of $`t_c/t_{mig}`$. Here we present results of simulations that explore how the eccentricity ratio $`e_2/e_1`$ depends on $`t_c`$ and $`t_{mig}.`$

Figure 3 shows the evolution of the semi–major axes and eccentricities for the run R1, whose model parameters are described in table 1. This figure shows the inward migration of the outer planet that subsequently locks to the inner planet as it reaches the 2:1 commensurability. The subsequent evolution is such that the two planets, now resonantly locked, migrate inwards. The eccentricities result from the balance between eccentricity driving through the resonant interaction, and eccentricitiy damping due to the disc interaction. As the planets approach their final semi-major axes, the effects of migration and eccentricity damping are removed on a time scale of $`t_{disp}=8\times 10^3`$ local orbits, causing them to cease migration at semi–major axes $`a_10.2`$ and $`a_20.126`$ which are values appropriate to GJ876. The subsequent evolution beyond a time of $`t0.45`$ P(1 AU) in figure 3 occurs in the absence of disc effects, and suggests a long–term stability of the system given that it remains stable for $`2\times 10^7`$ orbits of the outer planet at its final semi-major axis. Figure 4 shows the evolution of the resonant angles $`\varphi `$ and $`\psi `$, which librate about $`\varphi =\psi =0`$ in agreement with the results presented in section 3 for the full hydrodynamic simulations.

Figure 4 shows the results from run R2 described in table 1, and illustrates the effect of reducing $`t_{mig}`$ while keeping $`t_c`$ constant. As expected from equation 18, the eccentricity of the outer planet increases from $`e_10.1`$ to $`e_10.35`$ when $`t_{mig}`$ changes from $`3.33\times 10^3`$ to $`3.33\times 10^2`$ local orbital periods. However, the ratio $`e_2/e_1`$ has also changed from $`e_2/e_14`$ to $`e_2/e_12`$.

Calculation R3 illustrates the effect of keeping $`t_{mig}`$ constant while changing $`t_c`$. As expected from equation 18, reducing $`t_c`$ leads to a reduction in $`e_1`$, since the disc model damps the eccentricity more effectively. We also find that the ratio $`e_2/e_1`$ again changes, it now being $`e_2/e_16`$, as compared to $`e_2/e_14`$ for run R1.

Equation 18 predicts that keeping the ratio $`t_{mig}/t_c`$ constant, but changing both $`t_{mig}`$ and $`t_c`$ independently, should leave $`e_1`$ unchanged. Calculation R4 indicates that this is what happens, and also shows that the ratio $`e_2/e_1`$ remains unchanged.

Overall the results are entirely consistent with the analytic predictions presented in section 2 and with the hydrodynamical simulations presented in section 3. Furthermore, they indicate that the ratio $`e_2/e_1`$ scales rather weakly with $`t_{mig}/t_c`$. We comment that from section 2 equation 10 we expect the eccentricity ratio to reach $`11`$ as $`e_10`$. Long–term stabilty of two–planet systems that become locked due to disc–induced orbital migration is also indicated by our calculations. In particlular run R1 covers a time scale corresponding to $`2\times 10^7`$ orbits of the outer planet in its final configuration. We find that it is possible to arrange $`e_1`$ to match the observed value of the outer planet in the GJ876 system by fixing $`t_{mig}`$ and choosing $`t_c`$ appropriately, but it is difficult to then obtain a value for $`e_2`$ that matches the reported value of $`e_2=0.27`$.

### 4.2 Dependence on Disc Dispersal Time Scale

We have performed simulations to examine the effect of removing the disc on different time scales on the stability of the system. We find that the stability is largely unaffected by the rate at which the disc is removed. Figure 3 shows the evolution of the semi-major axes and eccentricities from a run in which the disc was removed on a time scale of $`t_{disp}=8\times 10^3`$ local orbital periods. In this case the disc dispersal was switched on once the outer planet semi-major axis was $`a_1<0.5`$ AU. The final semi-major axis of the planet is determined by the radius at which the disc dispersal is initiated, and the time scale over which the disc is removed. Calculation R5 was similar to R1 except that the disc was removed instantaneously. The results presented in table 1 show that this has little effect on the final outcome. A slight increase in the scatter of the temporal distribution of eccentricities was observed.

## 5 Summary and Discussion

In this paper we have considered two protoplanets gravitationally interacting with each other and a protoplanetary disc. The two planets orbit interior to a tidally maintained disc cavity while the disc interaction induces inward migration.

We have supposed that the previous evolution of the system results in the the planets getting into a configuration with an orbital separation just exceeding that required for a $`2:1`$ commensurability. This evolution is likely to have involved both accretion and migration. Given the set up we consider a natural evolution is towards both planets orbiting in a cavity outside of which orbits the protoplanetary disc. Tidal interaction results in the inner disc material either being expelled into the outer disc or accreting onto the central star. Subsequently the outer planet migrates towards the inner one as a result of interacting tidally with the exterior disc.

When the migration is slow enough, we found that the outer protoplanet approached and became locked into a $`2:1`$ commensurability with the inner one. This was maintained in subsequent evolution. We studied the nature of these interactions using a simple analytic model, hydrodynamic 2D simulations and longer time N body integrations. These all gave consistent results.

The magnitude of the stabilized eccentricities was found to be determined by the ratio of the migration rate to the circularization rate induced in the outer planet orbit by the external disc. The eccentricity ratio $`e_1/e_2`$ was found to vary with the magnitude of the ratio being more extreme for smaller eccentricities.

We have applied our results to the recently discovered resonant planets around GJ876. Simulation shows that a disc with parameters expected for protoplanetary discs causes trapping in the $`2:1`$ commensurability when the planets orbit in an inner cavity and that eccentricities in the observed range may be obtained. In this case the orbits were found to be aligned with both resonant angles librating about zero. The whole system then precessed in a retrograde sense. Finally when the disc is removed on a range of timescales the orbital configuration has been found to be stable for up to $`2\times 10^7`$ orbits subsequently.

###### Acknowledgements.

This work was supported by PPARC grant number PPARC GR/L 39094. We thank Udo Ziegler for making a FORTRAN Version of his code NIRVANA publicly available.


format: text

text:
# Itinerant Ferromagnetism for Mixed Valence Systems

## Abstract

We introduce a novel mechanism for the unusual itinerant ferromagnetism found in mixed valence systems like Ce(Rh1-xRux)3B2 , LaxCe1-xRh3B2 , US, USe, and UTe. With it we can provide an explanation for the long-unexplained large value of $`T_c`$ ($``$ 100K) value and the maximum in the magnetization below $`T_c`$ found experimentally. We also show that this novel itinerent ferromagnetism can be continuously connected with the localized case for which the energy scale is much smaller ($`J_{RKKY}`$ 1K).

Introduction. Understanding the mechanisms for ferromagnetism usually involves a reconciliation of a localized electron picture, traced back to Heisenberg , and an itinerant electron picture, traced back to Bloch . For many insulating materials, adding itinerant features, such as indirect exchange, to the localized picture brings satisfactory agreement with basic experimental features. For many metals, adding localized features, such as spin waves, to the itinerant picture has a similar effect. On the other hand it is now appreciated that for some novel classes of materials, such as the heavy Fermion and mixed valence materials, doctoring one picture or the other is often a questionable procedure. The electrons can display both localized and itinerant properties as these materials are neither good insulators nor good metals.

In this letter we propose a simple mechanism for ferromagnetism (FM) in heavy fermion mixed valence materials. The mechanism relies on several specific energy scales set by the band structure of these materials even though the electron correlation energy is the dominant energy scale. The simplicity and generality of the mechanism however should provide a useful framework for discussing ferromagnetism in broad classes of heavy fermion mixed valence materials often found in 4$`f`$ and 5$`f`$ materials and their compounds. In fact we will show that it seems to explain some of the very unusual long unexplained ferromagnetic properties of such Ce based compounds as Ce(Rh1-xRux)3B2  and the uranian monochalcogenides .

The mechanism is based upon a two-band structure. The lower band near the zone center is dispersive and nearly free of electron correlation effects. Away from the zone center, this band is very flat and doubly occupied states experience strong correlation effects. The upper band is separated from the lower band by a gap that is smaller than the correlation energy. When the chemical potential is at the flat band level, the system becomes unstable to the formation of a ferromagnetic state. There are compounds satisfying this mixed valence requirement under normal conditions or by either adding additional electrons (doping), applying pressure, or changing the temperature. This ferromagnetic state is a manifestation of a Hund-like rule among electrons in band states as opposed to usual case of electrons in localized states (atomic orbitals).

We will show that the periodic Anderson model (PAM) is a simple microscopic realization of this picture. Using quantum Monte Carlo simulations, we will demonstrate that the model admits a ferromagnetic ground state in the mixed valence regime with electron occupancies reflecting our physical picture. We then develop a mean field picture that reproduces the ground state properties of the PAM and also allows us to compute finite temperature properties. Under certain conditions we find the finite temperature mean-field approximation predicts the same highly unusual behavior as a peak, below $`T_c`$, in the temperature dependence of the magnetization  and the deviation of the inverse susceptibility from Curie-Weiss behavior .

Model Hamiltonian. The PAM is described by the Hamiltonian

$`H`$ $`=`$ $`t{\displaystyle \underset{i,j,\sigma }{}}(d_{i\sigma }^{}d_{j\sigma }+d_{j\sigma }^{}d_{i\sigma })+V{\displaystyle \underset{i,\sigma }{}}(d_{i\sigma }^{}f_{i\sigma }+f_{i\sigma }^{}d_{i\sigma })`$ (2)
$`+ϵ_f{\displaystyle \underset{i,\sigma }{}}n_{i\sigma }^f+{\displaystyle \frac{U}{2}}{\displaystyle \underset{i,\sigma }{}}n_{i\sigma }^fn_{i\overline{\sigma }}^f`$
where $`d_{i\sigma }^{}`$ and $`f_{i\sigma }^{}`$ create an electron with spin $`\sigma `$ in $`d`$ and $`f`$ orbitals at lattice site $`i`$ and $`n_{i\sigma }^f=f_{i\sigma }^{}f_{i\sigma }`$ is the number operator for the $`f`$-electrons of spin $`\sigma `$ at site $`i`$. The hopping amplitude $`t`$ between $`d`$-orbitals is only to nearest-neighbor sites. The hopping amplitude $`V`$ hybridizes different orbitals on the same site.

When $`U=0`$, the resulting non-interacting Hamiltonian $`H_0`$ has two dispersive bands:

$$E_\sigma ^\pm (𝐤)=\frac{1}{2}\left[e_𝐤+ϵ_f\pm \sqrt{(e_𝐤ϵ_f)^2+4V^2}\right]$$
(3)
separated by an hybridization gap. For a chain, the band structure is illustrated in Fig. 1a. Irrespective of spatial dimension, each band has regions where the $`𝐤`$-states are dominantly of $`d`$ or $`f`$-character. A mixed valence condition arises when the Fermi energy $`E_F`$ sits in a crossover region which is around the energy $`ϵ_f`$. Three energies are important: the hybridization gap $`\mathrm{\Delta }`$, the $`f`$-state dispersion $`\delta _f=E^{}/k|_{k_F}=\mathrm{}v_F`$ ($`v_F`$ is the Fermi velocity), and the Coulomb (correlation) energy $`U`$. $`\delta _f`$ is measuring the band dispersion close to the Fermi level. In the interacting problem we will assume that $`U>\delta _f,\mathrm{\Delta }`$. As we shall now discuss, the competition is between $`\delta _f`$ and $`\mathrm{\Delta }`$ which will be a competition between a ferromagnetic and paramagnetic state.

Mechanism for Ferromagnetism. In Fig. 1a-c we illustrate the physical mechanism for the ferromagnetic state which emerges when the electron filling $`\rho `$ is arround 1/4 (one electron per lattice site) and the energy $`ϵ_f`$ of the $`f`$-orbitals is close to $`E_F`$. The latter condition defines the intermediate valence regime.

To explain the mechanism, we first consider the non-interacting case (Fig. 1a). States in the lower band with mainly $`f`$-character have a small dispersion $`\delta _f`$ because of the absence of direct hopping between the $`f`$-orbitals. In the interacting case, because the Coulomb interaction only affects the electrons in the $`f`$-orbitals, the electrons which were doubly occupying the states in the lower band with mainly $`d`$-character are practically unaffected. On the other hand, the electrons in the many states which are close to the Fermi level and have mainly an $`f`$-character will be strongly affected. These electrons spread to higher unoccupied $`𝐤`$-states in the $`f`$-part of the band and polarize by an in-band Hund’s rule (Fig. 1b). To see the analogy to Hund’s rule more clearly, we consider the limiting case, represented in Fig. 1c, where $`\delta _f=0`$ and $`\mathrm{\Delta }0`$, reducing Fig. 1b to a two level system in momentum space with each level being strongly degenerate (we only plot the states with predominant $`f`$ character). If we add electrons to this two levels system it is easy to show that the ferromagnetic solution has the lowest energy and is therefore the ground state. By polarizing, the spatial part of the wave function becomes antisymmetric and there is no double occupancy in the real space. In this way, the Coulomb repulsion is reduced to zero. The kinetic energy has also the lowest possible value if the electrons are occupying the lower energy (degenerated) levels (see Fig. 1c).

We now describe the conditions for the stability of the ferromagnetic state, and we will see how this mechanism departs from the one of Stoner . To this end we state the condition $`U\mathrm{\Delta }`$. From Fig. 1b it is clear that the cost of the spreading of the ferromagnetic solution is an increase in the kinetic energy proportional to $`\delta _f`$. On the other hand, if we build up a nonmagnetic solution only using the states of the lower band, there will be a minimum delocalization for each particle due to the fact that a finite set of $`f`$-character states is in the upper band (see Fig. 1). This can be seen by constructing an $`f`$ Wannier function with the lower band states. This Wannier function will have a delocalization lenght which depends on the $`𝐤`$ wave vector where the two bands are crossing. Therefore, by localizing in real space to avoid the double occupancy, the electrons will have a finite probability of occupying the $`f`$-states in the upper band. The energy cost per electron of occupancy is proportional to $`\mathrm{\Delta }`$. If the hybridization gap is much larger than the dispersion of the states, that is, $`\mathrm{\Delta }\delta _f`$, the ferromagnetic state lies lower in energy than the nonmagnetic state. The energy of the excited nonmagnetic state is proportional to $`\mathrm{\Delta }`$. This can be easily seen from the limiting (and non-realistic) case $`\delta _f=0`$ ploted in Fig. 1c, for which the analogy with the atomic shells and the Hund’s rule is evident.

It is important to emphasize that the hybridization gap $`\mathrm{\Delta }`$ and $`f`$-state dispersion $`\delta _f`$ are the two basic ingredients of our mechanism: $`\mathrm{\Delta }`$ only appears if there is more than one band and $`\delta _f`$ is small only if the hybridization is weak. If $`U>\mathrm{\Delta },\delta _f`$, then $`\mathrm{\Delta }`$ is the energy scale of the paramagnetic state, while $`\delta _f`$ is the energy scale of the ferromagnetic solution. We can estimate the magnitude of these scales for the PAM: for $`Vt/2`$, we get $`\delta _fe_k/k|_{k_F}V^2/(e_{k_F}ϵ_f)^2`$ and $`\mathrm{\Delta }V^2/(2dtϵ_f)+\sqrt{(2dt+ϵ_f)^2+4V^2}2dt`$. Therefore, if $`ϵ_f`$ is close to the bottom of the lower band ($`ϵ_f+2dtt`$) and the Fermi level is close to the top ($`E_Fϵ_fV^2/(2dtϵ_f)`$), $`\mathrm{\Delta }`$ is considerably larger than $`\delta _f`$ and the ground state is ferromagnetic. We can see that the ferromagnetic solution is stable for comfortably realistic values of the parameters. It is clear from this analysis the particle density $`\rho `$ for a ferromagnetic solution must be close to quarter filling.

Quantum Monte Carlo Method. Using the constrained-path Monte Carlo method (CPMC), we computed the ground-state properties of the PAM on a square lattice. The CPMC method projects the ground state from an initial state $`|\psi _T`$ by converting the iterative procedure

$$|\psi _{i+1}=e^{\tau H}|\psi _i$$
(4)
into a branched random walk. The details of the method are described elsewhere  as are our prior uses of it on the square PAM .

The defining characteristic of the CPMC method is its elimination of the fermion sign problem by excluding random walkers $`|\varphi `$ that violate $`\psi _T|\varphi >0`$. If $`|\psi _T`$ were the exact ground state $`|\psi _0`$, this constraint would generate an exact elimination of the sign problem, and hence an exact solution. Extensive benchmarking indicates the CPMC method provides accurate estimates of the energy and various correlation functions.

We prepared the trial state in specific values of the total spin and z-component of spin . Because the Hamiltonian(2) conserves these quantum numbers, the iterative process (4) produces a ground state with the same $`S`$ and $`S_z`$. In the ferromagnetic regime, when the number of lattice sites $`N`$ was greater than the number of electrons $`N_e`$, the resulting energy as a function of $`S`$ almost always showed a minimum at a value of $`S`$ satisfying $`0<S<\frac{1}{2}(2NN_e)`$. Hence we typically found a partially polarized ferromagnetic ground state.

For square lattices we also computed dependence of the static spin structure factor and the electron occupancies both on position and wave number space . Most notable for the present work is the wave vector dependence of the different spin components of the upper and lower band electron occupancies shown in Fig. 2. It has the same quantitative features as Fig. 1b.

Mean Field Theory. The quantum Monte Carlo simulations are limited to chains and square lattices of relatively small sizes at zero temperature. However, we have found that their predictions are well described by a simple spin-polarized mean-field approximation. We can easily extend this approximation to larger systems sizes in higher dimensions and at finite temperatures. We can use large lattices in three-dimensions at finite temperatures to compare with the results of experiments.

The mean field Hamiltonian $`H_{MF}`$ is

$$H_{MF}=H_0+\frac{U}{2}\underset{i,\sigma }{}(n_{i\sigma }^fn_{i\overline{\sigma }}^f+n_{i\sigma }^fn_{i\overline{\sigma }}^fn_{i\sigma }^fn_{i\overline{\sigma }}^f)$$
(5)
$`n_i^f`$ and $`n_i^f`$ are the determined in a self-consistent way within the mean-field theory. The translational symmetry of $`H`$ implies that we can always find a translationally invariant ground state. This means that $`n_i^f`$ and $`n_i^f`$ will not depend on the site index $`i`$. Therefore there are two variational parameters to be determined from the mean-field equations: $`n_{}^f`$ and $`n_{}^f`$.

Even though the mean field approximation can give a good description of the FM ground state, it is well known that such approximations overestimate the energy of the paramagnetic phase because they improperly estimate the real space correlations that are very important for the paramagnetic solution. Therefore it is crucial to check that the ground state is the FM one by using a more accurate method to evaluate the energy of the PM state. To this end we calculated the ground state of the PAM with $`6\times 6`$ unit cells using the CPMC technique. We found excellent agreement with the energy and as seen in Fig. 2 the electron occupancies.

Consequences. Bolstered by this agreement we used the mean-field approximation to calculate various thermodynamic properties of the ferromagnetic phase for the three-dimensioanl cubic PAM. In Fig. 3a we show the temperature dependence of the magnetization for different values of $`ϵ_f`$. From the arguments given above it is clear that the magnetization comes from the electrons which are occupying the states with $`f`$-character. By increasing $`ϵ_f`$ we are reducing the number of $`f`$-electrons and therefore the magnitude of the zero temperature magnetization. Some critical value of $`ϵ_f`$ exists where the chemical potential starts departing from $`ϵ_f`$ and a new energy scale $`ϵ_fE_F`$ emerges in the system. This new energy scale is reflected in the appearance of a magnetization peak. (See the $`ϵ_f=t`$ case in Fig. 3a.) To understand the last statement one first has to realize that the zero temperature magnetization is small due to the reduction of the number of $`f`$-electrons. When the temperature is of the order of $`ϵ_fE_F`$, electrons are promoted electrons from the double occupied band states to the $`f`$-character states which have a large entropy (large density of states). The $`f`$-electrons will be polarized due to the energy considerations above explained. In this way we can explain the origin of the magnetization peak below $`T_c`$. It is important to note that the source for the large entropy is associated with charge and not with spin degrees of freedom. This fact explains why a state with larger magnetization has a higher entropy. The magnetization curves shown in Fig. 3a are in good qualitative agreement with the magnetization versus temperature data measured in Ce(Rh1-xRux)3B2 for different values of $`x`$ .

Another interesting aspect of this ferromagnetism is the deviation from linearity (Curie-Weiss behavior) for the inverse magnetic susceptibility above $`T_c`$ over a large temperature range (see Fig. 3b) . The behavior contrasts that predicted by theories based on localized magnetic moments.

From above it is also clear that the change in entropy from the magnetic to the paramagnetic phase depends on the number for $`f`$ electrons (spin degrees of freedom). This dependence is seen in Fig. 3d where we plot the calculated specific heat for different values of $`ϵ_f`$.

A characteristic of this ferromagnetic state that has consequences for the photoemission experiments is the implication that the quasi-particle dispersion should be close to that of the electrons in the non-interacting case. We recall that the ferromagnetic state described in Fig. 1 is very similar to a non-interacting polarized state. As in the Stoner mechanism, the main role of the Coulomb repulsion is to polarize states with well defined momentum. The hybridization gap of the noninteracting solution is most likely replaced by a pseudogap (see Fig. 3c). This pseudogap should be seen in the optical conductivity measurements.

Conclusions. Our quantum Monte Carlo simulations and mean-field calculations clearly show the existence of an itinerant ferromagnetic phase in mixed valence materials that is supported by the PAM. This new phase describes qualitatively many experimental features such as an unusually large value for $`T_c`$($``$ 100K) and the maximum in the magnetization below $`T_c`$  that are found in such ferromagnetic compounds as Ce(Rh1-xRux)3B2 , Ce(Rh1-xOsx)3B2 , and LaxCe1-xRh3B2 .

Many previous theories of 4$`f`$ and 5$`f`$ electron materials treated these materials as systems of localized moments in the $`f`$-orbitals. These theories are thus unable to describe an itinerant ferromagnetic phase and experimental consequences  like the peak in the magnetization below $`T_c`$ observed in some of these systems , the large value of $`T_c`$, the deviation of $`\chi ^1(T)`$ from the Curie-Weiss law above $`T_c`$ in the uranium monochalcogenides , and the mixed valence behavior of these compounds .

The physical picture just presented, when combined with our previous results , allows a reconciliation of the localized and delocalized ferromagnetism pictures painted by Heisenberg and Bloch seventy years ago in the sense that it is possible to go continuously from the mixed valence FM state where the $`f`$-electrons are delocalized to a FM state where there is one localized electron in each $`f`$-orbital . In our picture one can do this by decreasing $`ϵ_f`$ from the Fermi level to a value near or below the bottom of the valence band. We note that the energy scale of a localized ferromagnetic state is that of the RKKY interaction, which according to de Gennes’s rule  is of the order of 1K, while the scale for the itinerant case introduced here is that of the hybridization gap $`\mathrm{\Delta }`$, which is of the order of 100K for heavy fermion compounds . In going from the delocalized (mixed valence) to the localized regime we thus expect a strong reduction of $`T_c`$ accompanied by an increase of the zero temperature magnetization. (The $`T=0`$ magnetization is proportional to the $`f`$ occupation number.) This expectation is consistent with the observed behavior of the magnetization in LaxCe1-xRh3B2 as function of $`x`$ .

Acknowledgements. This work was sponsored by the US DOE under contract W-7405-ENG-36. We acknowledge useful discussions with A. J. Arko, B. Brandow, M. F. Hundley, J. J. Joyce, J. M. Lawrence, S. Trugman, and J. L. Smith. In particular, we thank J. M. Lawrence for pointing out the experimental work on the Ce compounds. J. B. acknowledges the support Slovene Ministry of Education, Science and Sport, Los Alamos National Laboratory and FERLIN.


text:
# Laying the groundwork at the AGS: Recent results from experiment E895

## 1 Introduction

Modern relativistic heavy ion physics focuses largely on the creation and study of the Quark Gluon Plasma (QGP), a state of matter in which partonic degrees of freedom are relevant over large length scales; potential signatures of QGP formation (see, e.g. ) are studied in collisions at AGS, SPS, and now RHIC. However, it is now well-accepted that only through a simultaneous study of several observables will a clear understanding of the collision history emerge. Furthermore, the most convincing evidence of “new” physics stems from a break in well-understood (in terms of hadronic models) systematics in observables which have been measured over a large range of conditions.

Using the EOS Time Projection Chamber (TPC)  at the Brookhaven AGS, we have studied Au+Au collisions at bombarding energies of 2, 4, 6, and 8 AGeV. Models suggest that collisions at these energies reach the highest baryon density , making such systems ideal to stringently test our understanding of nuclear matter under extreme conditions, as well as to search for possible QGP formation in the high-density limit .

We have measured excitation functions in four primary areas: (1) proton elliptic and directed flow , (2) intensity interferometry , (3) charged particle spectra , and (4) strange particle production . In addition to improving our understanding of hot hadronic systems under extreme conditions, these systematics provide a valuable baseline for measurements at higher energy.

Going further, we have combined the analyses in these areas to achieve novel insights into the collision dynamics, and serve as an example of how to optimize physics information extraction for data at RHIC. In , strange particle reconstruction is correlated with proton flow to extract information on in-medium potentials. In this report, we focus primarily on combined analyses involving two-particle correlations.

## 2 Directed and Elliptic Flow

Anisotropic flow refers to the momentum distribution of particles in azimuth about the beam axis. The reference azimuthal angle is the direction of the impact parameter vector $`𝐛`$, which must be determined on an event-by-event basis; see  for details.

Two types of anisotropic flow are typically discussed . (Positive) sideward flow refers to the preferential in-plane emission of particles with $`y>y_{cm}`$ in the direction of $`𝐛`$, and those with $`y<y_{cm}`$ opposite $`𝐛`$, while elliptical flow analyses focus on the difference between in- and out-of-plane emission. In a simple picture in which the emission pattern can be described by an ellipsoid in momentum space (in general, the situation is more complicated ), then, sideward flow analyses measure the tilt angle of the momentum-space ellipsoid, while elliptical flow probes the ratio between its semi-minor axes.

E895 analyses of proton directed and elliptic flow have already been reported . There, we found that AGS energies represent a transition region, in which the elliptic flow changes sign from out-of-plane (the “squeeze-out” signal observed at Bevalac/SIS energies) to in-plane elliptic flow seen at SPS energies, and directed flow reaches a maximum and falls off quickly with beam energy. Pion directed flow is oriented in the direction opposite proton flow– i.e. the $`\pi `$ momentum-space ellipse tilts away from $`𝐛`$ for $`p_z>0`$ , while elliptic flow changes sign at $`E_{beam}`$4 AGeV, similar to the protons .

Dynamic transport models reproduce the gross trends in flow measured by E895, but the details are strongly sensitive to the underlying physics (e.g. the Equation of State) assumed ; hence, these excitation functions provide a very stringent test of our understanding of the dynamics of heavy ion collisions. In particular, no single model is currently able to reproduce the detailed flow excitation function at AGS energies.

Clearly, fresh information on the physics of anisotropic flow is needed. It is commonly assumed that these flow signals result from space-momentum correlations generated in the collision. By measuring azimuthally-sensitive momentum distributions, therefore, we look at only half the story. In the next section, we describe how two-pion correlations provide a more complete picture.

## 3 Pion Intensity Interferometry

Intensity interferometry (a.k.a. HBT) techniques have been used extensively to probe the space-time structure of heavy ion collisions (see, e.g. ). Until recently, most HBT analyses implicitly assumed a cylindrically symmetric emission region, strictly valid only for $`|𝐛|=0`$ collisions. In this standard framework, E895 $`\pi ^{}`$ correlation systematics as a function of beam energy, $`m_T`$, centrality, and rapidity, have already been reported .

In addition to the (momentum-space) anisotropic flow studies discussed above, theoretical studies  suggest that coordinate-space anisotropies of the emitting source are equally interesting. Recently, a formalism was proposed to use azimuthally-sensitive intensity interferometry, correlated event-by-event with the reaction plane, to extract this information experimentally . We have performed the first experimental measurement of the full coordinate-space anisotropies of a hot nuclear source, for pions emitted from Au+Au collisions at 2-6 AGeV. Full details may be found in .

$`\pi ^{}`$ correlation functions were constructed in the Bertsch-Pratt (“out-side-long”) decomposition as measured in the Au+Au c.m. frame , with cuts on the pair angle with respect to the reaction plane $`\varphi =\mathrm{}(𝐊_{},𝐛)`$, where $`𝐊_{}=(𝐩_\mathrm{𝟏}+𝐩_\mathrm{𝟐})_{}`$ is the total momentum of the pair perpendicular to the beam.

The correlation functions for each of eight 45-degree-wide $`\varphi `$ bins, are fit with the standard Gaussian parameterization

$$C(𝐪,\varphi )=1+\lambda (\varphi )\mathrm{exp}\left[\underset{i,j=o,s,l}{}q_iq_jR_{ij}^2(\varphi )\right].$$
(1)
In contrast to $`\varphi `$-integrated analyses, all six HBT radii are relevant here .

For 4 AGeV collisions, the seven fit parameters $`R_{ij}^2`$ and $`\lambda `$ are plotted as a function of $`\varphi `$ in Figure 2. The $`\varphi `$-independence of $`\lambda `$ indicates that the fraction of $`\pi ^{}`$ from long-lived resonances is independent of emission angle with respect to the reaction plane, consistent with the relatively weak anisotropic flow at these energies , as well as the randomizing effect on the pion from decay of the parent. $`R_o`$ and $`R_s`$ show significant equal and opposite $`2^{nd}`$-order oscillations in $`\varphi `$, while $`R_{ol}^2`$ and $`R_{sl}^2`$ display equal-magnitude $`1^{st}`$-order oscillations.

From the azimuthal dependence of the HBT radii, one may extract the geometry of the source in the impact-parameter-fixed coordinate system . In particular, at midrapidity and low $`p_T`$, $`R_{ij}^2(\varphi )`$ are sensitive to the homogeneity lengths in x, y, z, and time, as well as the tilt of the source in the reaction plane. Thus, we determine the shape and orientation of the freeze-out ellipsoid in coordinate space. The lines in Figure 2 represent a simultaneous fit to all six $`R_{it}^2(\varphi )`$, using the formalism of , resulting in the following homogeneity lengths: $`\sigma _x=4.0\pm 0.1`$ fm, $`\sigma _y=5.2\pm 0.1`$ fm, $`\sigma _z=5.2\pm 0.1`$ fm, $`\sigma _t=3.7\pm 0.5`$ fm/c; the extracted tilt angle is $`\theta _s=37^{}\pm 4^{}`$. Finite reaction-plane resolution effects, which wash out oscillations in $`R_{ij}^2(\varphi )`$, are accounted for .

Similar shapes and orientations were found for collisions at 2 and 6 AGeV. At all energies, our results indicate a pion freezeout distribution as an ellipsoid whose major axis in the reaction plane is tilted with respect to the beam in the positive direction (i.e. in the direction of $`𝐛`$), and whose transverse axis perpendicular to the reaction plane is longer than the axis in the reaction plane. The extension in the temporal direction is consistent with observations at high energy .

The coordinate-space analog of elliptical flow, $`\sigma _y/\sigma _x`$, is shown as a function of beam energy in Figure 2. In contrast to the momentum space trends, there is no transition from out-of-plane to in-plane shape. Indeed, the observed transverse shape is reminiscent of the entrance channel geometry. In the simplest picture, pions are emitted from the overlap region of the two spherical Au nuclei. For impact parameter $`b=48`$ fm, the ratio of out-of-plane to in-plane RMS for this overlap region varies from 1.3-1.9 (indicated by dashed lines in Figure 2), consistent with the measured value $`1.3`$.

The tilt angle, shown in the lower panel of Figure 2, is the coordinate-space analog of directed flow, and is a geometric feature of the collision dynamics . Studies with the RQMD model indicate that most of the low-$`p_T`$ pions arise from $`\mathrm{\Delta }`$ decay, so it is not surprising that the tilted freeze-out distribution resembles the baryonic distribution in coordinate space calculated in transport codes .

Experimental access to this level of geometric detail on the freezeout distribution is unprecedented, and represents an exciting new opportunity to study the dynamical response of hot nuclear matter to compression. An azimuthally-sensitive interferometry analysis  on pions generated by the RQMD  model displays considerable sensitivity on the dynamical effect of the nuclear meanfield. Although at these energies the model underpredicts the spatial and overpredicts the temporal scale (noted already for central collisions ), qualitatively, the transverse shape and large positive tilt angles are reproduced remarkably well; transverse asymmetries and spatial tilt angles are shown in Figure 2. Since RQMD better describes proton directed flow (momentum-space tilt) when the meanfield is included in the calculation , it is interesting to note that the spatial tilt angles are better described when the meanfield is off. Clearly, probing asymmetries in both momentum and coordinate space results in a more stringent test of the dynamics.

These coordinate-space anisotropies represent new and independent information. Momentum-space tilt angles (flow angles ) at these energies are only a few degrees , and, as mentioned above, for pions point in the negative direction , opposite the coordinate-space tilt $`\theta _s`$. Experimental information on the interplay between coordinate and momentum space anisotropies should help resolve theoretical issues, such as the coexistence of flow and antiflow components .

## 4 Six-dimensional Pion Phasespace Density

We have measured an excitation function of charged particle spectra over a wide region of phasespace. Figure 3 shows $`m_T`$ distributions for $`\pi ^{}`$ from central Au+Au collisions . Assuming a thermal picture, these spectra may be combined with those of heavier particles, to extract freezeout temperatures and radial flow values .

As a test of this thermal assumption, it has been observed  that for a thermally equilibrated pion source, both the shape and the amplitude of the six-dimensional phasespace occupation density as a function of spatial and momentum coordinates, $`f(x,p)`$, are determined by the temperature. Experimentally, we may probe the spatially-averaged phasespace density by combining spectral and interferometry information, according to

$$f(p_T,y)=\sqrt{\lambda (p_T,y)\pi ^3}\frac{\frac{1}{E}\frac{d^3N}{dyp_Tdp_Td\varphi }}{R_o(p_T,y)R_s(p_T,y)R_l(p_T,y)},$$
(2)
(For noncentral collisions or away from midrapidity, “cross-term” HBT radii also enter the denominator.) The factor $`\sqrt{\lambda }`$ from the HBT analysis is included in Eq. 2 to account for pions from long-lived decays .

For heavy ion collisions at maximum AGS and SPS energies, a remarkable “universal” phasespace density $`f(p_T,y)`$ has been observed , independent of collision energy or system, though multiplicities vary by more than an order of magnitude. This suggests that $`\pi \pi `$ cross-sections dominate the freeze-out density of high-energy collisions, and it is natural to look for a break-down of the universal behavior at lower energy.

Figure 5 shows the $`\pi ^{}`$ average phasespace density $`f(p_T,y=0)`$ determined according to Equation 2 for central Au+Au collisions at 2-8 AGeV. Clearly, the universal behavior has broken down; $`f`$ grows monotonically with collision energy at the AGS. As at higher energies, $`f1`$, indicating that multiparticle symmetrization effects are weak . Lines indicate expectations for “pure thermal” (no flow) sources with $`T100`$ MeV; The extracted values appear roughly consistent with the shape and magnitude of a thermal source, with a temperature increasing with $`E_{beam}`$. In detail, however, the observed $`f(p_T)`$ is flatter than the “pure thermal” expectation, consistent with a high-$`p_T`$ enhancement due to radial flow ; preliminary fits (not shown) with a parameterization that includes radial flow  reproduce the distributions well, and indicate that both temperature and flow increase with beam energy.

Also shown on Figure 5 are values from Pb+Pb collisions at 158 AGeV  and preliminary data from STAR (Au(130 AGeV)+Au(130 AGeV)) , demonstrating the persistence of the universal behavior to RHIC energies. While the phasespace densities achieved at SPS and RHIC are higher than those at AGS at higher $`p_T`$, presumably due to the greater explosive flow generated, it is remarkable that the “universal” value at low $`p_T`$ is reached already for collisions below maximum AGS energy.

## 5 $`\mathrm{\Lambda }`$-p Correlations

In the baryon-rich system generated at the AGS, it is interesting to study the baryon freezeout geometry. Two-proton correlations may be used to extract this information ; E895 proton interferometry results indicate a proton freezeout volume independent of beam energy . It has been suggested that $`\mathrm{\Lambda }p`$ correlation measurements could offer more sensitive determinations of baryon source sizes than pp interferometry, particularly for large sources. These measurements may also shed light on the $`\mathrm{\Lambda }`$-p interaction, which is not fully understood .

For semi-central ($`b7`$ fm) Au+Au collisions, $`120,000`$ $`\mathrm{\Lambda }`$s were identified in the TPC through their decay topology into $`p+\pi ^{}`$, using a neural network . The $`\mathrm{\Lambda }`$-p correlation function, as a function of the relative momentum in the pair rest frame, is shown in Figure 5. A significant enhancement is seen at low relative momentum, consistent with an attractive $`\mathrm{\Lambda }`$-p final state interaction.

Also shown are expectations for the correlation function for various-sized static sources, using a parameterized final state potential . The shape of the experimental correlation function does not match these expectations for any source size. However, since the effects of the reconstructed $`\mathrm{\Lambda }`$ purity ($`80\%`$) and of $`\mathrm{\Lambda }`$ feed-down from the electromagnetic decay of $`\mathrm{\Sigma }^0`$s (estimated at $`2530\%`$) have not yet been accounted for, it is premature to draw conclusions regarding source size or $`\mathrm{\Lambda }`$-p potential.

## 6 Summary

The E895 Collaboration has performed extensive studies of spectra, interferometry, strangeness, and flow at AGS energies, providing valuable energy systematics useful for identification of “new” physics at SPS and RHIC. Further, we have combined these analyses to extract new physics information. The azimuthal dependence of pion HBT radii reveals the anisotropic transverse shape of the emitting source, and its tilted orientation in the reaction plane, imaging proton and pion flow in coordinate space for the first time. The 6-dimensional $`\pi ^{}`$ phasespace occupation density is consistent with emission from a flowing, thermal source; $`f_\pi ^{}`$ grows with beam energy, approaching the “universal” limit observed at SPS and RHIC energies. First measurements of the $`\mathrm{\Lambda }`$-p correlation function display an enhancement at low relative momentum; more detailed study of this correlation should provide valuable information on the baryon freezeout geometry and the $`\mathrm{\Lambda }`$-p potential.

## Acknowledgements

This work supported by the U.S. Department of Energy under contracts DE-AC03-76SF00098 and DE-AC02-98CH10886, and grants DE-FG02-89ER40531, DE-FG02-88ER40408, DE-FG02-87ER40324 and DE-FG02-87ER40331; by the U.S. National Science Foundation under grants PHY-98-04672, PHY-97-22653, PHY-96-01271, PHY-96-05207, and INT-92-25096; by the University of Auckland Research Committee, NZ/USA Cooperative Science Programme CSP 95/33; and by the National Natural Science Foundation of P.R. China under grant number 19875012.


text:
# Anticipated synchronization in coupled chaotic maps with delays

## 1 Introduction

Time-delay systems have attracted a lot of attention in recent years, in part due to the fact that multistability, i.e. the coexistence of multiple attractors, is a common occurrence when the delays are large –typically, much larger than the response time of the system . Interest in multistability arises because multistable systems play a key role in pattern recognition processes  and memory storage devices. By choosing appropriate initial conditions, prescribed periodic solutions can be stored as oscillatory patterns of a time-delay system . Synchronization of chaotic time-delay systems has also received attention, since it has potential applications to secure communications . Perez and Cerdeira  have shown that, in low-dimensional chaotic systems, a hidden message can be unmasked by the dynamical reconstruction of the chaotic signal using nonlinear dynamical methods. Encrypting a message in the chaotic output of a time-delay system has the advantage that the dynamics is in this case high-dimensional (the dimension increases linearly with the delay ) but, in spite of this fact, synchronization can be achieved by transmitting a single scalar signal. Unfortunately, this method is not as secure as initially expected, since it has been shown that by using a special embedding space, the delay time can be identified, and the message can be successfully unmasked . Coupled oscillators with time delays in the coupling, which represent interactions being transmitted at finite speed, have been extensively studied as well .

Recently, a new effect of delayed coupling was reported by Voss , who showed the existence of an anticipating synchronization regime. In this regime, the slave system becomes synchronized to the chaotic future state of the master system. Anticipation occurs when the coupling is delayed, and results from the interplay of memory effects and relaxation mechanisms. Numerically, this regime was found in coupled semiconductor lasers with optical feedback .

In this paper we investigate the existence and stability of anticipated synchronization in coupled time-delay maps, which have the advantage of allowing for analytical calculations. In the next section, we introduce a system of two time-delay coupled maps in a master-slave configuration, and present analytical results on the stability of anticipated and retarded synchronization for generic maps. In Section 3, we apply the results to a delay map that arises from the discretization of the Ikeda delay-differential equation. Section 4 presents a summary and the conclusions.

## 2 Master-slave coupled delay maps

We consider a one-dimensional map of the form

$$x_{n+1}=bx_n+f(x_{nn_1}).$$
(1)
For $`|b|<1`$, the first term in the right-hand side represents a relaxation mechanism. Under its sole action $`x_n`$ would asymptotically vanish. This relaxation, however, competes with the effect of the nonlinear function $`f(x)`$, which has the form of a time-delayed feedback with delay $`n_1`$.

The (master) map (1) is used to partially drive the evolution of a new (slave) system $`y_n`$. The dynamics of this system is, in principle, the same as for $`x_n`$, except that a part of the nonlinear component is replaced by the evolution of $`x_n`$, namely

$$y_{n+1}=by_n+(1\eta )f(y_{nn_1})+\eta f(x_{nn_2}).$$
(2)
The parameter $`\eta [0,1]`$ measures the strength of the driving action. Note that $`x_n`$ enters the dynamics of $`y_n`$ with a delay $`n_2`$. This form of unidirectional time-delay coupling is not only an extension to maps of the configuration proposed by Voss for continuous-time systems , but also generalizes the master-slave interaction by introducing a parameter ($`\eta `$) that controls its strength. Voss’s configuration corresponds to the extreme value $`\eta =1`$.

The structure of Eqs. (1) and (2) suggests that, if driving is effective enough, $`y_n`$ may become synchronized to the driving signal, in the form

$$y_n=x_{nn_2+n_1}.$$
(3)
If this synchronization is realized, the dynamics of $`y_n`$ coincides with that of $`x_n`$ up to a time shift $`n_1n_2`$. If $`n_1>n_2`$ the slave system anticipates the evolution of the master system, and we have anticipated synchronization. On the other hand, if $`n_1<n_2`$ the synchronization is retarded.

Linear stability of the synchronized state (3) can be analyzed in the usual way. Taking $`y_n=x_{nn_2+n_1}+\delta y_n`$, the evolution of the deviation $`\delta y_n`$ in the limit of infinitesimal deviations is readily derived from (2):

$$\delta y_{n+1}=b\delta y_n+(1\eta )f^{}(x_{nn_2})\delta y_{nn_1},$$
(4)
with $`f^{}df/dx`$. Note that, if $`|b|<1`$ and for $`\eta =1`$, the deviation decreases exponentially with time and synchronization is therefore stable. In view of well-known properties of synchronization in coupled maps without delays, it is possible to advance that, if the dynamics of $`x_n`$ is chaotic, the synchronized state (3) will be stable above a certain threshold $`\eta _c<1`$. On the other hand, if $`x_n`$ is nonchaotic, synchronization should be stable for any $`\eta >0`$.

For $`n_10`$, Eq. (4) can be formally integrated by introducing a linear $`(n_1+1)`$-dimensional map for a variable $`𝐳_n=(z_n^0,z_n^1,\mathrm{},z_n^{n_1})`$, with $`z_n^k=\delta y_{nk}`$. This equivalent map is given by

$$𝐳_{n+1}=M_n𝐳_n,$$
(5)
where the elements of the matrix $`M_n`$ read

$$M_n^{i,j}=\{\begin{array}{cc}b\hfill & \text{if }i=j=0\text{,}\hfill \\ (1\eta )f^{}(x_{nn_2})\hfill & \text{if }i=0,j=n_1\text{,}\hfill \\ 1\hfill & \text{if }i=j+1\text{ }(j=0,\mathrm{},n_11)\text{,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(6)
The deviations given by Eq. (4) can then be obtained from the solution of the equivalent map,

$$𝐳_n=U_n𝐳_0=M_{n1}M_{n2}\mathrm{}M_1M_0𝐳_0.$$
(7)
Thus, synchronization is linearly stable if all the eigenvalues of the evolution matrix $`U_n`$ vanish for $`n\mathrm{}`$. For $`\eta =1`$ the eigenvalues are $`(b^n,0,\mathrm{},0)`$. For $`\eta <1`$, on the other hand, $`U_n`$ and its eigenvalues must be calculated numerically, since they involve the successive states of the variable $`x_n`$, cf. Eq. (6). It is interesting to point out that stability is determined by the asymptotic properties of $`U_n`$ in the limit $`n\mathrm{}`$, where the driving system $`x_n`$ has exhaustively explored the accessible region of phase space. In this limit, thus, the delay $`n_2`$ of Eq. (2) plays no role in the determination the eigenvalues of $`U_n`$. In consequence, $`n_2`$ is irrelevant to the stability of the synchronized state.

## 3 Application to the Ikeda delay map

In this section we present numerical results corresponding to maps (1) and (2), with $`f(x)=a\mathrm{sin}x`$. With this choice of $`f(x)`$, we refer to Eq. (1) as the Ikeda delay map, since it arises from the discretization of the Ikeda delay-differential equation ,

$$\dot{x}=x+\alpha \mathrm{sin}x(t\tau ).$$
(8)
This equation was introduced to describe the dynamics of an optical bistable resonator, where the finiteness of the speed of light is relevant, and therefore require a description that incorporates explicitly the round-trip time of light in an optical cavity. Physically, $`x`$ is the phase lag of the electric field across the resonator, $`\alpha `$ is the laser intensity injected into the system, and $`\tau `$ is the round trip time of the light in the resonator.

Writing $`\dot{x}[x(t+\mathrm{\Delta }t)x(t)]/\mathrm{\Delta }t`$, $`tn\mathrm{\Delta }t`$, and $`\tau n_1\mathrm{\Delta }t`$, the Ikeda delay map is obtained from Eq. (8) with $`b1\mathrm{\Delta }t`$, $`a\alpha \mathrm{\Delta }t`$, and $`x_nx(t)`$. The Ikeda delay map is not to be confused with the discrete Ikeda map , which is obtained from (8) by discretizing time in units of $`\tau `$ in the singular limit where the delay-to-response time ratio diverges.

In the following, we focus the attention on the case $`b=0.9`$. Our numerical calculations are restricted to random initial condition in $`(\pi ,\pi )`$ both for the master and the slave system. In Fig. 1 we illustrate anticipated synchronization, for $`n_1=15`$ and $`n_2=0`$. The slave system (dashed line) anticipates in $`15`$ steps the state of the master system (solid line). After a transient, the difference $`x_{n+n_1n_2}y_n`$ (dotted line) decays to zero. We have verified that this behavior is independent of the value of $`n_2`$.

The critical value $`\eta _c`$ of the coupling strength, above which the synchronized state is stable, has been calculated from the analytical results of the previous section by means of a numerical evaluation of the matrix $`U_n`$ in Eq. (7). Since the delay $`n_2`$ is irrelevant to this calculation, we can take $`n_2=0`$ in Eq. (6). From a given initial condition for the master system, the matrix $`U_n`$ is evaluated for successive values of $`n`$. In order to avoid the calculation of its eigenvalues, which is specially troublesome for large $`n_1`$, $`U_n`$ is multiplied at each step by a randomly generated vector, $`𝐮_n`$, of unitary modulus. If, for a given value of $`\eta `$, the modulus of the product remains above a certain upper threshold $`H`$ for sufficiently large $`n`$, $`|U_n𝐮_n|>H`$ for $`n>n_{\mathrm{max}}`$, synchronization is considered unstable and $`\eta _c>\eta `$. If, on the other hand, $`|U_n𝐮_n|<h`$ for $`n>n_{\mathrm{max}}`$, where $`h`$ is a suitable lower threshold, synchronization is considered stable and $`\eta _c<\eta `$. With this criterion, $`\eta _c`$ can be found by decimation within the interval $`(0,1)`$ up to a certain previously fixed precision. In our calculations, we have taken $`H=10^{10}`$, $`h=10^{10}`$, and $`n_{\mathrm{max}}10^4`$.

Figure 2 shows $`\eta _c`$ as a function of the delay $`n_1`$ and for several values of the parameter $`a`$. Note that for small $`n_1`$ –specifically, for $`n_1=1`$– the critical value $`\eta _c`$ can vanish, which indicates that the master evolution is nonchaotic. Notice, furthermore, that for some values of $`a`$ and $`n_1=1`$ we have plotted more than one value of $`\eta _c`$. This is due to the fact that, in such cases, the master system is multistable –with, typically, two chaotic attractors. The value of $`\eta _c`$ depends thus on the initial conditions for both $`x_n`$ and $`y_n`$. Whereas the behavior is quite irregular for small $`n_1`$, for larger delays the critical coupling becomes practically independent of $`n_1`$. Though it has not been possible to prove this analytically, we conjecture that $`\eta _c`$ approaches a finite limit below unity as $`n_1\mathrm{}`$. The dependence on $`a`$ in the investigated interval is also moderately weak. As expected, since the dynamics of the Ikeda delay map becomes more irregular as $`a`$ grows –i.e. the Lyapunov exponent is higher– $`\eta _c`$ increases accordingly.

The degree of anticipated of retarded synchronization can be quantified by calculating the similarity function $`S_j`$, defined as :

$$S_j^2=\frac{[x_{n+j}y_n]^2}{[x_n^2y_n^2]^{1/2}}.$$
(9)
If $`x_n`$ and $`y_n`$ are independent time series with similar mean value and dispersion, the average of their square difference is of order $`2x_n^22y_n^2`$, and thus $`S_j\sqrt{2}1.4`$. If, on the other hand, there is perfect anticipated or retarded synchronization, the difference $`x_{n+n_1n_2}y_n`$ vanishes, and $`S_{n_1n_2}=0`$. In Fig. 3, the similarity function is shown for the same parameters as in Fig. 1.

Figure 4 shows $`\mathrm{min}(S)`$ –i.e. the minimum value of $`S_{n_1n_2}`$ recorded during sufficiently long realizations of the evolution for a given set of parameters– in the $`(n_1,\eta )`$-plane, with all the other parameters as in Fig. 1. The region where anticipated synchronization occurs, $`\mathrm{min}(S)0`$, is clearly visible for large values of $`\eta `$. The position of its boundary is in good qualitative agreement with the values of $`\eta _c`$ shown in Fig. 2 for $`a=3`$.

## 4 Conclusion

We have extended the study of anticipated synchronization, advanced by Voss for unidirectionally coupled differential equations with time delays , to delayed coupled chaotic maps. While the nature of anticipated synchronization of maps and differential equations is the same, delay discrete-time dynamics admits an analytical treatment which cannot be carried out for continuous-time systems. In fact, ordinary differential equations with finite time delays constitute an infinite-dimensional problem. On the other hand, since times delays in maps must be discrete, the dimensionality of the problem remains finite.

Taking advantage of this situation, we have analytically studied the stability of anticipated and retarded synchronization in a generic master-slave configuration. In the absence of coupling, master and slave dynamics are identical and involve an intrinsic delay $`n_1`$. Coupling consists in the replacement of a part of the slave dynamics by that of the master system, with a delay $`n_2`$. We have shown that the stability of synchronization is independent of $`n_2`$. The structure of the linearized problem, Eqs. (4-7), suggest meanwhile a strong –though not transparent– dependence on the Lyapunov exponent of the master system, as expected. In practice, the linearized problem has to be treated numerically, but it only involves the realization of the master system and the successive application of the linear-evolution operator, being thus a purely algebraic process.

These results have been applied to the Ikeda delay map, which derives from the application of the Euler integration scheme to the Ikeda delay-differential equation. We have calculated the critical coupling intensity above which synchronization is stable, as a function of the delay $`n_1`$ and of a parameter $`a`$ that controls the chaotic dynamics of the map. It is found that, whereas for small values of $`n_1`$ the critical coupling can vary considerably –due to the irregular appearance and disappearance of chaotic and nonchaotic attractors– the dependence for large $`n_1`$ is much smoother. In fact, as $`n_1`$ grows, the critical coupling seems to approach a constant value. The dependence with the dynamical parameter $`a`$ is also moderate in the considered range.

In this work we have focused the attention on exact anticipated synchronization. However, it has been previously shown that approximate anticipated synchronization is possible in coupled differential equations, even in the absence of intrinsic delays . The study of approximate anticipated synchronization in coupled maps constitutes therefore a line for future work. In particular, it would be interesting to investigate in detail the connection between the degree of synchronization, and the irregularity of the dynamics –as measured by the Lyapunov exponents. For delay differential equations, the number of positive Lyapunov exponents and the fractal dimension increase linearly with the delay, while the metric entropy remains roughly constant . We therefore conjecture that the metric entropy might be a good indicator of the possibility of synchronizing with anticipation, and thus predicting, chaotic dynamics.

## Acknowledgements

This work was supported by Proyecto de Desarrollo de Ciencias Básicas (PEDECIBA) and by Comisión Sectorial de Investigación Científica (CSIC), Uruguay.