text:
# Are There Unstable Planetary Systems Around White Dwarfs?

## 1 Introduction

The discovery of $``$ 80 planets, and counting, around solar type stars suggests that successful planet formation is quite common. The wealth of systems so unlike the Solar System leads one to conclude that many aspects of planetary system formation and dynamical evolution have yet to be fully explored. One particularly interesting area is the long term evolution of planetary systems specifically in the presence of post-main-sequence evolution of the central star. Observations of planets around post-main-sequence stars may provide additional information about the formation and evolution of planetary systems around main-sequence stars, and can inform us about the long term future of the solar system.

While planets at distances similar to the outer planets in the solar system will persist through post-main-sequence evolution (Duncan & Lissauer, 1998), it is unlikely that close Jovian companions to such stars survive. As the star evolves it expands, engulfing anything up to $``$ 1 AU (Sackmann et al., 1993; Siess & Livio, 1999a, b). Outwards of an AU, up to $``$ 5 AU any planet’s orbit will decay through tidal transfer of angular momentum and become consumed within the envelope of the star (Rasio et al., 1996). Anything with less mass than a brown dwarf will not survive in the stellar envelope (Livio & Soker, 1984; Soker et al., 1984).

Planets may still be observed in close orbits around white dwarfs if their orbits are significantly changed by some process that occurs after the AGB phase. If the planets become unstable to close approaches with each other, their interaction would result in a planet close to the central star, similar to scenarios proposed for the formation of close Jovian planets around main sequence stars (Rasio & Ford, 1996; Weidenschilling & Marzari, 1996; Lin & Ida, 1997). The results of close encounters between two Jovian planets have been studied in detail, with three possible outcomes: the two planets collide leaving a large planet, one planet is ejected, or both planets remain in a new stable configuration (Ford et al., 2001). For planets starting out with semi-major axes $`>`$ 5 AU, $``$8% of unstable pairs will collide, the remaining $``$92% will not. Of the systems which avoid collision, roughly 40% will settle into a configuration with a planet in a significantly closer and more eccentric orbit than in the initial system. Thus the onset of instability after post-main-sequence mass loss may create white dwarf systems with planets at obital radii scoured clear of the original inner planets during the star’s giant phase.

Duncan & Lissauer (1997), simulating the Uranian satellite system, found that Hill stable systems can become unstable with an increasing mass ratio for satellites orbiting a central massive object. This important work led to one of the few systematic studies of the post-main-sequence evolution of planetary systems dynamically similar to the solar system (Duncan & Lissauer, 1998). In that paper, Duncan & Lissauer found that the time to unstable close approaches for the planets followed a power-law relationship with the ratio of planetary mass to stellar mass, as an increasing fraction of mass was lost from the central star. At the level of individual planetary orbits, resonances may also play an important role in an adiabatically changing system, enhancing stability or creating instability. In this paper we generalize the specific results of Duncan & Lissauer (1998) to a range of system parameters applicable to a wider range of situations, such as those like the multi-planet extrasolar systems recently discovered.

Most likely, extrasolar planetary systems also possess Oort cloud analogues as a direct result of planet formation (Oort, 1950; Weissman, 1999). For comets in the outer Oort cloud orbital time scales are comparable to the more rapid late stages of post-main-sequence evolution, and the mass loss of a star is not adiabatic in the context of AGB evolution. If the mass loss is fairly symmetric, many of these comets survive the evolution and can later provide a steady flow of comets which impact the white dwarf directly, or break up due to tidal strain and populate the inner system with dust, causing photospheric metal contamination (Stern et al., 1990; Alcock, Fristrom, & Siegelman, 1986; Parriott & Alcock, 1998). However, if the planet systems become unstable to close approaches after the AGB phase, the entire system becomes dynamically young and new collisions and encounters can occur between surviving comets and planets. Many scenarios lead to a period of enhanced “late bombardment” as cometary orbits are perturbed and the flow of comets into the inner system is enhanced. In this paper we will investigate whether this mechanism can explain the observed IR excess around the white dwarf G29-38, attributed to a disk of dust whose extent is comparable to the tidal radius of the white dwarf (Graham, et al., 1990; Zuckerman & Becklin, 1987). Graham, et al. (1990) estimated that approximately 1015 g/yr of metal rich material would rain upon the white dwarf, if the infrared emission is due to a dust cloud at about half solar radius, requiring a steady and high rate of replenishment of the putative dust. Our model may resolve the problem of those DAZ white dwarfs which cannot be easily explained by either isolated cometary impacts or ISM accretion (Zuckerman & Reid, 1998).

We will show that the mass lost from a central star is sufficient to destabilize systems of two or more planets in previously stable orbits and cause them to suffer close approaches, producing several observable signatures. In section 2 we will develop the stability of planetary systems against close encounters in the presence of adiabatic mass loss, describe our numerical methods for testing our analytical estimates in section 3, present our results and discuss relevent observational signatures in section 4, and discuss the implications of these findings in section 5.

## 2 Stability for Planetary Systems Under Mass Evolution

### 2.1 Two Planet Systems

The stability of two planets against close approaches depends primarily on the masses of the planets relative to the central star and the separation of the two orbits. This separation is measured as $`\mathrm{\Delta }=(a_2a_1)/a_1`$ where $`a_1`$ and $`a_2`$ are the inner and outer semi-major axes respectively. A critical Hill separation, $`\mathrm{\Delta }_c`$ is then the minimum separation between two planets which ensures a lack of close approaches over all time (Hill, 1886). A full treatment of the Hill stability of two planets in the case of static masses can be found in Gladman (1993). Several approximations can be made that simplify the full treatment, such as equal planetary masses and small eccentricities. The criterion is then given by:

$$\mathrm{\Delta }_c\sqrt{\frac{8}{3}\left(e_1^2+e_2^2\right)+9\mu ^{2/3}}$$
(1)
where $`\mu `$ is the ratio of the planets’ mass to the central star, $`e_1`$, $`e_2`$ are the eccentricities of planet 1 and 2, and $`\mathrm{\Delta }_c`$ is in units of the inner planet’s semi-major axis $`a_1`$.

If either the mass of the planets or the mass of the star changes, the critical Hill radius will change as well. An increase in planet mass or a decrease in stellar mass will cause $`\mu `$ to become larger, increasing the width of the zone in which orbits are unstable to close approaches. During post main-sequence mass loss, the orbits of planets will widen as the central star loses mass. As long as this process is adiabatic, the planets will simply conserve their angular momentum and widen their orbits proportionally to the mass lost: $`a_{new}=a_{old}(M_i/M_f)`$. However, since the orbits widen together by the same factor, $`\mathrm{\Delta }`$ remains the same. Thus, while the critical separation at which the two planets will become unstable widens, their relative separation remains unchanged. Orbits that are initially marginally stable, or close to being unstable, will become unstable to close planet–planet approaches as a consequence of the mass loss from the central star. In the case of planetary mass accretion in a protoplanetary disk, the orbits of the two planets will remain the same while $`\mathrm{\Delta }_c`$ increases, creating the same effect as if the star were losing mass.

The opposite case of stellar mass accretion or planetary mass loss works to make previously unstable regions more stable. However, since close approaches generally happen within a few tens of orbits, objects likely would be cleared out of an unstable region quicker than the region could shrink.

### 2.2 Multiple Planet Systems

We expect that multiple planet systems should be common, e.g. the solar system, PSR 1257+12, and $`\upsilon `$ Andromedae (Wolszczan, 1994; Butler et al., 1999), and thus it is useful to develop an idea of how these systems remain stable. Chambers et al. (1996) found a relation between the separation of a system of planets and the time it would take for the system to suffer a close encounter:

$$\mathrm{log}t=b\delta +c$$
(2)
where $`b`$ and $`c`$ are constants derived through numerical simulations. The symbol $`\delta `$ is related to $`\mathrm{\Delta }`$ but is defined in a slightly different way. Here, $`\delta `$ is the separation between two planets ($`a_{i+1}a_i`$) in units of mutual Hill radii (Ri) defined as such:

$$R_i=\left(\frac{1}{3}\mu \right)^{1/3}\frac{\left(a_{i+1}+a_i\right)}{2}$$
(3)
where $`i`$ can be from 1 to $`N_{planets}1`$ and we assume the planets have equal masses and initially circular orbits. If the parameter $`\delta `$ is the same for each pair of adjacent planets the separations in units of AU will be different. For example, if we took three Jovian mass planets ($`\mu 10^3`$) with $`\delta =6.5`$ and the innermost Jovian at 5.2 AU from the central star, the next two planets would be at 9.4 and 16.7 AU. Compare to actual orbital radii of 9.6 and 19.1 AU for Saturn and Uranus respectively, in our presumably stable for several billion years solar system. We add the obvious caveat that Saturn and Uranus are significantly less massive than Jupiter and have correspondingly weaker mutual interactions.

Adiabatic mass evolution will have the effect of shortening the time it takes for orbits to suffer close approaches. The knowledge of this has long been used to speed up numerical calculations (Duncan & Lissauer, 1997, and references therein). However, this fact also leads to the hypothesis that planetary systems on the edge of stability for 1010 yr will be affected by mass loss. In general the new time to close aproaches for an initial $`\delta `$ with a change in mass is given by:

$$\mathrm{log}\frac{t_f}{t_i}=(b^{}b)\delta +(c^{}c).$$
(4)
We would expect $`c`$ to have little or no change with a change of mass since it represents the timescale for two planets at $`\delta `$0 to suffer a close approach. Mass loss will increase the mutual Hill radii of the planets which in turn will change $`b`$ to a new value we will define as $`b^{}`$:

$$b^{}=\left(\frac{\mu _f}{\mu _i}\right)^{\frac{1}{3}}b$$
(5)
where $`\mu _f`$ and $`\mu _i`$ are the final and initial mass ratios respectively. Such behavior suggests that bodies that are stable over the lifetime of a planetary system will become unstable over a timescale several orders of magnitude smaller than their original timescale for instability, when the central star becomes a white dwarf, assuming the relation of Equation 2 holds for large $`\delta `$. It has been found that for the case of three planets with $`\mu `$=10-7, the parameters $`b1.176`$ and $`c1.663`$ (Chambers et al., 1996). If the three planets are each separated from their neighbor by $`\delta `$=6, they will experience close encounters after $``$105 orbits of the inner planet. For comparison, three planets with the same mass ratio and separated by a $`\delta =8`$ will experience close encounters after $`6\times 10^7`$ yr. Assuming the central star loses half of its mass, the timescale to close encounters will shorten by an order of magnitude for the first case and two orders of magnitude for the second.

### 2.3 When is Mass Evolution Adiabatic?

The question of whether mass evolution is adiabatic needs to be addressed. In the case of mass loss for solar mass stars, roughly half the central star’s mass will be lost on the order of $`10^8`$ yr. A majority of the mass is lost at the tip of the AGB branch during a period of $``$106 yr. Even the quickest rate of mass loss is much longer than one orbital period of a planet inwards of 100 AU, the general region where planets are believed to have formed. Stars heavier than a solar mass probably have superwinds which will cause significant mass loss on the order of a few hundred or thousand years (Vassiliadis & Wood, 1993; Schröder et al., 1999). Whether this is important or not will be the subject of further study. Objects very far away from the central star, such as Oort cloud object analogues, have orbital timescales comparable to the mass loss timescale and will also not follow the adiabatic case. It should be noted that for Kuiper and Oort cloud distances that the mass loss by the star would become adiabatic if the asymptotic wind velocity were orders of magnitude smaller than the escape velocity at the surface of the star since the crossing time of the wind would then be larger than the orbital timescale of the comets.

Mass gain of stars and planets is much slower than the orbital timescale of a planet. Accretion rates for protostars are on the order $`10^6`$ M yr-1 (Shu et al., 1987). The formation of giant planets through runaway gas accretion takes $``$107 yr, the rough lifetime of gaseous protoplanetary disks (Pollack et al., 1996). If some giant planets are formed more quickly by more efficient runaway accretion, gravitational collapse (Boss, 2000), or seeding through the formation of other planets (Armitage & Hansen, 1999), they would not be described by the adiabatic case.

## 3 Numerical Methods

In order to test the hypothesis that adiabatic mass evolution should change the stability of planetary systems, we ran several numerical simulations of two planet and multi-planet systems in circular orbits around a central star losing mass. The equations of motion were integrated using a Bulirsch-Stoer routine (Stoer & Bulirsch, 1980; Press et al., 1992). Since the case of mass loss of the central star and mass gain of the planet is the same, mass loss can be modeled in two ways. Either the star’s mass can be decreased, or the planets’ masses can be increased. If the planets’ masses are increased the time coordinate must be scaled to reflect the fact that the orbits are widening. To keep our investigations scalable, we chose the units of time to be orbits of the inner planet. We chose to increase the mass of the planets over a period of 1000 orbits. In the absence of mass evolution, energy and angular momentum were conserved to better than 1 part in 106 for 105 orbits. Since changng mass makes this a non-conservative system, energy and angular momentum could not be used as a test of accuracy. However, since the simulations were integrated until a close approach and then terminated, any error is similar to the case of no mass evolution. Several simulations without mass evolution were run with stable results. A close approach was defined by an encounter separated by a radius of $`<`$ 2$`\mu ^{2/5}`$ (Gladman, 1993, and references therein). This radius was chosen because at separations smaller than this the planet-planet system is dominant and the star becomes a perturbation. Other authors have chosen different criteria (Chambers et al., 1996), but the results are insensitive to the exact choice.

In the two planet case, we started simulations at the critical separation predicted by Equation 1 assuming no mass loss, and increased the separation between the two planets at regular intervals in $`\mathrm{\Delta }`$. We integrated the equations of motion until a close approach or for 105 orbits. We increased $`\mathrm{\Delta }`$ until it was 25% greater than what would be predicted in the presence of mass loss. These simulations were run an order of magnitude longer than Gladman (1993), and in the no mass loss case were consistent with what they found. The two planets were initially started with true anomalies separated by 180. Our separations are then lower limits for the critical separation and thus truly reflect the minimum possible separation between orbits that remain stable. For multiple planets, $`\delta `$ was started at 2.2 and raised until several consecutive separations did not experience close encounters for 107 orbits. Here random phases in the orbits where chosen with the restriction that adjacent orbits were separated by at least 40. Three separate runs with different random initial phases were performed to improve the statistics for each mass, as there is significant scatter in the actual time to a close approach for each separation.

## 4 Results

### 4.1 Two Planets

We looked at a wide range of planetary masses for a solar mass star, from a sub terrestrial-sized planet ($`\mu =10^7`$) to a Jovian planet ($`\mu =10^3`$). Figure 1 shows the border for onset of instability in two planet systems after mass loss. The dashed line represents the initial critical Hill radius for no mass loss. The solid line which goes through the points is the critical Hill radius for $`\mu `$ equal to twice that of the initial system, corresponding to the planets doubling in mass or the central star losing half of its mass. Several of the higher mass points are greater than that predicted by the solid curve, an indication of higher order $`\mu `$ terms becoming important. It should be noted that these results are general to any combination of planet and stellar mass that have these ratios.

In a few cases, separations predicted to become unstable after mass loss by the Hill criterion were stable for the length of our simulations. Particularly in the $`\mu `$=10-3 case, there was a large region in which the two planets suffered no close encounters (See Figure 2). These orbits corresponded to a range of $`\mathrm{\Delta }`$ from .32 to .37, which were predicted to be unstable under mass evolution from the simple scaling of the equation for $`\mathrm{\Delta }`$. It is interesting to note that all of these orbits are close to the 3:2 resonance (See Fig. 3). For the same reason that the Hill radius will not change, these orbits will retain the ratio of their periods. The reason for the stability around the 3:2 resonance may be due to those separations being near but not in a region of resonance overlap (Wisdom, 1980; Murray & Holman, 2001), clearly this conjecture needs to be confirmed but that goes beyond the scope of this particular paper.

### 4.2 Multiple Planets

Figures 4-6 show the results for three different runs, looking at three planet systems in circular orbits. We looked at the mass ratios $`\mu `$=10-7, 10-5, and 10-3. The results are compared to simulations without mass loss, and the difference between the two is quite noticeable for the whole range of mass. It is important to note that for separations whose time to close approach is comparable to the mass loss timescale show little change in behavior between the two cases. This is because the change in the time to close approach is smaller than the scatter in the simulations. Least squares fitting of the static and mass loss cases were performed to get the coefficients $`b`$, $`c`$, and $`b^{}`$. To test our assumption of $`c`$ not changing under mass evolution, we also measured $`c^{}`$, the intercept for the mass loss case. Planets with initial separations in $`\delta `$ that were less than $`2\sqrt{3}`$, the two planet stability criterion in units of $`R_i`$, were discarded. For the mass loss case, points where the timescale of close approaches were comparable to the mass loss timescale were also discarded. Once the coefficients were determined they were compared to what was predicted from Equation 2. Similarly, $`b`$ and $`c`$ from the $`\mu `$=10-7 case without mass loss were compared with the results of Chambers et al. (1996). Table 1 shows that within the uncertainties, $`c`$ indeed does not change with mass evolution and the slopes are consistent with predictions. Additionally, our results for the static case with $`\mu `$=10-7 are consistent with Chambers et al. (1996)’s values for $`b`$ and $`c`$.

As mass increases, the presence of strong resonances becomes more important. This is due to our choice of equal separations and equal masses, many of these resonances would disappear with small variations in mass, eccentricity, and inclination (Chambers et al., 1996), aspects that will be tested with further study. The presence of resonances is most easily seen in Figure 6 where $`\mu `$=10-3. In the range of $`\delta `$=4.4 to 5.2, the points greatly depart from the predicted curve. the spike at $`\delta `$=5.2 corresponds to the first and second, as well as the second and third planets being in 2:1 resonances. This particular example shows that the basic dynamics of a system undergoing adiabatic mass evolution favor stability near strong resonances. Such a process potentially could augment the current ideas about how resonant extrasolar planets such as those around GJ876 formed (Snellgrove et al., 2001; Armitage et al., 2001; Murray et al., 2001; Rivera & Lissauer, 2000).

### 4.3 Observational Implications

These simulations have several observational implications which can be broadly separated into two categories–the character and the signature of planetary systems around white dwarfs.

Surviving planets that are marginally stable will suffer close approaches soon after the star evolves into a white dwarf, or possibly as early as the AGB phase. There are three possible end states for planets that suffer close approaches, ejection, collision, and a settling into a different and more stable configuration for all planets. The case of two planets has been studied carefully, and for two Jovian mass planets with one planet starting at $``$5 AU, the probability of collision is roughly 8%, ejection 35%, and rearrangement 57% (Ford et al., 2001). We naively assume that these results hold similarly for multiple planets as well, since collisions have been shown to hold for multiple planet systems (Lin & Ida, 1997), while ejections and rearrangements should have similar probability (certainly to within a factor of 2 or so). Ejections will leave planets that are closer to the white dwarf, while often a rearrangement will leave one or two planets with larger semi-major axes (up to $``$ 103 times greater) and one with a smaller semi-major axis (as close as .1 times smaller). Collisions are potentially more exciting because as the two planets merge they essentially restart their cooling clock and as such will be anomalously luminous by 2 orders of magnitude for 108 yr (Burrows et al., 1997).

To estimate how many white dwarfs might have planets that collided ($`F_c`$) we can take the fraction of white dwarfs that have marginally stable planets and multiply them by the fraction of marginally stable planets that have collisions:

$$F_c=f_{pl}f_{ms}f_c$$
(6)
where $`f_{pl}`$ is the fraction of white dwarfs with planets, $`f_{ms}`$ is the fraction of marginally stable planet systems, and $`f_c`$ is the fraction of marginally stable systems that suffer a collision. We can estimate the number of Jovian sized planets around white dwarfs by looking at the number of young stars that still have significant disks after 1 Myr, the approximate time to form a Jovian planet. This has been found to be about 50% of young stars in nearby clusters (Haisch, Lada, & Lada, 2001). Several numerical simulations (Barnes & Quinn, 2001; Laughlin & Adams, 1999; Quinlan, 1992; Rivera & Lissauer, 2000; Stepinski et al., 2000, for example,) point to a high frequency of marginally stable systems around stars as well as the discovery of the marginally stable planetary systems around GJ876 and HD 82943 (Marcy et al., 2001). However, factors such as multiple planets with widely different mass ratios could greatly change the effects of stability and need to be studied further. We estimate this fraction to be about 50% as well, although a large uncertainty is associated with this estimate. Taking the results above, we estimate then that $`2(f_{ms}/0.5)`$% of young white dwarfs should have the product of a recent planet-planet collision in orbit. Thus we predict that observations of young ($`\tau 10^9`$ yr) white dwarfs should reveal $`2`$% have overluminous planet mass companions, some in orbits with semi-major axis smaller than the minimum (5 AU) expected to survive the AGB phase. These planets would be detectable through their significant IR excess, and should be distinguishable from brown dwarf companions.

A natural byproduct of the formation of Jovian planets is the existence of a large cloud of comets at large heliocentric distances (Oort, 1950; Weissman, 1999). The survival of such a cloud through post main-sequence evolution has been closely studied in the context of accounting for observed water emission in AGB stars and an explanation for metals in DA white dwarfs (Stern et al., 1990; Alcock, Fristrom, & Siegelman, 1986; Parriott & Alcock, 1998). The general result to date is that comets at semi-major axes greater than a few hundred AU survive the AGB phase. Then massive comets are predicted to strike the central white dwarf at a rate of $`10^4\mathrm{yr}^1`$, depositing fresh metals in the white dwarf photosphere. Such a cometary influx can account for the DAZ phenomenon, but has difficulty explaining some of the strongest metal line systems, because of the short predicted residence times of metals in the photosphere of these white dwarfs. An alternate explanation for the origin of metals in white dwarfs is ISM accretion, where a steady drizzle of metal rich dust is spherically accreted from the ambient ISM. Both scenarios have difficulty explaining the frequency of DAZ white dwarfs, and accounting for the systems with strongest metal lines and shortest residence times, due to the high mass accretion rate required to sustain those systems.

Recent observations of the DAZ phenomenon do not seem to be consistent with either scenario (Zuckerman & Reid, 1998). In one DAZ, G238-44, the diffusion time for metals is 3 days, which means that neither previous scenario can explain the high observed metal abundances nor the stability of the metal lines (Holberg, Barstow, & Green, 1997). Another white dwarf, G29-38, has a high abundance as well as an infrared excess, possibly from a dust disk at small orbital radii (Zuckerman & Becklin, 1987; Koester, Provencal, & Shipmann, 1997).

The evolution of a planetary system after post-main sequence mass loss coupled with the presence of an Oort cloud may provide an alternative explanation for the DAZ phenomenon and in particular the peculiarities of G29-38 and G238-44.

#### 4.3.1 Cometary dynamics

The mass loss during post-main-sequence is near impulsive for Oort cloud comets. Previous work has shown that a significant fraction of any Oort cloud like objects will survive the mass loss phase, even in the presence of mildly asymmetric mass loss (Alcock, Fristrom, & Siegelman, 1986; Parriott & Alcock, 1998). The immediate result of the mass loss phase is to leave the remaining bound objects on orbits biased towards high eccentricity, but with similar initial periastrons. Orbital time scales are on the order of $`10^6`$ yr.

The number and typical size of Oort cloud objects is poorly constrained, canonical estimates scale to 1 km sized comets, mostly composed of low density ices and silicates, with masses of $`10^{16}`$ gm each, with order $`10^{12}`$ objects per star. Clearly there is a range of masses, and it is possible the true numbers and masses of Oort cloud comets vary by several orders of magnitude from star to star. Dynamical effects also lead to a secular change in the amount of mass in any given Oort cloud.

External perturbations ensure a statistically steady flux of comets from the outer Oort cloud into the inner system, where interactions with Jovian planets lead to tidal disruption of comets (and direct collisions), scattering onto tightly bound orbits restricted to the inner system, ejection from the system, and injection into central star encountering orbits. For the solar system, the flux of comets into orbits leading to collision with the Sun is of the order $`10^2`$ per year, of these a significant fraction undergo breakup before colliding with the Sun, with individual fragments colliding with the Sun over many orbital periods (Kreutz sungrazers). SOHO detects $`10^2`$ such objects per year in the solar system, or one every 3 days or so on average. A single 1 km comet can fragment into $`10^4`$ fragments with sizes of order 50 m, consistent with those observed by SOHO, and consistent with the collision rates estimated both for the parent comets and the fragments. Each fragment then deposits about $`10^{12}`$ gm into the Solar photosphere. Note that if the typical comet were 20 km rather than 1 km, the deposition rate would be about $`10^{16}`$ gm every three days.

A white dwarf has a radius about 0.01 of the solar radius. Due to gravitational focusing, the cross-section for collision for comets scattered into random orbits in the inner system is linear in radius, so the collision rate expected for a white dwarf with a solar-like Oort cloud is $`10^4\mathrm{yr}^1`$. However, the perturbation of the outer orbits due to AGB mass loss, combined with the expansion of the outer planet orbits will drastically change this rate, leading to a new, late “heavy bombardment” phase with significantly higher rates of comet influx into the inner system. If one of the outer (Jovian mass) planets is scattered into a large ($`a_{fin}a_{in}`$) eccentric orbit after the onset of instability, as we expect to happen in about 2/3 of the cases, then there will be strong periodic perturbations to the outer Kuiper belt and inner Oort cloud. About 10% of those systems will lead to the outermost bound planet being placed on very wide ($`a_{fin}>10^3a_{in}`$) highly eccentric orbits, with orbital time scales comparable to the cometary orbital timescales. Perturbations on the Oort cloud from these planets lead to a persistent high flux of comets to the inner system, until the Oort cloud is depleted of comets.

The net effect of the dynamical rearrangement of the post-main-sequence planetary system is a greatly enhanced rate of cometary influx into the inner system, starting $`10^7`$-$`10^8`$ years after the mass loss phase, tapering off gradually with time on time scales of $`10^8`$$`10^9`$ years, leading to enhanced metal deposition to the white dwarf photosphere, and increased dust formation in the inner system for some white dwarfs, depending on the final configuration of the outer planets.

Several processes affect the comet bombardment rate:

* a fraction of the previously stably orbiting outer Kuiper belt objects, that survived the AGB phase, are injected into the inner system by newly established dynamical resonances with the outer planets over $`10^8`$ years;
* planets ejected to the outer Oort cloud by planet-planet perturbations will randomise the orbits of a small ($`4(m/M)^2`$) fraction of the Oort cloud comets, some of these will enter the inner system providing an enhanced flux of the normal Oort comet infall over $`10^9`$ years;
* surviving inner planets, scattered to the smaller orbital radius, will trap the comets injected into the inner system, providing both direct tidal disruption at a few AU, and providing a much higher influx of comets to very small radii where they are tidally disrupted by the white dwarf (or in rare cases collide directly);
* dust from tidally disrupted cometary debris will be driven to the white dwarf surface by PR drag, while larger debris will be dragged in through the Yarkovsky affect, both on a timescale shorter than the WD cooling time.

We expect the Kuiper belt to be severely depleted by the post-main-sequence phase (Stern et al., 1990; Meinick et al., 2001), however, a substantial population of volatile depleted rocky bodies may survive the AGB phase in the outer belt. A substantial fraction of these burnt out comets will become vulnerable to resonant perturbations by the surviving outer planets, now in new, wider orbits. The outer belt objects have orbital periods ($`10^4`$-$`10^5`$ years) that are becoming comparable to the most shortest AGB mass loss timescales, and therefore will not generally expand adiabatically in proportion to the expansion of the planetary orbits. The solar Kuiper belt is inferred to have $`10^5`$ objects with size above 100 km, assuming a mass function with approximately equal mass per decade of mass, characteristic of such populations, we infer a population of $`10^{11}`$ Kuiper belt objects with size of about 1 km, at an orbital radius of order $`10^3`$ AU. Order 1% of those will be vulnerable to the new dynamical resonances after the AGB phase, allowing for evaporative destruction and ejection, we estimate $`10^8`$ Kuiper belt objects will enter the inner system in the $`10^7`$-$`10^8`$ years after the AGB phase. The rate will peak at $`10^8`$ years and then decline as the reservoir of cometary bodies in orbits vulnerable to the new planetary resonances declines.

If there are multiple surviving Jovian planets, then the post-AGB planet-planet interactions will typically leave the inner planets on eccentric orbits, leading to broader resonances and a larger fraction of perturbed Kuiper belt objects. We expect in $`2/3`$ of the cases where there were multiple, initially marginally stable Jovian planets in the outer system, for the final configuration to have an eccentric outer planet, and a more tightly bound inner planet.

Some of the comets injected into the inner system will be tidally disrupted by the surviving Jovian planets, some will be ejected, and some will be injected into the inner system and be tidally disrupted by the white dwarf (and about 1% of those will directly impact the white dwarf). Dynamical time scales in the inner system are $`10^2`$ years, and the probability of ejection or disruption per crossing time is of the order of $`10^2`$ per crossing time, assuming there is an inner planet, scattered inward of 5 AU, matching the outer planet scattered to wider orbital radius. So at any one time $`>10^3`$ Kuiper belt objects are in the inner system. The rate for tidal disruption by the surviving innermost Jovian planets is $`10^6\mathrm{yr}^1`$ per comet. Tidal disruption rates due to close approaches to the white dwarf may be as high as $`>10^4\mathrm{yr}^1`$ per comet, the rates are uncertain because of the possibility of non-gravitational processes breaking up the comet and deflecting debris. With each comet massing about $`>10^{16}`$ gm, by hypothesis, we get a flux of disrupted cometary material, from the Kuiper belt remnant, assuming solar system like populations, of $`10^{14}`$$`10^{16}`$ gm per year. This is volatile depleted material, by hypothesis. Given our assumed mass function, disruption of rarer more massive comets can sustain mass accretion rates an order of magnitude higher still for time scales comparable to the inner system dynamical time scales, in a small fraction of systems.

We can now compare our mechanism to the accretion rates needed to explain the constant, detected metal lines in G29-38 and G238-44, two DAZ white dwarfs with the highest measured abundances of Ca. An accretion rate has already been quoted in the literature for G238-44, where $`3\times 10^{17}\mathrm{g}\mathrm{yr}^1`$ of solar abundance ISM would need to be accreted (Holberg, Barstow, & Green, 1997). In a volatile depleted case, only metals would be present converting to only $`6\times 10^{15}\mathrm{g}\mathrm{yr}^1`$ for cometary material. G29-38 also has a roughly estimated value of $`1\times 10^{19}\mathrm{g}\mathrm{yr}^1`$ corresponding to $`2\times 10^{17}\mathrm{g}\mathrm{yr}^1`$ in the volatile depleted case (Koester, Provencal, & Shipmann, 1997).

Both of these estimates were based on calculations made by Dupuis, Fontaine, & Wesemael (1993); Dupuis, Fontaine, Pelletier, & Wesemael (1993, 1992), who uses the ML3 version of mixing length theory. In fact, there are several other methods that can be used to model the convective layer of white dwarfs, including using other efficiencies of the ML theory and the CGM model of convection (Althaus & Benvenuto, 1998). The calculations differ by up to four orders of magnitude on the mass fraction $`q`$ of the convection layer’s base in white dwarfs with Teff similar to G29-38. Taking values of $`q`$ for the base of the convection layer from Figures 4 and 6 in Althaus & Benvenuto (1998) and getting values for the diffusion timescale from tables 5 and 6 in (Paquette, et al., 1986) one can estimate what steady state accretion rate G29-38 requires for the different models. The smallest rate came from ML1 theory and the largest from ML3 theory, with CGM having an intermediate value, giving a range of $`2\times 10^{13}\mathrm{g}\mathrm{yr}^1`$ to $`4.4\times 10^{17}\mathrm{g}\mathrm{yr}^1`$. We favor the CGM value of $`10^{15}\mathrm{g}\mathrm{yr}^1`$, which is consistent with the estimate based on observations conducted by Graham, et al. (1990). The rate for G238-44 may be more robust due to the fact that convective models converge for hotter white dwarfs.

Both rates are consistent with our scenario if either white dwarf has two Jovian mass planets, one in a $`>10`$ AU, eccentric orbit, and another in a $`<5`$ AU orbit, after rearrangement. Alternatively, their progenitors had an order of magnitude richer Kuiper belt population than inferred for the Solar system. With a post-AGB age of $`6\times 10^8`$ years, and a mass of $`0.7\mathrm{M}_{}`$, the original main sequence star of G29-38 was most likely more massive than solar and a more massive planetary and cometary system is not implausible. G238-44 is almost 108 years old and would represent an object close to the peak of predicted cometary activity.

Our scenario may provide a consistent picture for the presence of DAZ white dwarfs and their anomalous properties (Zuckerman & Reid, 1998). We don’t expect all white dwarfs to have metal lines. Only about 2/3 of those which possessed marginally stable planetary systems containing two or more Jovians at orbital radii greater than $`5`$ AU will be able to generate significant late cometary bombardment from the outer Kuiper belt and inner Oort cloud. Following a similar estimate as in Equation 6, we predict about 14% of white dwarfs will be DAZs. Further, the rate will peak after $`10^8`$ years, after the planet-planet perturbations have had time to act, and then decline as the resevoir of perturbable comets is depeleted. The convective layer of the white dwarf will also increase by several orders of magnitude over time which would create a sharp drop of high abundance DAZs with decreasing Teff. The drop would be greatest between 12000K and 10000K where the convective layer has its steepest increase (see fig. 4 of Althaus & Benvenuto (1998)). Zuckerman & Reid (1998), in conducting their survey of DAZ white dwarfs, estimated that $``$20% of white dwarfs were DAZ and that metal abundance dropped with Teff sharply between 12000K and 8000K.

We expect DAZ white dwarfs to have potentially detectable (generally) multiple outer Jovian planets, whose orbits will show dynamical signatures of past planet-planet interaction, namely an outer eccentric planet, and an inner planet inside the radius scoured clean by the AGB phase.

## 5 Discussion

Using the above results we can compare the greatest fractional change in stability for two Jovian planets around different stars $`>`$ 1 M that produce white dwarfs. We took the initial-final mass relation of (Weidemann, 2000) and calculated $`\mathrm{\Delta }_c`$ without mass loss and with mass loss ($`\mathrm{\Delta }_c^{}`$). As can be seen in Figure 7, the higher the initial mass star, the greater the fractional change. This is expected, since higher mass stars lose more mass to become white dwarfs. The best candidates for unstable planetary systems would be higher mass white dwarfs, if planet formation is equally efficient for the mass range considered here. The scheme conjectured in this paper, provides a method for identifying and observing the remnant planetary systems of intermediate mass stars, which might otherwise be hard to observe during their main sequence life time.

One can predict the change in the critical separation where two planets will remain stable based on the change in $`\mu `$ over time, simply by differentiating $`\mathrm{\Delta }_c`$ with time:

$$\frac{d\mathrm{\Delta }_c}{dt}=\mu ^{2/3}\frac{d\mu }{dt}.$$
(7)
In the general case $`d\mu /dt`$ depends on two factors, the change in mass of the central star and the change in mass of the planets, given by

$$\frac{d\mu }{dt}=\mu \left(\frac{d\mathrm{ln}M_{pl}}{dt}\frac{d\mathrm{ln}M_{}}{dt}\right).$$
(8)
For the critical separation to widen, $`\mu `$ must increase with time. Putting the two equations together gives the rate of change in $`\mathrm{\Delta }_c`$:

$$\frac{d\mathrm{\Delta }_c}{dt}=\mu ^{\frac{1}{3}}\left(\frac{d\mathrm{ln}M_{pl}}{dt}\frac{d\mathrm{ln}M_{}}{dt}\right).$$
(9)
The results of our multi-planet simulations are scalable to many situations, but for planet systems surviving around white dwarfs we are interested in timescales of $``$1010 yr for solar type stars to $``$108 yr for higher mass stars. The highest $`\delta `$ we studied for $`\mu `$=10-3 was roughly 5.2, which by Equation 2 corresponds to a timescale to close approaches of 107 orbits of the inner planet. After the central star loses half of its mass, the timescale shortens to $``$2000 orbits. For a planetary system with $`\delta `$=5.2 to be stable over the main sequence lifetime of the star, the minimum semimajor axis of the innermost planet for a higher mass star (say, 4M) would be 8.2 AU and 100 AU for a solar-type star. Longer integrations need to be performed to investigate the behavior of systems with larger values of $`\delta `$. We expect our results for time scales for onset of instability to scale to larger $`\delta `$, the initial computational effort we made here limited the exploration of slowly evolving systems with large $`\delta `$ in exchange for a broader exploration of the other initial condition parameters. It will also be instructive to model systems with unequal mass planets, to explore the probability of ejection and hierarchical rearrangement as a function of planetary mass ratio.

By performing numerical integrations of two planet and multiple planet systems we have shown that the stability of a system changes with mass evolution. In the specific case of mass loss as the central star of a planetary system becomes a white dwarf, we have found that previously marginally stable orbits can become unstable fairly rapidly after the mass loss process. Coupled with our knowledge of the survival of material exterior to outer planets such as Kuiper Belt and Oort Cloud analogues (Stern et al., 1990), a picture of the evolution of circumstellar material over the latter stages of a star’s lifetime becomes clear.

As a star reaches the RGB and AGB phases, inner planets are engulfed by both the expanding envelope of the star and through tidal dissipation. The surviving planets move slowly outwards, conserving their angular momentum as the star loses its mass over several orbital periods of the planets. The planets may sculpt the resulting wind of the giant star (Soker, 2001) and if they are on the very edge of stability undergo chaotic epsiodes during the AGB phase, creating some of the more exotic morphologies in the resulting planetary nebula. When the star becomes a white dwarf, two planet systems that are marginally stable will become unstable and suffer close approaches, while for three or more planets the timescale to close approaches shortens by orders of magnitude. There are three possible outcomes once the planets start suffering close approaches: the planets collide, one planet is ejected, or the two planets remain but are in highly eccentric orbits (Ford et al., 2001). One major open question is how many marginally stable systems there are, but there are indications that many, if not most, general planetary systems should be close to instability (Barnes & Quinn, 2001; Laughlin & Adams, 1999; Quinlan, 1992; Rivera & Lissauer, 2000; Stepinski et al., 2000) Rocky material in the inner edge of the Kuiper Belt, which is defined by the last stable orbits with respect to the planetary system, will follow the same fate as marginally stable planets, suffering close approaches with the planetary system and becoming scattered into the inner system, which increases the rate of close encounters with the planets or the central white dwarf. The surviving material at outer Kuiper Belt and Oort Cloud distances will have orbital periods comparable to the timescale of the central star’s mass loss. These objects have their eccentricity pumped up by the effectively instantaneous change in the central star’s mass, and then through interactions with planets create a new dust disk around the white dwarf, and contaminate the white dwarf photosphere to an observable extent.

The sensitivity of stability to changes in mass has implications for planet formation as well. Further research on the migration of Hill stable regions while the planet/star mass ratio evolves may illuminate further the general issue of how Jovian planets in the process of formation become unstable to close encounters and gross changes in orbital parameters (Ford et al., 2001). One possibility is that the mass accretion of the planets occurs at a rate fast enough that d$`\mu /dt>0`$. Other factors would need to be considered, in particular the interplay between the onset of rapid mass accretion by the planet, and the accretion rate from the protoplanetary disk onto the central protostar. Gas drag and stellar mass accretion could work to stabilize orbits if the planets are embedded in a circumstellar disk, while orbital migration would change the relative separations of proto-planets. Since the stability of multi-planet systems is also sensitive to changes in mass ratio, this could help solve problems of isolation for planetary embryos and speed up the timescale for the production of giant planet cores.

The dependence of stability on both the mass of the planet and the mass of the central star suggests that stars of different masses may be more efficient at producing a certain size planet. This is exemplified by the fact that $`\mu `$ for a Jovian planet can change by an order of magnitude in either direction over the mass range of stars that might have planetary companions. For larger mass stars, planets can be more tightly spaced and still be mutually dynamically stable, which suggests that when planets are forming it is easier for them to become dynamically isolated in disks around more massive protostars. For lower mass stars, there is a wider annulus in which material is unstable to planetary gravitational perturbations, and so forming planets would have a larger reservoir of material to draw from. Other factors that need further research would play into this result as well and may dominate over this scenario, such as a star’s temperature and radiation pressure. However, such effects will tend to reinforce the conclusion that less massive stars should be more efficient at creating more massive planets while higher mass stars will produce more, lighter planets if they are capable of forming planets at all. This prediction will be testable as many space and ground based programs are devoting a great deal of effort to look for planetary companions to stars.

We wish to thank Hans Zinnecker, Bill Cochran, and Brad Hansen for fruitful discussions. We would also like to extend our special thanks to Dr. James Liebert for helpful suggestions regarding the DAZ/planet connection. J.D. acknowledges funding by a NASA GSRP fellowship under grant NGT5-119. S.S. acknowledges funding under HST grant GO-8267.


format: text

text:
# Reconstruction of motional states of neutral atoms via MaxEnt principle

## I Introduction

A reconstruction of states of quantum systems from experimental data represents an important tool for a verification of predictions of quantum theory. It also allows us to check the fidelity of quantum state preparation as well to study the fidelity of processing of information encoded in states of quantum systems. Complementary, discrepancies between estimated (reconstructed) states based on the measured data and theoretical predictions can serve as an indicator of various noise sources which occur during the quantum information processing or in the measurement of quantum states. Without a priori assumptions about the character of physical processes and properties of reconstruction schemes the reconstruction cannot distinguish between imperfections related to an incoherent quantum state processing and non-ideal measurements . Determination of limits for coherent control of quantum degrees of freedom or identification of sources of decoherence are essential for systems which are considered for quantum computing and information processing .

In atomic optics a highly coherent control of motional degrees of freedom has been achieved for trapped ions  and recently also for neutral atoms . Cold atoms can be cooled into specific quantum states within micro-wells of an optical lattice which is induced by laser beams. Cold neutral atoms in optical lattices represent a promising system for quantum information processing. To verify a degree (fidelity) of coherent control over motional degrees of freedom of neutral atoms a reconstruction of their motional quantum states from measured data has to be considered. We develop a reconstruction procedure based on the Jaynes principle of Maximum Entropy (MaxEnt)  to achieve this goal. This scheme allows us to perform a state reconstruction from the experimental data obtained at the ENS in Paris.

In Section II we present a brief description of the reconstruction procedure of a density operator of a quantum system based on the Jaynes principle of Maximum Entropy. In Section III we utilize the MaxEnt principle for development of the reconstruction scheme of motional states of atoms and we perform numerical tests of our approach in Section III. Our reconstruction scheme is applied to the experimental data in Section IV.

## II MaxEnt principle and reconstruction of density operators

Let us assume a set of observables $`\widehat{G}_\nu `$ $`(\nu =1,\mathrm{},n)`$ associated with the quantum system under consideration. This system is prepared in an unknown state $`\widehat{\rho }`$. Let us assume that from a measurement performed over the system meanvalues $`\overline{G}_\nu `$ of the observables $`\widehat{G}_\nu `$ are found. The task is to determine (estimate) the unknown state of the quantum system based on the results of the measurement. Providing the set of the observables $`\widehat{G}_\nu `$ is not equal to the quorum (i.e., the complete set of system observables ), then the measured meanvalues do not determine the state uniquely. Specifically, there is a large number of density operators which fulfill the conditions

$`\mathrm{Tr}\widehat{\rho }_{\{\widehat{G}\}}`$ $`=`$ $`1,`$ (1)
$`\mathrm{Tr}(\widehat{\rho }_{\{\widehat{G}\}}\widehat{G}_\nu )`$ $`=`$ $`\overline{G}_\nu ,\nu =1,2,\mathrm{},n,`$ (2)
that is, the the normalization condition and the constraints imposed by the results of the measurement. To estimate the unknown density operator in the most reliable way we utilize the Jaynes principle of maximum entropy (MaxEnt principle)  according to which among those operators which fulfill the constraints (2) the most reliable reconstruction (estimation) $`\widehat{\rho }_r`$ is the one with the maximal value of the von Neumann entropy $`S(\widehat{\rho })=\mathrm{Tr}(\widehat{\rho }\mathrm{ln}\widehat{\rho })`$:

$`S(\widehat{\rho }_r)=\mathrm{max}[S(\widehat{\rho }_{\{\widehat{G}\}});\widehat{\rho }_{\{\widehat{G}\}}].`$ (3)
As shown by E.Jaynes  The operator which fulfills the constraints (2) and simultaneously maximizes the von Neumann entropy can be expressed in the generalized canonical form

$`\widehat{\rho }_r={\displaystyle \frac{1}{Z_{\{\widehat{G}\}}}}\mathrm{exp}({\displaystyle \underset{\nu }{}}\lambda _\nu \widehat{G}_\nu );`$ (4)
where

$`Z_{\{\widehat{G}\}}(\lambda _1,\mathrm{},\lambda _n)=\mathrm{Tr}[\mathrm{exp}({\displaystyle \underset{\nu }{}}\lambda _\nu \widehat{G}_\nu )],`$ (5)
is the generalized partition function and $`\lambda _\nu `$ are the Lagrange multipliers. The Lagrange multipliers $`\lambda _\nu `$ are chosen so that the density operator (4) fulfills the constraints (2) imposed by the results of the measurement. It is then obvious that the meanvalues $`\overline{G}_\nu `$ which determine the density operator are related to the Lagrange multipliers via the derivatives of the partition function

$`\overline{G}_\nu =\mathrm{Tr}(\widehat{\rho }_r\widehat{G}_\nu )={\displaystyle \frac{}{\lambda _\nu }}\mathrm{ln}Z_{\{\widehat{G}\}}(\lambda _1,\mathrm{},\lambda _n).`$ (6)
If we solve the last equation with respect to the Lagrange multipliers we can express them in terms of the measured meanvalues

$`\lambda _\nu =\lambda _\nu (\overline{G}_1,\mathrm{},\overline{G}_n).`$ (7)
When we substitute the Lagrange multipliers (7) into the expression for the generalized canonical density operator (4) we obtain the explicit expression for the reconstructed (estimated) density operator.

As a typical example of the application of the MaxEnt principle we can consider a measurement of a single-mode electromagnetic field, modeled as a harmonic oscillator. Imagine, that as a result of the measurement we know the mean photon number $`\overline{n}`$ in the given field mode. Certainly, there are (infinitely) many quantum states of a single-mode electromagnetic field (e.g., a Fock state, a coherent state, a squeezed state, etc.) with the given mean photon number. So the question is: Which is the best (most reliable) estimation of the measured state given the mean photon number is known? The mean photon number is in some sense the least available information about the measured state. Consequently, the state is least determined. On the other, pure states are completely determined which is reflected by the fact that they have a zero von Neumann entropy. Therefore we expect that the most reliable reconstruction in the given case is a statistical mixture which is determined just by a single parameter - the mean photon number. It is well known that a statistical mixture which is parametrized just by a single parameter is a thermal state. This statistical mixture of Fock states is characterized by a temperature, or the corresponding mean photon number. In addition, for a given temperature (mean photon number) the thermal state exhibits the largest von Neumann entropy. Consequently, from the Jaynes principle of the maximum entropy it follows that if from a measurement only a mean photon number is known that the most reliable estimation of the measured state is the thermal state.

The MaxEnt principle is not the only criterion how to choose an appropriate density operator among those $`\widehat{\rho }_{\{\widehat{G}\}}`$ which fulfill the constraints (2). Based on an intuition or some addition a priori knowledge one can apply other criteria. For example, the maximum likelihood principle has been adopted successfully for estimation of quantum states . Although this reconstruction scheme can result in nonphysical estimations (e.g., density operators which are not properly normalized, etc). In general, a consistent reconstruction scheme has to avoid nonphysical results (e.g. occurrence of negative probabilities). In the MaxEnt reconstruction a physical estimate is guaranteed by the canonical form of the density operator (4). The MaxEnt principle is the most conservative assignment in the sense that it does not permit one to draw any conclusions not warranted by the data. From this point of view the MaxEnt principle has a very close relation (or can be understood as the generalization) of the Laplace’s principle of indifference which states that where nothing is known one should choose a constant valued function to reflect this ignorance. Then it is just a question how to quantify a degree of this ignorance. If we choose an entropy to quantify the ignorance, then the relation between the Laplace’s indifference principle and the Jaynes principle of the Maximum Entropy is transparent, i.e. for a constant-valued probability distribution the entropy takes its maximum value.

The MaxEnt reconstruction has been applied for various quantum systems, such as light field mode, spin systems . In what follows we adopt it for the reconstruction of vibrational states of neutral atoms. We assume the experimental setup realized by the group of C. Salomon at the ENS in Paris .

## III Reconstruction of motional states of neutral atoms

Recently, experimental manipulations of motional quantum states of neutral atoms have been reported by the group of C. Salomon in Paris . Cold Cs atoms can be cooled into specific quantum states of a far detuned 1D optical lattice. The optical lattice is induced by the interference of two laser beams. Along the vertical $`z`$ axis a periodic potential of “harmonic” micro-wells is produced with a period of 665 nm and with an amplitude of about 0.2 $`\mu `$K . The vertical oscillation frequency in a micro-well at the center of the trap is $`\omega _z/2\pi =85`$ kHz. The corresponding ground state has the rms size $`\mathrm{\Delta }z_0=\sqrt{\mathrm{}/2m\omega _z}21`$ nm and $`\mathrm{\Delta }p_0/m=\sqrt{\mathrm{}\omega _z/2m}11`$ mm/s is its rms velocity width. The trapped cloud of neutral Cs atoms has a nearly Gaussian shape with a vertical rms size $`\mathrm{\Delta }\xi _0=53\mu `$m. With the help of deterministic manipulations the neutral atoms can be prepared in non-classical 1D motional states along the vertical axis such as squeezed states, number states, or specific superpositions of number states . The measurement of the prepared quantum state $`\widehat{\rho }`$ is performed as follows: The system is evolved within the harmonic potential during the time $`\tau `$. Then the lasers are turned off and the system undergoes the ballistic expansion (BE). After the time of flight $`T=8.7`$ ms a 2D absorption image of the cloud is taken in 50 $`\mu `$s with a horizontal beam . Integration of 2D absorption images in the horizontal direction gives us the spatial distribution along the vertical $`z`$ axis. Therefore we will consider only 1D quantum-mechanical system along the vertical axis.

To confirm that a desired quantum state has been obtained (engineered) one can compare the spatial distributions along the vertical axis with the predicted ones. The coincidence of these spatial distributions is a necessary but not the sufficient requirement. A complete verification of the fidelity of the preparation of desired quantum states requires a quantum state reconstruction procedure. In order to perform this task we adopt the MaxEnt principle . To do so we utilize a close analogy between quantum homodyne tomography  and the BE absorption imaging for the case of the point-like cloud (with the rms size equal to zero).

### A Quantum tomography via MaxEnt principle

Quantum tomography is based on the inverse Radon transformation of the the measured probability density distributions $`w_{\widehat{\rho }}(x_\theta )`$ for rotated quadratures $`\widehat{x}_\theta =\frac{1}{\sqrt{2}}\left(\widehat{a}\mathrm{e}^{i\theta }+\widehat{a}^{}\mathrm{e}^{i\theta }\right)`$ . These distributions can be represented as a result of the measurement of the continuous set of projectors $`|x_\theta x_\theta |`$. Based on the measurement of the distributions $`w_{\widehat{\rho }}(x_\theta )`$ for all values of $`\theta [0,\pi ]`$ we can formally reconstruct the density operator according to the formula

$`\widehat{\rho }_r={\displaystyle \frac{1}{Z_0}}\mathrm{exp}\left[{\displaystyle _0^\pi }𝑑\theta {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x_\theta |x_\theta x_\theta |\lambda (x_\theta )\right],`$ (8)
where the Lagrange multipliers $`\lambda (x_\theta )`$ are given by an infinite set of equations

$`w_{\widehat{\rho }}(x_\theta )=\sqrt{2\pi }x_\theta |\widehat{\rho }_r|x_\theta ;x_\theta (\mathrm{},\mathrm{}).`$ (9)
If the distributions $`w_{\widehat{\rho }}(x_\theta )`$ are measured for all values of $`x_\theta `$ and all angles $`\theta `$ then the density operator $`\widehat{\rho }_r`$ is reconstructed precisely and is equal to density operator obtained with the help of the inverse Radon transformation or with the help of the pattern functions (for more details see ).

In practical experimental situation (e.g., see the experiments by Raymer et al.  and by Mlynek et al. ) it is impossible to measure the distributions $`w_{\widehat{\rho }}(x_\theta )`$ for all values of $`x_\theta `$ and all angles $`\theta `$. What is measured are distributions (histograms) for finite number $`N_\theta `$ quadrature angles $`\theta `$ and the finite number $`N_x`$ of “bins” for quadrature operators. This means that practical experiments are associated with an observation level specified by a finite number of observables

$`\widehat{F}_{jk}=|x_{\theta _k}^{(j)}x_{\theta _k}^{(j)}|`$ (10)
with the number of quadrature angles equal to $`N_\theta `$ and the number of bins for each quadrature equal to $`N_x`$. We can therefore assume that from the measurement of the observables $`\widehat{F}_{jk}`$ the meanvalues $`\overline{F}_{jk}`$ are determined (these meanvalues correspond to “discretized” quadrature distributions). In addition it is usually the case that the mean excitation number of the state is known (measured) as well.

The operators $`\widehat{F}_{jk}`$ together with $`\widehat{n}`$ form a specific observation level corresponding to the incomplete tomographic measurement. In this case we can express the generalized canonical density operator in the form

$$\widehat{\rho }_r=\frac{1}{Z}\mathrm{exp}\left(\lambda _n\widehat{n}\underset{j=1}{\overset{N_x}{}}\underset{k=1}{\overset{N_\theta }{}}\lambda _{j,k}\widehat{F}_{jk}\right)$$
(11)
The knowledge of the mean photon number is essential for the MaxEnt reconstruction because it formally regularizes the MaxEnt reconstruction scheme (the generalized partition function is finite in this case).

### B Motional states of atoms via MaxEnt principle: Formalism

In the quantum homodyne tomography the probability distributions are measured for the rotated quadrature operators $`\widehat{x}_\theta `$. The annihilation and creation operators of motional quanta, $`\widehat{a}`$ and $`\widehat{a}^{}`$, are related to the position and momentum operators, $`\widehat{z}`$ and $`\widehat{p}`$, via expressions $`\widehat{z}=\frac{1}{\sqrt{2}}(\widehat{a}+\widehat{a}^{})`$ and $`\widehat{p}=\frac{1}{\sqrt{2}}i(\widehat{a}\widehat{a}^{})`$, respectively. The angle $`\theta `$ of the quadrature operator corresponds to $`\omega _z\tau `$ and vertical “cuts” of the absorption images (taken after the BE) can be associated with quadrature probability distributions. However, for a real physical situation with a nonzero rms size of the cloud the vertical “cuts” of absorption images correspond to a coarse-grained quadrature probability distributions. In particular, the vertical cuts of measured absorption images (taken in 2D) give us (after integration along the horizontal direction) the spatial distribution along the vertical axis. The spatial distribution along the vertical $`z`$ axis can be expressed as

$`\overline{F}_\tau (z)=T^1{\displaystyle F_0(\xi _0)P_\tau ((z\xi _0)/T)𝑑\xi _0},`$ (12)
where $`F_0(\xi _0)`$ is the initial spatial distribution of the cloud in the $`z`$-direction (i.e., a Gaussian distribution with the rms size $`\mathrm{\Delta }\xi _0`$). The function $`P_\tau (v)`$ denotes the velocity probability distribution of the measured quantum state which has been evolved for time $`\tau `$ in the harmonic potential before the BE, i.e.

$`P_\tau (v)=|v|\psi (\tau )|^2,|\psi (\tau )=\widehat{U}(\tau )|\psi (0).`$ (13)
Here $`\widehat{U}(\tau )=\mathrm{exp}(i\widehat{H}\tau /\mathrm{})`$ represents the time-evolution operator for the harmonic oscillator with the Hamiltonian $`\widehat{H}=\widehat{p}^2/2m+m\omega _z^2\widehat{z}^2/2`$. Now we can treat the measured “cuts” as the mean values of specific observables: $`\overline{F}_\tau (z)=\text{Tr }[\widehat{\rho }\widehat{F}_\tau (z)]`$. In practice just a few discrete times $`\tau _j`$ ($`j=1,\mathrm{},N_\tau `$) are considered and the $`z`$ coordinate is discretized into the bins $`z_k`$ ($`k=N_z,\mathrm{},N_z`$) of a given resolution $`\mathrm{\Delta }z`$. The set of operators which enters the equation (11) for the MaxEnt reconstruction then takes the form

$`\widehat{F}_{\tau _j}(z_k)=T^1{\displaystyle F_0(\xi _0)\widehat{U}^{}(\tau _j)|\frac{z_k\xi _0}{T}\frac{z_k\xi _0}{T}|\widehat{U}(\tau _j)𝑑\xi _0}(j=1,\mathrm{},N_\tau ;k=N_z,\mathrm{},N_z).`$ (14)
We have already commented that the operator of mean phonon number $`\widehat{n}`$ is added to the set of observables $`\{\widehat{F}_{\tau _j}(z_k)\}`$. Knowledge of the mean excitation number $`\overline{n}`$ is essential in the case of an incomplete set of observables . Knowledge of the mean excitation number leads to a natural “truncation” of the Hilbert space. The inclusion of the mean phonon number into the MaxEnt reconstruction scheme does not represent its limitation as the mean energy represents one of basic characteristics of any system which should be inferred from the measurement.

The experimental “cuts” of the BE absorption images \[$`\widehat{F}_{\tau _j}(z)`$\] can be taken at few selected times, for example $`\omega _z\tau _j=0,\pi /4,\pi /2,3\pi /4`$ ($`N_\tau =4`$). To perform the reconstruction we have to determine the Lagrange multipliers $`\{\lambda _{j,k}\}`$ and $`\lambda _n`$ associated with $`\{\widehat{F}_{\tau _j}(z_k)\}`$ and $`\widehat{n}`$, respectively, in the expression for the generalized canonical density operator (11). The Lagrange multipliers can be determined via the minimization of a deviation function $`\mathrm{\Delta }F`$ with respect to the measured data, i.e.

$`\mathrm{\Delta }F`$ $`=`$ $`{\displaystyle \underset{j,k}{}}w_{j,k}\left\{\overline{F}_{\tau _j}(z_k)\mathrm{Tr}\left(\widehat{\rho }_r\widehat{F}_{\tau _j}(z_k)\right)\right\}^2`$ (15)
$`+`$ $`w_{\overline{n}}\left\{\overline{n}\mathrm{Tr}\left(\widehat{\rho }_r\widehat{n}\right)\right\}^2.`$ (16)
Here $`\{w_{j,k}\}`$ and $`w_{\overline{n}}`$ represent positive weight factors for particular observables. Without any prior knowledge about the state we can take for simplicity $`w_{i,j}=1`$. The weight factor $`w_{\overline{n}}`$ associated with the mean phonon number can be chosen according to our preference either to fit better the “cuts” of the BE images or the mean phonon number. In the case of the perfect measurement and the complete reconstruction the result has to be independent of the choice of the weight factors (in this case we can take $`w_{\overline{n}}=1`$). The weight factors could be also associated with the prior information about the dispersion of the measured observables. In particular, the weight factors can be taken as $`w_\nu \sigma _\nu ^2`$ to reflect the knowledge of variances $`\sigma _\nu `$ for the measured observables $`\widehat{G}_\nu `$. When the mean values of the observables for the MaxEnt estimate $`\widehat{\rho }_r`$ fit within desired interval $`\overline{G}_\nu \pm \sigma _\nu `$ then contributions of the observables to the deviation function $`\mathrm{\Delta }F`$ are of the same order ($`1`$). However, in our case we do not assume the knowledge of variances for the measured discretized probability distributions (taking $`w_\nu =1`$).

Once the Lagrange multipliers are numerically fitted, the result of the reconstruction – the generalized canonical density operator $`\widehat{\rho }_r`$ – can be visualized, for example, via the corresponding Wigner function  which can be defined as a particular Fourier transform of the density operator $`\widehat{\rho }`$ of a harmonic oscillator expressed in the basis of the eigenvectors $`|q`$ of the position operator $`\widehat{q}`$:

$`W_{\widehat{\rho }}(q,p){\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\zeta q\zeta /2|\widehat{\rho }|q+\zeta /2\mathrm{e}^{ip\zeta }.`$ (17)
### C Numerical simulation

To test our reconstruction procedure let us consider the reconstruction of the Wigner function of the motional quantum state $`|\psi (0)=(|0+|1)/\sqrt{2}`$ of Cs atoms trapped in the optical lattice. This kind of states has been demonstrated in recent experiments . We assume the following setup parameters: $`\omega _z/2\pi =80`$ kHz, the rms size of the ground state $`\mathrm{\Delta }z_0=22`$ nm, the rms velocity width $`\mathrm{\Delta }p_0/m=11`$ mm/s and the rms width of the cloud of the atoms about 60 $`\mu `$m. Before BE (with BE time $`T=8.7`$ ms) the atoms evolve within the harmonic trapping potential for $`\tau =0,1.6,3.2,4.8`$ $`\mu `$s. As the input for the reconstruction via the MaxEnt principle four vertical “ideal” cuts of the BE absorption images are taken as shown in insets of Fig. 1. In addition, for the phonon number operator $`\widehat{n}`$ which is included in the set of measured observables (see discussion above) we assume the mean value $`\overline{n}=0.5`$. The result of the ideal reconstruction is shown in Fig. 1. The fidelity of the measured and the reconstructed quantum states is close to unity which means a perfect reconstruction with $`\mathrm{\Delta }F=10^{10}`$, entropy $`S=10^7`$, $`\mathrm{\Delta }\rho =10^8`$ has been achieved. Here $`\mathrm{\Delta }\rho =_{m,n}|(\widehat{\rho }\widehat{\rho }_r)_{mn}|^2`$ denotes a deviation of the original and reconstructed density operators.

Obviously, in a real measurement the measured values are always fluctuating around the exact ones due to an experimental noise. Therefore we simulate a non-ideal measurement introducing random fluctuations to the measured values of observables. It means that instead of the ideal values $`\overline{F}_{\tau _j}(z_k)`$ we use for the MaxEnt reconstruction procedure the fluctuating (“noisy”) values

$`\overline{F}_{\tau _j}^{}(z_k)=\overline{F}_{\tau _j}(z_k)+\eta \xi _{j,k}\left(\overline{F}_{\tau _j}(z_k)\right)^{1/2}.`$ (18)
Here $`\eta `$ is a relative-error parameter which characterizes the quality of the measurement and $`\{\xi _{j,k}\}`$ represents a Gaussian noise for observables. The result of the reconstruction is shown in Fig. 2 for $`\eta =0.1`$. Noisy mean values of the observables are shown in the insets. Despite of a significant relative error the reconstruction is almost perfect with the fidelity of the measured and the reconstructed states still close to one ($`\mathrm{\Delta }F=0.16`$, entropy $`S=0.01`$, $`\mathrm{\Delta }\rho =0.05`$). The minimum value of the deviation function $`\mathrm{\Delta }F=0.16`$ can serve also as a measure of the imperfection of the given measurement (due to a technical noise) .

A typical non-classical state which we can utilize for a further test is the even coherent state $`𝒩_e(|\alpha +|\alpha )`$ which is a superposition of two coherent states with opposite phases . For the amplitude $`\alpha =\sqrt{2}`$ we obtained $`\mathrm{\Delta }F=10^8`$, the entropy $`S=0.026`$ and $`\mathrm{\Delta }\rho =10^4`$ (under assumption that the exact mean phonon number $`\overline{n}=1.928`$ is known). In the case of the imperfect measurement with $`\eta =0.1`$ the reconstruction leads to $`\mathrm{\Delta }F=0.14`$, entropy $`S=0.13`$, and $`\mathrm{\Delta }\rho =0.06`$ for $`\overline{n}=2.09`$. The fidelity of the reconstructed and the measured states is in this case also close to one.

In order to model a technical noise in the measurement we have been considered Gaussian fluctuations proportional to the square root of the mean values. It means that tails of the “cuts” of BE images do not introduce a significant error (compare insets in Fig. 1 and Fig. 2). However, in the current measurements the situation seems to be different and the fluctuations do not decrease with the amplitude of the expected values.

The fundamental question in the context of the MaxEnt reconstruction of states from incomplete tomographic data is whether the quality of the reconstruction can be improved using additional data from subsequent time moments $`\tau `$ and how many such time moments $`\tau `$ are required for the complete reconstruction of the unknown state $`\widehat{\rho }`$. We have shown recently  that for the quantum tomography just three quadrature distribution are sufficient for a complete reconstruction using the MaxEnt principle (in the case of the perfect measurement). This corresponds to the ideal case without the spatial dispersion of the cloud of atoms, i.e. the choice with $`\omega _z\tau _j=0,\pi /4,\pi /2`$ ($`N_\tau =3`$) is sufficient for $`\mathrm{\Delta }\xi _00`$. Obviously, in experiments with neutral atoms the spatial size of the atomic cloud is nonzero. However, in the case of the ideal measurement three BE absorption images associated with three “rotations” $`\omega _z\tau _j`$ are still sufficient for a complete reconstruction of tested examples of quantum states. On the other hand, it seems that for higher mean phonon numbers the spatial distributions along the vertical axis which are directly determined from absorption images should be known with improving precision (and on a wider interval of values as well). In the above examples we have considered for convenience BE images for four “rotations” ($`N_\tau =4`$) which results in a very good reconstruction.

## IV Reconstruction from experimental data

In what follows we will apply the MaxEnt reconstruction scheme to the data obtained at the ENS in Paris . Firstly we note that the unknown quantum state should belong to a Hilbert subspace which can be determined easily. Thus we can limit ourselves to the subspace spanned by Fock (number) states $`|0,|1,\mathrm{}|N1`$. The upper bound on the accessible phonon number $`N`$ is given by experimental limitations such as, for example, a feasible depth of micro-wells of the optical lattice and the validity of the harmonic potential approximation. For recent experiments $`N`$ has been typically of the order of 10. This value is large enough to demonstrate the preparation of many non-classical states but on the other hand excludes highly squeezed states from a coherent processing.

Let us consider the experimental arrangement used in Paris  with the parameters: $`\omega _z/2\pi =80`$ kHz, the rms size of the ground state $`\mathrm{\Delta }z_0=22`$ nm, the rms velocity width $`\mathrm{\Delta }p_0/m=11`$ mm/s, the rms width of the cloud of the atoms about 60 $`\mu `$m and BE time $`T=8.7`$ ms. Initially the atoms are prepared in a well defined motional state $`|\psi _0`$ (e.g. in the vacuum state $`|0`$). Then the optical lattice is switched off for the time period $`t_1`$ during which the atoms evolve freely towards the state $`|\psi _1=\mathrm{exp}(it_1\widehat{p}^2/2m)|\psi _0`$. Next, the optical lattice is again switched on for the time $`\tau `$ during which the atoms evolve within the harmonic trapping potential. The measurement is performed after the BE time. The first two stages can be considered as the preparation of the state $`|\psi _1`$. After its “rotation” by $`\omega _z\tau `$ (within the phase space of the harmonic oscillator) and the subsequent BE the absorption images are taken.

The considered data are for the initial vacuum state which means that under ideal conditions a squeezed state $`|\psi _1=\mathrm{exp}(it_1\widehat{p}^2/2m)|0`$ should be prepared. The vertical spatial distributions obtained from the measured 2D-absorption images are discretized into pixels (bins) with the pixel width 5.45 $`\mu `$m. The optical density of each pixel is averaged in the horizontal direction in which the absorption images are divided into 50 rows, each 3.9 $`\mu `$m wide (these rows cover the size of the cloud in the horizontal direction). For the reconstruction via the MaxEnt principle four vertical spatial distributions for “rotation” times $`\tau =0,1.6,3.2`$ and $`4.8\mu `$s are taken. The selected times roughly correspond to rotations within the phase space by $`\omega _z\tau =0,\pi /4,\pi /2`$ and $`3\pi /4`$, respectively. Unfortunately, the mean excitation number $`\widehat{n}`$ for measured state $`|\psi _1`$ was not measured explicitly in the experiment, therefore we have to estimate it as follows: During the free expansion period the rms size of the cloud increases by $`\mathrm{\Delta }x=p_0\tau _1/m`$. The corresponding increase of the potential energy $`\frac{1}{2}m\omega _z^2(\mathrm{\Delta }x)^2`$ in units $`\mathrm{}\omega _z`$ gives us the increase of the number of excitation quanta with respect to the initial state $`|\psi _0`$. For $`\tau _1=4`$ $`\mu `$ and the initial vacuum it means $`\overline{n}1`$. Experiments can be realized also for higher $`\tau _1`$. For example, $`\tau _1=8`$ $`\mu `$ leads to $`\overline{n}4`$. However, as mentioned above, such “squeezed” states with a significant contribution of higher phonon number states violate the underlying harmonic approximation for the potential. To keep a coherent control an an-harmonic part of the potential has to be taken into account.

The result of the reconstruction via the MaxEnt principle is shown in Fig. 5. The deviation of the fitted and measured values is $`\mathrm{\Delta }F=0.09`$ and the entropy of the reconstructed state $`S=1.0`$. It means that the reconstructed state is a statistical mixture. We see a two peak structure which suggests that there is a mixture of two squeezed states coherently displaced from each other. It is caused by the fact that the vertical center of the cloud was not fixed in the experiment and it has to be determined by our fit for each measured BE absorption image separately. Assuming a priori knowledge that the Wigner function has a symmetric shape with respect to the origin of the phase space (i.e. there is no coherent amplitude) a Gaussian fit can be used to determine the center of the cloud for each vertical distribution. For states with a non-zero coherent amplitude the center of the cloud should be fixed already in the experiment.

It turns out that the reconstruction results do not describe the squeezed vacuum state as was originally expected . The main reason is that the mean phonon number was not measured directly in the experiment. It can be inferred only indirectly from the ideal case without any incoherence during preparation or measurement. As we discussed above, it is essential to include the information about the mean number of vibrational quanta into the MaxEnt reconstruction scheme. In optical tomography the analogous information about mean photon number can be obtained from distributions of two “orthogonal” quadratures. In our case it could correspond to two absorption images such that $`\omega _z(\tau _j\tau _k)=\pi /2`$. However, it would require a precise timing of the evolution within the harmonic trapping potential. Therefore the mean number of vibrational quanta should be determined in an independent measurement.

Another problem arises from a slow convergence of anti-squeezed spatial distributions which are derived directly from the measured absorption images. In particular, the convergence of tails is too slow for those “rotations” which corresponds to anti-squeezed phases, i.e. $`\tau =1.6`$, $`4.8\mu \text{s}`$ (see insets of Fig. 5). The slow convergence is reflected by the presence of non-negligible backgrounds for Gaussian fits to these spatial distributions. If we eliminate (subtract) these backgrounds from the measured distributions the MaxEnt reconstruction gives almost the same Wigner function as in Fig. 5 but with a highly reduced deviation function $`\mathrm{\Delta }F=0.02`$ (comparing to $`\mathrm{\Delta }F=0.09`$ in Fig. 5). Such background in these absorption images can be caused by an incoherence associated with a violation of the the harmonic approximation. In fact, in our analysis we have neglected the change of the oscillation frequency along the $`z`$-axis. In recent experiments, the oscillation frequency decreases 10% from $`\omega _z`$ for micro-wells at the edge of the initial cloud.

## V Conclusions

We have presented a very efficient reconstruction scheme for the reconstruction of motional states of atoms based on the MaxEnt principle. The main advantage of the scheme is that it always results in reconstructions which are physical states (unlike in the case of the maximum likelihood estimation which can results in nonphysical estimations). Moreover, the scheme is very efficient in a sense that it requires just a small number of tomographic phase-space cuts.

We have applied the scheme for a reconstruction of motional quantum states of neutral atoms. As an example have we analyzed the experimental data obtained by by the group of C. Salomon at the ENS in Paris and we reconstruct the Wigner function of motional quantum states of Cs atoms trapped in the optical lattice. In our analysis we have neglected the change of the oscillation frequency along $`z`$ axis in recent experiments. The dispersion of the oscillation frequency is of the order of a few percent. This source of errors can significantly affect the quality of a quantum state preparation and its reconstruction. In addition, only up to first 10 bound states of micro-wells of the optical lattice can be approximated by a harmonic potential. It implies limits on coherent manipulations of quantum states. It means that states with a significant contribution of higher number (Fock) states cannot be prepared and manipulated in a controlled way.

###### Acknowledgements.

This work was supported by the European Union projects QUBITS and QUEST under the contracts IST-1999-13021 and HPRN-CT-2000-00121, respectively. We thank Christophe Salomon and Isabelle Bouchoule for providing us with their experimental data and for helpful discussions and correspondence.


text:
# Competition Among Companies: Coexistence and Extinction

## I Introduction

During the last few years we have witnessed a wealth of work on the application of methods of statistical physics to the study of economic problems configuring what some authors called econophysics uno ; dos ; tres ; cuatro . Within this framework a great deal of effort was dedicated to the analysis of economic data ranging from stock exchange fluctuations uno , production models, size distribution of companies, the appearance of money, effects of control on the market to market critical properties cinco . Another problem that attracted enormous interest was the origin of power (Paretto) laws, and lognormal distribution with power law tails, for the income of individuals, wealth distribution, debt of bankrupt companies SS . An interesting source of several mathematical descriptions and models used in economical and sociological contexts can be found in haag .

Here, our interest is the study, in a deterministic way, of aspects of the competition and coexistence of agents or companies in a common market. We present a simple “toy” model describing, in analogy to some ecological problems, a situation of competition among several companies. Within a “Malthusian-like” model we analyze the effect of a kind of intracompetitive contribution on the possibility of companies coexistence in a certain market, and the changing leadership role (measure through some wealth parameter) between, according to the standard ecological definition, “strong” and “weak” companies.

According to ecological studies, starting with Volterra’s first results on the mathematical theory of competition vol , the problem of competition and coexistence between species has been analyzed and resumed within the Competitive Exclusion Principle (or Ecological Theorem), that states: $`N`$ species that compete for $`n(<N)`$ food resources, cannot coexist mur . Several aspects of this problem have been analyzed by different authors, emphasizing, for instance, the conflict between the need to forage and the need to avoid competition; effects of diffusion-mediated persistence. Generally, the system describing competition between species can be represented by a set of $`N`$ differential equations for the species and $`n`$ for the resources. As an example, with only one food resource we find only two stable stationary solutions: the trivial one (extinction of all he species), and that corresponding to the survival of only one species, the “strongest” one. There are also studies of the problem related to the possibility of coexistence in the form of wave–like solutions sch ; wkvh .

Here we adapt the model used in Refs.sch ; wkvh  for the case of a homogeneous competitive market with a unique wealth resource and several firms. In the next Section we introduce the model and some particular, instructive, solutions. Due to the difficulties of finding analytical solutions for the general case, in Section III we focus on a representative system with a small number (here 4) of firms, and analyze its behavior by numerical methods. In the last Section we draw some conclusions.

## II The Model

We start describing the model we use, which is related to the one used in Refs.sch ; wkvh  for the study of coexistence in an ecological framework in the form of wave–like solutions. Such a model has been adapted to the problem of competition of $`N`$ agents or companies for a unique common wealth resource. We indicate with $`n_j`$ the “size” or wealth-parameter representing the welfare of the $`j`$-company ($`j=1,\mathrm{},N`$) and with $`M`$ the total amount of “wealth”. The set of differential equations that we use to describe the behavior (for the homogeneous case) of such a complex system includes a “Malthusian-like” birth-death equation mal  for each company ($`\beta _jM(t)n_j(t)`$ corresponding to benefits coming from the wealth’s share and $`\alpha _jn_j(t)`$ for the “standard” losses–or costs– of the $`j`$-th company). We also include a contribution that corresponds, in ecological language, to taking into account the existence of some kind of intracompetition, that is if the $`j`$-th the company is alone and can get all the wealth $`M`$, it can only grow up to a maximum bounded size, a behavior that can be modelled by a Verhulst-like term mur . This takes into account the in ecological contexts so called carrying capacity of the economic environment ccap ; yad . In economic terms it could correspond to an increase of the company efficiency, for instance, through the reduction of operative internal costs, improved management, avoiding of competition between different branches of the same firm, instrumentation of new technologies, etc.

In addition, instead of assuming a constant finite wealth resource as in the so called “Total Wealth Conserved model” piav , we consider that $`M`$ has its own dynamics. For the equation for $`M`$, the economic wealth accessible to the companies that, in order to simplify this initial analysis we assume is unique, we also consider a “Malthusian-like” behavior including: its production (new resources and technologies, harvest and grain production, etc) that we assume has a constant rate $`Q`$, and its disappearance due to: (a) natural degradation or rotting of crops, technologies becoming old, some resources being exhausted, that we assume has a rate proportional to the total wealth amount $`GM`$ (only a certain portion of $`M`$ disappears), (b) the share of each company given by $`\beta _kn_k(t)M(t)`$.

The set of equations is ($`j=1,\mathrm{},N`$)

$`{\displaystyle \frac{d}{dt}}n_j(t)`$ $`=`$ $`[\beta _jM(t)\alpha _j]n_j(t)\gamma _j{\displaystyle \frac{n_j^2}{M}}`$
$`{\displaystyle \frac{d}{dt}}M(t)`$ $`=`$ $`Q(t)[G+{\displaystyle \beta _kn_k(t)}]M(t).`$ (1)
Similarly to what was discussed in wkvh  for the case of only two species, we can here define a hierarchy from the “strongest” (largest ratio between wealth share and standard losses, i.e. largest $`\frac{\beta _i}{\alpha _i}`$) to the “weakest” (smallest ratio) companies. Assuming the following hierarchical order

$$\frac{\beta _1}{\alpha _1}>\frac{\beta _2}{\alpha _2}>\mathrm{}>\frac{\beta _N}{\alpha _N},$$
(2)
we have that $`n_1`$ is the strongest company while $`n_N`$ is the weakest one. It is worth noting in passing that Eqs.(II) resembles the form of multimode laser systems lasers , making it possible to transfer some results from one system to the other.

The stationary solutions results from taking $`\frac{d}{dt}M=0`$ and $`\frac{d}{dt}n_j=0`$ ($`j=1,\mathrm{},N`$). We found

$`M_s`$ $`=`$ $`{\displaystyle \frac{Q(t)}{G+\beta _kn_k^s}}`$ (3)
$`0`$ $`=`$ $`\left([\beta _jM_s\alpha _j]\gamma _j{\displaystyle \frac{n_j^s}{M_s}}\right)n_j^s.`$ (4)
The last equation implies one of two possibilities

$$n_j^s=0,\text{or}n_j^s=\frac{[\beta _jM_s\alpha _j]M_s}{\gamma _j}.$$
(5)
It is clear that for large $`N`$, to find the solution of this system is not easy. In order to fix ideas we consider the simplified case where, instead of the above indicated hierarchy, we have that all companies are equivalent, that is

$$\beta _j=\beta ,\alpha _j=\alpha ,\gamma _j=\gamma ,$$

implying

$$n_j^s=n_sj.$$

In this case we have

$$n_s=\frac{[\beta Q(\alpha /Q)(G+N\beta n_s)]Q^2}{\gamma (G+N\beta n_s)^2},$$
(6)
that can be rewritten as

$$\frac{\gamma N^2\beta }{Q^2}n_s^3\frac{2\gamma GN\beta }{Q^2}n_s^2\left(\frac{\alpha \beta N}{Q}+\frac{\gamma G^2}{Q^2}\right)n_s+(\beta \frac{\alpha G}{Q})=0.$$
(7)
It is possible to find under which conditions at least one solution of Eq. (7) is real. However, it is more instructive to look for the behavior at small $`n_s`$ ($`n_s0`$) as all the coefficients of $`n_s^\nu `$ with $`\nu >0`$ are negative. Hence, we can easily obtain that a solution, given by

$$n_s\frac{\beta Q^2\alpha GQ}{\alpha \beta NQ+\gamma G^2},$$

exists (is positive) if $`\beta Q>\alpha G.`$ In this case, as one of the associated eigenvalues is zero, a linear stability analysis does not give a clear information about the stability of the solution and it is necessary to resort to a more refined analysis.

Another instructive case is to consider

$$\frac{\beta _1}{\alpha _1}>\frac{\beta _2}{\alpha _2}=\mathrm{}=\frac{\beta _N}{\alpha _N}.$$

Here we reduce to essentially the same situation studied in sch ; wkvh . In particular, it is coincident with the situation studied in KW , but now having an “effective” weak species given by $`(N1)n_j`$, $`j=2,N`$. As in KW , and as discussed in detail latter for the case of several firms, it is possible to find a coexistence region when $`\gamma _1`$ overcomes some threshold value. In this case a linear stability analysis shows a change in the stability of these solutions.

As a general analytical study of our system even for $`N`$ not too large is far beyond our interest, in the next section we focus on a numerical approach for a case with a small, however representative, value of $`N`$ analyzing some relevant situations.

## III Numerical Results

As indicated before, here we focus in a case with $`N`$ small (in fact $`N=4`$) that shows all the relevant aspects we can expect in the large $`N`$ situation. Defining the ratio $`\rho _i=\frac{\beta _i}{\alpha _i}`$ we consider several situations:
a) when the different companies are in hierarchical order, that is $`\rho _1>\rho _2>\rho _3>\rho _4`$;
b) when we have $`\rho _1>\rho _2=\rho _3=\rho _4`$;
c) when $`\rho _1=\rho _2>\rho _3=\rho _4`$.

Among all the possible scenarios we have chosen those showing regions of coexistence of all the species, thus in case (a) we investigate the situation with fixed $`\gamma _j`$, with $`j=2,3,4`$ and varying $`\gamma _1`$. That allows us to have a control parameter, but we recall that the choice is arbitrary. In this case the usual strong and weak concepts make us consider the species 1 as the strongest and the 4 as the weakest. With the inclusion of the new term we find not only the possibility of coexistence but also that the original hierarchical order can be permuted several times as the parameter values are varied. Thus, we find that the ranking of companies suffers many changes, with companies interchanging roles several times. Some examples of this case are shown in Fig.1. Besides these new features we observe the classical extinction of companies as predicted by the exclusion theorem. In all the cases the extinction occurs by one species at a time. Figure 1 shows two typical results for the stationary values reached by $`n_j`$ as function of $`\gamma _1^1`$. It is apparent that according to the different values, the relative status of each company can change and even reversed with several crossing among them. In this figure $`\gamma _1`$ varies continuously from 10 to 0. The ratios $`\rho _i=\frac{\beta _i}{\alpha _i}`$ are $`\rho _1=5`$, $`\rho _2=2.5`$, $`\rho _3=1.5`$, $`\rho _4=1`$. In Fig.1a we have $`\gamma _2=1`$, $`\gamma _3=0.5`$ and $`\gamma _4=0.1`$, while in Fig 1.b $`\gamma _2=1`$, $`\gamma _3=0.1`$ and $`\gamma _4=0.5`$.

For the case (b) above, we considered as our variable $`\gamma _2`$. In this case we have one strong species with $`\rho _1=4`$, and three similar species, $`\rho _2=\rho _3=\rho _4=2`$, that could coexist if there was not a strongest one. We want to note that though $`\rho _2=\rho _3=\rho _4`$ the same is not true for $`\alpha `$ and $`\beta `$ values. Once again, we observe coexistence between species and a reorganization of the company or agents ranking, not accordingly to the original concept of strength but depending of the values of $`\gamma _i`$. The coexistence is achieved within certain parameter region. The extinction is gradual and governed mainly by $`\gamma `$ values. The results are shown in Fig. 2, where again we depict the stationary values of $`n_j`$ as functions of $`\gamma _x^1`$. In Fig 2.a we have $`\gamma _1=0.25`$, $`\gamma _3=0.5`$, $`\gamma _4=0.1`$, while in Fig 2.b $`\gamma _1=1.`$, $`\gamma _3=0.5`$, $`\gamma _4=0.1`$. It is apparent that the most relevant parameter when considering competition is the value of $`\gamma _i`$. When species are equally strong or weak, we observe that different stationary density levels are reached according to $`\gamma _i`$. If $`\gamma _i`$ is of the same order the coexistence is granted. On the contrary, a species with a high $`\gamma _i`$ will not survive even if competing with similar species. As an additional feature we observe that a usual strong species $`j`$ will not survive even if competing with weaker species if $`\gamma _j`$ is much higher than that of the other species.

In case (c) we considered two situations: a first one varying $`\gamma _1`$; and a second varying $`\gamma _3`$. The results are shown in Fig. 3, where again we depict the stationary values of $`n_j`$ as functions of $`\gamma _1^1`$ and $`\gamma _3^1`$ respectively. In both cases 3 $`\rho _1=\rho _2=4`$, $`\rho _3=\rho _4=2`$, while in a) $`\gamma _2=0.25`$, $`\gamma _3=0.001`$, $`\gamma _4=1.5`$ and in b)$`\gamma _1=0.5`$, $`\gamma _2=1.5`$, $`\gamma _4=0.01`$ . We confirm that the coexistence can be achieved for proper $`\gamma `$ values. At the same time we observe that the usual concept of strongest species cannot be applied as in both previous cases, species 2, one of the strongest if $`\gamma _i`$ were zero remains as the weakest species due to a high $`\gamma _2`$ value.

## IV Conclusions

The results shown above indicates how the study of simplified models could help in the understanding of the role played by $`\gamma `$, an internal company parameter, associated to the company efficiency, in situations where the complexity of the economic reality makes very hard to obtain a complete model. In our case, the model so far studied could give some hints on the behavior of systems of companies in competition and the possibility of coexistence and the way that one company can use to “eliminate” the competing ones adopting a policy tending to adequately change its $`\gamma `$. Here we have analyzed the effect of explicitly including the carrying capacity of the environment within our toy model for describing the coexistence of species in competition. The results put in evidence the role played by the term associated to $`\gamma `$ in the possibility of coexistence. We recall that in the absence of such term the exclusion principle is valid and only the species (one or more) with the highest $`\rho `$ survive. It is also this term that, within certain parameter region, governs the company or agents ranking. A model written in the same terms can clearly also be applicable to ecological situations. But in this case, rather than a deterministic control of $`\gamma `$, some fluctuations or cyclic changes in this parameter should be considered. This is the subject of a work of us in progress.

Acknowledgements: Partial support from CONICET and ANPCyT, both Argentina agencies, as well as to Fundación Antorchas, are greatly acknowledged. HSW wants to thank to Iberdrola S.A., Spain, for an award within the Iberdrola Visiting Professor Program in Science and Technology, and to the IMEDEA and Universitat de les Illes Balears, Palma de Mallorca, Spain, for the kind hospitality extended to him.