format: text

text:
# Chaos and flights in the atom-photon interaction in cavity QED

## I Introduction

Cavity quantum electrodynamics is a rapidly developing field of physics studying the interaction of atoms with photons in high-finesse cavities in a wide range of electromagnetic spectrum, from microwaves to visible light, in such conditions under which both atoms and fields may manifest their quantum nature (for reviews on quantum and atom optics, see B ; ASM ). Modern experiments in cavity QED have achieved the exceptional circumstance of strong atom-field coupling for which the strength of the coupling exceeds both the atomic dipole decay and the cavity field decay providing manipulations with single atoms and photons Y99 ; HLD ; EMK ; MF99 ; PFM . Trapped atoms and ions, interacting with laser fields in the regime of the strong coupling, have been used not only to study fundamentals of quantum mechanics W98 ; MM96  but also for applications in the rapidly growing fields of quantum computing, quantum communication (see, for example, CZ95 ; MM95 ), quantum chaos Z81 ; H91 , and decoherence Z91 ; GJ96 .

A special comment should be related to the notion of chaos that we use in the paper. The real system is quantum and the quantum chaos per se does not exist. The notion “quantum chaos” just refers to the behavior of quantum systems whose classical counterparts behave chaotically. Through the paper the notion “chaos” will be used only in that meaning since all our consideration will be semiclassical. The key issue of quantum chaos is the correspondence between the classical and quantum pictures of chaotic dynamics. Chaotic dynamics in the atom-photon physics appeared in the paper Haken , where the semiclassical model of a single-mode, resonant, and homogeneously broadened laser, considered as an open dissipative system, has been shown to be equivalent to a Lorentz-type strange attractor, and the paper BZT , where Hamiltonian semiclassical chaos has been shown to arise in the Dicke model with non-resonant terms describing the interaction of an ensemble of identical two-level atoms with their own radiation field in an ideal resonant cavity. The fully quantum version of the latter model has been considered in BBZ . It was shown in  that in a parameter range, corresponding to classical chaos, the evolution of the system becomes essentially quantum after the so-called “breaking time” $`\tau _{\mathrm{}}`$ BZ78  which obeys in this regime the logarithmic law $`\tau _{\mathrm{}}\lambda ^1\mathrm{ln}(\mathrm{const}/\mathrm{})`$, where $`\lambda `$ is the maximal Lyapunov exponent. All results presented here as chaotic dynamics are valid up to $`t\tau _{\mathrm{}}`$.

The physical mechanism of semiclassical chaos in the Dicke model is tied to virtual transitions which are described by non-resonant (or counter-rotating) terms in the Dicke Hamiltonian. They are small under usual conditions. Trying to find another mechanism of local instability in the Hamiltonian atom-photon dynamics without pump and losses, the authors of PK97 ; PKK99  proposed the semiclassical model with atoms moving through a standing-wave cavity in a direction along which the cavity sustains a space-periodic field. The standing wave modulates the atom-field coupling providing in a certain range of the system’s parameters intermittent Rabi oscillations. Another way to change the atom-field coupling is modulating the cavity length, studied in IKP . New effects in the model with moving atoms may arise beyond the simple semiclassical approximation. Chaotic vacuum Rabi oscillations, a new kind of reversible spontaneous emission, have been shown P00 ; PK00  to occur in the model with interatomic quantum correlations.

Experiments to study the quantum dynamics of classically chaotic systems for the atom-photon interaction in cavities and traps have been intensively studied with cold atoms in a phase-modulated standing wave MR94  and in an amplitude-modulated standing wave HT  (following to the proposition of GSZ ), and in a pulsed standing wave MR95 ; AG . Cold sodium or cesium atoms, which are kicked by a periodically pulsed standing wave of far-detuned light, is an excellent experimental realization of a paradigm model of quantum chaos, a $`\delta `$-kicked quantum rotor. Dynamical localization, that was observed in the atomic momentum distribution, is the quantum suppression of the classical momentum diffusion. A typical underlying phase-space structure of classical chaotic systems consists of stochastic webs, islands, and chains of islands embedded in the stochastic sea. The chaotic motion occupies a certain area in phase space. Because of islands and their boundaries, typical behavior of chaotic systems can be intermittent with long (quasi)regular oscillations (the so-called Lévy flights) interrupted by the chaotic pieces of trajectories. This intermittency leads to the anomalous diffusion with Lévy distribution functions or a similar one, which have power-wise tails (for reviews on Lévy processes in physics see SZK , where the term “strange kinetics” was coined, and Z98 ). It was found with kicked cold cesium atoms K0  that for certain pulses amplitudes, where the respective classical analog may exhibit anomalous diffusion, the momentum distributions were not exponentially localized for the time of observation (see also BS99 ; AI99 ; RA01 ). It should be mentioned that anomalous diffusion and Lévy flights have been even earlier found with cold atoms and employed in a subrecoil laser cooling scheme BB94 ; RB95 ; SL99 .

In experiments with $`\delta `$-kicked atoms, the detuning between the optical and atomic transition frequencies is large (relative to the natural line-width), so the probability is small to find an atom, initially prepared in the ground state, in the excited state, and the excited state amplitude can be adiabatically eliminated GSZ . In this approximation, an effective Hamiltonian is that of a driven nonlinear oscillator with 3/2 degrees of freedom. Generally speaking, the atom-photon interaction in a high-finesse cavity is, mainly, the interaction between internal (electronic) and external (motional) atomic degrees of freedom and the cavity field. A corresponding one-dimensional model, including the interaction of all those degrees of freedom, has been introduced in papers PK01 ; PS01  in the context of Hamiltonian chaos. In this paper, we use a slightly generalized version of that model with three degrees of freedom to study the effects of chaos and Lévy flights in the strongly-coupled atom-field system.

In Sec. II the semiclassical basic equations in the Heisenberg representation are derived in the form of six coupled nonlinear equations with two control parameters, the normalized detuning $`\delta `$ and atomic recoil frequency $`\alpha `$. Taking into account that the frequency of the atomic (and field) Rabi oscillations is much more larger than the frequency of translational oscillations, reduced Bloch-like equations of the atom-field internal motion are derived and solved in Sec. III. The motion of the center of the atom mass is governed by a single equation of a parametrically perturbed pendulum. This dynamics generates a stochastic layer with an exponentially small width. In Sec. IV we simulate the basic set of equations at the fixed value of the normalized recoil frequency $`\alpha =10^3`$ corresponding to a light atom in a microcavity with realistic parameters Y99 ; HLD ; EMK ; MF99 ; PFM . The maximal Lyapunov exponents and Poincaré sections are calculated and it is found that the detuning $`\delta `$ is a crucial parameter in transition to chaos. Statistical properties of the atom-photon interaction are considered in Sec. V. We give an evidence of two types of Lévy flights of an atom, one corresponding to an almost linear dependence of the atomic position $`x`$ on time (superdiffusive and superballistic regime), and another one corresponding to small regular oscillations of the atom in the potential well. The Lévy flights influence strongly such statistical properties of atoms and photons as distribution of Poincaré recurrences and moments of $`x`$. The distribution of recurrences and time-evolution of the moments depend on the value of the detuning $`\delta `$ demonstrating different regimes from an almost normal diffusion to a superdiffusion. In Sec. VI we discuss briefly ways to manipulate the atomic motion by varying control parameters and initial internal atomic states. It is possible, in particular, to cool and trap atoms by adjusting the detuning.

## II Basic equations

The basic model of interaction of radiation with matter describes the energy exchange between a two-level atom and a single mode of the quantized radiation field in an ideal lossless cavity JC . In general, this interaction should involve not only the internal atomic transitions and field states but also the center-of-mass motion of the atom. With the recoil effect to be included into consideration, the standard Jaynes-Cummings Hamiltonian can be extended as follows:

$$H=\frac{\widehat{p}^2}{2m}+\mathrm{}\omega _a\widehat{s}_z+\mathrm{}\omega _f\widehat{a}^+\widehat{a}\mathrm{}\mathrm{\Omega }_0(\widehat{a}^+\widehat{s}_{}+\widehat{a}\widehat{s}_+)\mathrm{cos}k_f\widehat{x},$$
(1)
where $`\widehat{x}`$ and $`\widehat{p}`$ are the atomic position and momentum operators, respectively. Transitions between two electronic states, separated by the energy $`\mathrm{}\omega _a`$, are described by the spin operators with the commutation relations, $`[\widehat{s}_+,\widehat{s}_{}]=2\widehat{s}_z`$ and $`[\widehat{s}_z,\widehat{s}_\pm ]=\pm \widehat{s}_\pm `$. The photon annihilation and creation operators with the commutation rule $`[\widehat{a},\widehat{a}^+]=1`$ characterize a selected mode of the radiation field of the frequency $`\omega _f`$ and the wave number $`k_f`$ in a lossless cavity of the Fabry-Perot type. The parameter $`\mathrm{\Omega }_0`$ is the amplitude value of the atom-field dipole coupling which depends on the position of an atom inside a cavity. As it is usually adopted in cavity QED, we write down the Heisenberg equations for the external atomic operators, $`\widehat{x}`$ and $`\widehat{p}`$, and for slowly-varying amplitudes of the field and spin operators $`\widehat{a}(t)=\widehat{a}\mathrm{exp}(i\omega _ft),\widehat{a}^+(t)=\widehat{a}^+\mathrm{exp}(i\omega _ft),\widehat{s}_\pm (t)=\widehat{s}_\pm \mathrm{exp}(\pm i\omega _ft)`$ and $`\widehat{s}_z(t)=\widehat{s}_z`$

$$\begin{array}{ccc}\frac{d}{dt}\widehat{x}=\frac{\widehat{p}}{m},\hfill & & \\ \frac{d}{dt}\widehat{p}=\mathrm{}k_f\mathrm{\Omega }_0(\widehat{a}^+\widehat{s}_{}+\widehat{a}\widehat{s}_+)\mathrm{sin}k_f\widehat{x},\hfill & & \\ \frac{d}{dt}\widehat{s}_+=i(\omega _f\omega _a)\widehat{s}_++2i\mathrm{\Omega }_0\widehat{a}^+\widehat{s}_z\mathrm{cos}k_f\widehat{x},\hfill & & \\ \frac{d}{dt}\widehat{s}_{}=i(\omega _f\omega _a)\widehat{s}_{}2i\mathrm{\Omega }_0\widehat{a}\widehat{s}_z\mathrm{cos}k_f\widehat{x},\hfill & & \\ \frac{d}{dt}\widehat{s}_z=i\mathrm{\Omega }_0(\widehat{a}^+\widehat{s}_{}\widehat{a}\widehat{s}_+)\mathrm{cos}k_f\widehat{x},\hfill & & \\ \frac{d}{dt}\widehat{a}^+=i\mathrm{\Omega }_0\widehat{s}_+\mathrm{cos}k_f\widehat{x},\hfill & & \\ \frac{d}{dt}\widehat{a}=i\mathrm{\Omega }_0\widehat{s}_{}\mathrm{cos}k_f\widehat{x}.\hfill & & \end{array}$$
(2)
To avoid cumbersome notations in Eqs. (2) we use for the amplitudes the same notations as for the respective whole operators.

In order to derive a tractable closed set of equations for expectation values from the Heisenberg operator equations (2), we use the semiclassical approximation. It means that all the operators and their products in Eqs. (2) are averaged over an initial quantum state, which is supposed to be a product state of the translational, electronic, and field states. The expectation values of all the operator products are factorized to the products of the respective expectation values, e.g., $`(\widehat{a}^+\widehat{s}_{}+\widehat{a}\widehat{s}_+)\mathrm{sin}k_f\widehat{x}`$ $`=(\widehat{a}^+\widehat{s}_{}+\widehat{a}\widehat{s}_+)\mathrm{sin}(k_f\widehat{x})`$. By choosing the following dimensionless expectation values: $`x=k_f\widehat{x}`$, $`p=\widehat{p}/\mathrm{}k_f`$, $`a_x=\frac{1}{2}\widehat{a}+\widehat{a}^+`$, $`a_y=\frac{1}{2i}\widehat{a}\widehat{a}^+`$, $`s_x=\frac{1}{2}\widehat{s}_{}+\widehat{s}_+`$, $`s_y=\frac{1}{2i}\widehat{s}_{}\widehat{s}_+`$, as dynamical variables, we finally get from Eqs. (2) a nonlinear dynamical system

$$\begin{array}{ccc}\dot{x}=\alpha p,\hfill & & \\ \dot{p}=2(a_xs_x+a_ys_y)\mathrm{sin}x,\hfill & & \\ \dot{s}_x=\delta s_y+2a_ys_z\mathrm{cos}x,\hfill & & \\ \dot{s}_y=\delta s_x2a_xs_z\mathrm{cos}x,\hfill & & \\ \dot{a}_x=s_y\mathrm{cos}x,\hfill & & \\ \dot{a}_y=s_x\mathrm{cos}x,\hfill & & \end{array}$$
(3)
that describes the interaction between three degrees of freedom in the strongly coupled atom-field system, translational $`(p,x)`$, electronic $`(s_x,s_y)`$, and the field $`(a_x,a_y)`$ ones. The dot in Eqs. (3) denotes the derivative with respect to the normalized time $`\tau =\mathrm{\Omega }_0t`$. The control parameters are the normalized recoil frequency $`\alpha =\mathrm{}k_f^2/m\mathrm{\Omega }_0`$ and the normalized detuning between the frequencies of the field mode and the electronic transition, $`\delta =(\omega _f\omega _a)/\mathrm{\Omega }_0`$. The system (3) conserves the energy

$$W=\frac{\alpha }{2}p^22(a_xs_x+a_ys_y)\mathrm{cos}x\delta s_z,$$
(4)
and it possesses two additional first integrals

$$s_x^2+s_y^2+s_z^2=S^2,a_x^2+a_y^2+s_z=N.$$
(5)
The first one is simply the conservation of the atomic probability, and the second one is a conserved total number of excitations, which is known to be a constant in the rotating-wave approximation. The equation of motion for the atomic inversion $`s_z`$ is easily derived with the help of the integral $`N`$

$$\dot{s}_z=2(a_xs_ya_ys_x)\mathrm{cos}x.$$
(6)
The semiclassical approximation, we used for (3), means that the atom as a classical particle with external and internal states moves in a self-consistent classical radiation field.

From the dynamical systems point of view equations (3) represent a system with 3 degrees of freedom (one degree of freedom per a canonically conjugate pair of the generalized momentum and coordinate) in 6-dimensional phase space. Indeed, we reinsert Eq.(6) for $`\dot{s}_z`$ into (3). After that the system (3) describes 3 degrees of freedom: the atom external coordinates ($`x,p`$), the atom internal coordinates ($`s_x,s_y`$), and the field coordinates ($`a_x,a_y`$). There are 2 constrains: energy integral $`E=H(p,x;s_x,s_y;a_x,a_y)`$ and the spin integral $`S`$. The number of excitations $`N`$ should be used to determine $`s_z`$ as a function of other variables. That means that full dynamics is defined in 4-dimensional hyperspace and should have domains of chaos due to its nonintegrability. One can say that the location of the domains of chaotic motion, islands of regular dynamics, set of stationary points, and boundaries define a topology of the system’s flow. The topology is two-parametric $`(E,S)`$ and very complicated. Its description needs a separate investigation.

## III Reduced dynamics and the estimation of the stochastic layer width

We can simplify further the basic equations (3) introducing the combined atom-field variables

$$u=2(a_xs_x+a_ys_y),v=2(a_ys_xa_xs_y)$$
(7)
and using the integrals (5). As a result, one arrives at the closed five-dimensional dynamical system

$$\begin{array}{ccc}\dot{x}=\alpha p,\hfill & & \\ \dot{p}=u\mathrm{sin}x,\hfill & & \\ \dot{u}=\delta v,\hfill & & \\ \dot{v}=\delta u+2(2Ns_z3s_z^2+S^2)\mathrm{cos}x,\hfill & & \\ \dot{s}_z=v\mathrm{cos}x,\hfill & & \end{array}$$
(8)
which generalizes the corresponding equations of the paper PS01  (see Eqs. (3) therein which were derived in the limit of large $`N`$). It is obvious from Eqs. (8) that at exact resonance, $`\delta =0`$, the slow translational variables $`x`$ and $`p`$ are separated from the fast atom-field $`u,v`$, and $`s_z`$, and the system (8) becomes integrable. At $`\delta =0`$, the atom moves in a spatially periodic optical potential $`U=u(0)\mathrm{cos}x`$ with $`u(0)=u(\tau =0)`$ = const, and its center-of-mass motion satisfies the pendulum equation $`\ddot{x}+\alpha u(0)\mathrm{sin}x=0`$. It is easy to find that the dynamics of the internal atomic variable $`s_z`$ satisfies the following equation:

$$\dot{s}_z=\pm 2\sqrt{s_z(s_z^2Ns_zS^2)+C}\mathrm{cos}x,$$
(9)
where $`C`$ is an integration constant, and $`x(\tau )`$ is a solution of the pendulum equation mentioned above. The equation (9) can be integrated in terms of elliptic Jacobian functions with a solution which converges to the well-known Jaynes-Cummings semiclassical solution JC  in the limit $`x`$ = const.

Out of resonance, at $`\delta 0`$, the system (8) exhibits chaotic dynamics. In order to clarify the origin of chaos, consider Eqs. (8) in the limit of large number of excitation $`N`$ and large detunings $`\delta `$ comparing to $`S^2`$. Taking into account that the normalized Rabi frequency is of the order of $`\sqrt{N}>1`$ and is much more larger than the frequency of small translational oscillations, $`\sqrt{\alpha }1`$, the equations for the fast atom-field oscillations are reduced to the Bloch-like form

$$\begin{array}{ccc}\dot{u}=\delta v,\hfill & & \\ \dot{v}=\delta u+4Ns_z\mathrm{cos}x,\hfill & & \\ \dot{s}_z=v\mathrm{cos}x,\hfill & & \end{array}$$
(10)
where the function $`\mathrm{cos}x`$ may be considered approximately as a constant $`c`$ over a period of time of many Rabi oscillations. In this approximation, the quantity

$$u^2+v^2+N(2s_z)^2=R$$
(11)
plays a role of the length of Bloch vector, and the general solution of the Bloch-like equations (10) can be found

$$\begin{array}{c}u=u(0)\left[N\left(\frac{2c}{\mathrm{\Omega }_N}\right)^2+\left(\frac{\delta }{\mathrm{\Omega }_N}\right)^2\mathrm{cos}\mathrm{\Omega }_N\tau \right]+\hfill \\ \hfill +\frac{\delta }{\mathrm{\Omega }_N}v(0)\mathrm{sin}\mathrm{\Omega }_N\tau +\frac{4N\delta c}{\mathrm{\Omega }_N^2}s_z(0)(1\mathrm{cos}\mathrm{\Omega }_N\tau ),\end{array}$$

$$\begin{array}{c}v=\frac{\delta }{\mathrm{\Omega }_N}u(0)\mathrm{sin}\mathrm{\Omega }_N\tau +v(0)\mathrm{cos}\mathrm{\Omega }_N\tau +\hfill \\ \hfill +\frac{4Nc}{\mathrm{\Omega }_N}s_z(0)\mathrm{sin}\mathrm{\Omega }_N\tau ,\end{array}$$
(12)
$$\begin{array}{c}s_z=u(0)\frac{c\delta }{\mathrm{\Omega }_N^2}(1\mathrm{cos}\mathrm{\Omega }_N\tau )\frac{c}{\mathrm{\Omega }_N}v(0)\mathrm{sin}\mathrm{\Omega }_N\tau +\hfill \\ \hfill +s_z(0)\left[\left(\frac{\delta }{\mathrm{\Omega }_N}\right)^2+N\left(\frac{2c}{\mathrm{\Omega }_N}\right)^2\mathrm{cos}\mathrm{\Omega }_N\tau \right],\end{array}$$

where the quantity

$$\mathrm{\Omega }_N=\sqrt{\delta ^2+(2c)^2N}$$
(13)
is similar to the Rabi frequency.

Since the function $`\mathrm{cos}x`$ varies in time slowly comparing to the fast oscillating $`u,v`$, and $`z`$, the atom-field variable $`u`$ may be considered approximately as a spatially independent frequency- and amplitude-modulated signal that parametrically excites the translational motion:

$$\ddot{x}+\alpha u(\tau )\mathrm{sin}x=0,$$
(14)
that follows from the first two equations of the system (8). The modulation has especially simple form for initial conditions $`u(0)=v(0)=0`$ and $`s_z(0)=S`$, that corresponds to the atom prepared initially in the upper state while the field may be initially at any state, and $`c=1`$:

$$u(\tau )=\frac{4N\delta S}{\mathrm{\Omega }_N^2}(1\mathrm{cos}\mathrm{\Omega }_N\tau ).$$
(15)
The Eq. (14) is derived from the following effective classical Hamiltonian

$$=\frac{1}{2}\dot{x}^2\omega ^2\mathrm{cos}x+\omega ^2\mathrm{cos}\mathrm{\Omega }_N\tau \mathrm{cos}x=_0+V,$$
(16)
where $`_0`$ is the unperturbed Hamiltonian of a free pendulum with the following frequency of small oscillations:

$$\omega =\frac{2}{\mathrm{\Omega }_N}\sqrt{\alpha N\delta S}.$$
(17)
Rewriting the perturbation $`V=\omega ^2\mathrm{cos}\mathrm{\Omega }_N\tau \mathrm{cos}x`$ in the form

$$V=\frac{\omega ^2}{2}\left[\mathrm{cos}(x+\mathrm{\Omega }_N\tau )+\mathrm{cos}(x\mathrm{\Omega }_N\tau )\right],$$
(18)
one may consider (16) as the Hamiltonian of a particle moving in the field of three plane waves in a frame moving with the phase velocity of the first wave, while the phase velocities of the second and third waves are $`\mathrm{\Omega }_N`$ and $`\mathrm{\Omega }_N`$, respectively.

As it follows from the general theory of perturbed motion of Hamiltonian systems with 3/2 degrees of freedom ZS91 , the Hamiltonian (16) induces chaotic dynamics in the so-called stochastic layer that appears due to the separatrix splitting. Let us consider the motion in the neighborhood of unperturbed separatrix of the pendulum (16). Consider the Poincaré-Melnikov integral

$$\begin{array}{c}\mathrm{\Delta }E=_{\mathrm{}}^{\mathrm{}}\{_0,V\}𝑑\tau =\hfill \\ \hfill =\omega ^2_{\mathrm{}}^{\mathrm{}}\dot{x}\mathrm{sin}(x\mathrm{\Omega }_N\tau )𝑑\tau ,\end{array}$$
(19)
where $`\{_0,V\}`$ is the Poisson bracket. This integral describes changes of the atomic translational energy at the separatrix $`_0=E_s=\omega ^2`$. To estimate (19) for the dynamics near the separatrix, one can use for $`x`$ and $`\dot{x}`$ their known unperturbed separatrix solutions

$$\begin{array}{c}x_s=4\mathrm{arctan}\mathrm{exp}[\pm \omega (\tau \tau _n)],\\ \dot{x}_s=\pm \frac{2\omega }{\mathrm{cosh}[\omega (\tau \tau _n)]},\end{array}$$
(20)
where $`\tau _n`$ is introduced as an initial condition. Using the solutions (20), we get

$$\mathrm{\Delta }E_n=\pm 2\omega ^2_{\mathrm{}}^{\mathrm{}}\frac{dt}{\mathrm{cosh}t}\mathrm{sin}\left(x\frac{\mathrm{\Omega }_N}{\omega }t\varphi _n\right),$$
(21)
where the new time $`t=\omega (\tau \tau _n)`$ and phase $`\varphi _n=\mathrm{\Omega }_N\tau _n`$ were introduced. The integral (21) has been calculated to give

$`\mathrm{\Delta }E_n=\mathrm{\Delta }E_s\mathrm{sin}\varphi _n,`$ (22)
$`\mathrm{\Delta }E_s=2\pi \mathrm{\Omega }_N^2{\displaystyle \frac{\mathrm{exp}(\pi \mathrm{\Omega }_N/2\omega )}{\mathrm{sinh}(\pi \mathrm{\Omega }_N/\omega )}}.`$ (23)
The oscillating function $`\mathrm{\Delta }E_n`$ has simple zeroes that implies transversal intersections of stable and unstable heteroclinic manifolds of saddle points known as a complicated heteroclinic structure. On the basis of general properties of motion near the separatrix, the separatrix map can be introduced ZS91 .

$`E_{n+1}=E_n+\mathrm{\Delta }E_s\mathrm{sin}\varphi _n,`$ (24)
$`\varphi _{n+1}=\varphi _n+{\displaystyle \frac{\mathrm{\Omega }_N}{\omega }}\mathrm{ln}{\displaystyle \frac{32E_s}{|E_{n+1}E_s|}}.`$ (25)
The condition

$$\begin{array}{c}K=\mathrm{max}\left|\frac{\delta \varphi _{n+1}}{\delta \varphi _n}1\right|\hfill \\ \hfill \frac{2\pi \mathrm{\Omega }_N^3}{\omega }\frac{\mathrm{exp}(\pi \mathrm{\Omega }_N/2\omega )|\mathrm{sin}\varphi _n|}{|E_{n+1}E_s|\mathrm{sinh}(\pi \mathrm{\Omega }_N/\omega )}1\end{array}$$
(26)
defines the stochastic layer width which can be estimated in the case of the large parameter $`\mathrm{\Omega }_N/\omega 1`$ (see the respective estimations with real atoms in the concluding section) as follows:

$$\delta E_s|E_{n+1}E_s|\frac{8\pi \mathrm{\Omega }_N^3}{\omega }\mathrm{exp}\left(\frac{\pi \mathrm{\Omega }_N}{2\omega }\right).$$
(27)
The dimensionless width of the stochastic layer is finally given by

$$\frac{\delta E_s}{\omega ^2}8\pi \left(\frac{\mathrm{\Omega }_N}{\omega }\right)^3\mathrm{exp}\left(\frac{\pi \mathrm{\Omega }_N}{2\omega }\right),$$
(28)
where the large parameter

$$\frac{\mathrm{\Omega }_N}{\omega }=\frac{\delta ^2+4N}{2\sqrt{\alpha N\delta S}}$$
(29)
under the conditions $`N,\delta S1`$ and $`\delta ^24N`$ is estimated as

$$\frac{\mathrm{\Omega }_N}{\omega }2\sqrt{\frac{N}{\alpha \delta }}.$$
(30)
For the considered case the width of the stochastic layer of the reduced atom-field dynamics is exponentially small in (28), multiplied by a large parameter. Due to (30) the final width depends on the control parameters $`N`$, $`\alpha `$ and $`\delta `$, and the formula (28) is useful in estimating the ranges of the control parameters where one may expect chaotic motion.

The estimation (28) provides the lower bound for the the width of the stochastic layer that appears due to the simplest harmonic modulation (15). Small changes in energy produce comparatively small changes in frequency of oscillations. Nearby the bottom of potential wells and high over potential hills (where the energy is much less and much greater than $`E_s`$), small changes in frequency give rise to respectively small changes in phase during the period of oscillations. Nearby the unperturbed separatrix, where the period goes to infinity, even small changes in frequency lead to dramatic changes in phase. This is the reason of exponential instability of the parametric oscillator (14) and (15) which models chaotic motion of the atom moving through a periodic standing wave.

## IV Lyapunov exponents and Poincaré sections

In this section, we present numerical simulations with the basic set of equations (3) and the integrals of motion (4) and (5) with $`S^2=3/4`$ (the actual value of $`S^2`$ has no importance since it always can be renormalized to 1) and $`N=10`$. The system (3) has two control parameters, the normalized detuning $`\delta `$ between the atomic transition and cavity frequencies and the normalized recoil frequency $`\alpha `$. As it will be estimated in the concluding section, $`\alpha `$ is in the range from $`10^5`$ to $`10^2`$ for real atoms. We choose $`\alpha =10^3`$ in simulations throughout the paper.

The detuning $`\delta `$, as it was shown in PK01 ; PS01 , is the crucial parameter in transition to chaos in the atom-field system with the center-of-mass motion. It is obvious from the set (8), which is equivalent to the basic one (3), that at exact resonance, $`\delta =0`$, the motion is regular. At large detunings, $`\delta N`$, the motion is expected to be quasiregular since the nonlinear term in the fourth equation of the set (8) is small, compared to the linear term of the same equation. With far-detuned light, one does not expect pronounced atomic Rabi oscillations. In order to find the range of the detunings, where the motion is expected to be chaotic, we compute the dependence of the maximal Lyapunov exponent $`\lambda `$ on the detuning $`\delta `$.

Lyapunov exponents characterize the behavior of close trajectories in phase space. Consider a trajectory, some sequence of time instants $`\tau _0`$, $`\tau _1`$, $`\tau _2`$, … with equal intervals $`\mathrm{\Delta }\tau `$, and a ball around the initial point of the trajectory. Lyapunov numbers $`\mathrm{\Lambda }_j`$ ($`j=1,\mathrm{},M`$; M is a number of variables) show per-interval $`\mathrm{\Delta }\tau `$ changes of the axes of the “ellipsoid” of the deformed ball ASY . In our case $`M=6`$ after exclusion of $`s_z`$. The $`k`$-th Lyapunov exponent is defined as $`\lambda _k=\mathrm{ln}\mathrm{\Lambda }_k`$. Typically $`\mathrm{\Lambda }_k`$ depends on time and $`\lambda _k`$ should be replaced by their mean values ASY . In Hamiltonian systems, due to the phase volume conservation, $`\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{}\mathrm{\Lambda }_M=1`$ and $`\lambda _1+\mathrm{}+\lambda _M=0`$. For integrable system all $`\mathrm{\Lambda }_k`$ are pure imaginary and they make pairs: $`\lambda _2=\lambda _1`$, $`\lambda _4=\lambda _3`$, $`\lambda _6=\lambda _5`$, since the number $`M`$ of equations is even. This result follows from the so-called Liouville-Arnold theorem A88 . In our case, for the reduced system of 6 variables ($`p,x;a_x,a_y;s_x,s_y`$) and two constrains (integrals of motion), we have two imaginary pairs, say $`\lambda _{1,2}=\pm i\sigma _1`$, $`\lambda _{3,4}=\pm i\sigma _2`$ ($`\sigma _{1,2}`$ real), and $`\lambda _5`$, $`\lambda _6`$ that satisfy the condition $`\lambda _5+\lambda _6=0`$, and therefore $`\mathrm{\Lambda }_5\mathrm{\Lambda }_6=1`$, i.e. $`\mathrm{\Lambda }_6=1/\mathrm{\Lambda }_5`$. The chaos means that $`\lambda _{5,6}`$ are real LL ; ASY . If, say $`\mathrm{\Lambda }_5<1`$ ($`\lambda _5<0`$), then $`\mathrm{\Lambda }_6>1`$ ($`\lambda _6>0`$) and $`\lambda _6`$ is called maximal Lyapunov exponent. It has a nice physical meaning: the maximal Lyapunov exponent measures a rate of the separation of initially close trajectories, and typically for the practical goal, the mean value $`\lambda =\overline{\lambda _6}`$ over time is used. To compute $`\lambda `$, we use the standard algorithm LL

$$\lambda =\underset{n\mathrm{}}{lim}\frac{1}{n}\underset{k=1}{\overset{n}{}}\mathrm{ln}\frac{\mathrm{\Delta }(\tau _k)}{\mathrm{\Delta }(\tau _{k1})},$$
(31)
where $`\mathrm{\Delta }(\tau _{k1})`$ is a distance between two close trajectories at time $`\tau _{k1}`$, and the value $`\mathrm{\Delta }(\tau _k)/\mathrm{\Delta }(\tau _{k1})`$ shows the level of separations of the trajectories during the interval $`(\tau _{k1},\tau _k)`$.

In the case that separation doesn’t go exponentially, $`\lambda =0`$. This happens when $`\delta =0`$ since the system becomes integrable, and the exponential separation disappears.

The corresponding results of computing maximal Lyapunov exponent $`\lambda `$ are presented in FIG. 1 with three different initial values of the atomic population inversion $`s_z(0)=0.863`$; 0; 0.863, respectively. The other initial conditions are the following: $`a_x(0)=s_x(0)=0`$, $`x(0)=3.14`$, $`p(0)=2`$, and $`a_y(0)`$ and $`s_y(0)`$ are found from the Eqs.(4) and (5) with given $`S=\sqrt{3}/2`$ and $`N=10`$. The value $`s_z(0)=0.863`$ corresponds to the atom initially prepared closely to its ground state for which $`s_z=\sqrt{3}/2`$. Note that the unusual amplitude values of the atomic population inversion we have are the result of the chosen normalization $`S^2=3/4`$. The atom with $`s_z(0)=0.863`$ is prepared closely to its excited state. In both the cases, the initial components of the transition electric dipole moment, $`s_x(0)`$ and $`s_y(0)`$, are almost zero with $`|s_z(0)|=0.863`$. The atom with $`s_z(0)=0`$ has a maximal electric dipole moment.

As it was expected, at exact resonance $`(\delta =0)`$, the maximal Lyapunov exponent is exactly equal zero in all the cases assuming a regular motion. As it follows from the results of previous section, a stochastic layer appears with infinitesimally small values of detuning, but its width decreases fast with increasing $`\delta `$ (see (28) and (29)). We find $`\lambda 0`$ with $`|\delta |3`$. In physical terms, it means that at exact resonance an atom will periodically exchange excitation with the field, whereas far off resonance its internal states will not (almost) be affected by the field. This interplay results in a maximum of the $`\lambda (\delta )`$ dependence with almost the same maximal values for $`s_z(0)=\pm 0.863`$ and 0. The results, however, are different in the range $`3\delta 3`$ with different initial values of the population inversion $`s_z`$. It is easy to understand why it is. As it follows from the Bloch-like solution (12) for $`s_z`$, the atom starting, say, in its ground state, $`s_z(0)=\sqrt{3}/2`$ and $`u(0)=v(0)=0`$, could reach the upper state $`s_z=\sqrt{3}/2`$ only with $`\delta =0`$. The same is valid with the other initial values of $`s_z(0)`$: the atom starting with $`s_z(0)=\sqrt{3}/2`$ (or with $`s_z(0)=0`$) will not reach $`s_z=\sqrt{3}/2`$ (or $`s_z=\pm \sqrt{3}/2)`$ except for the case of exact resonance, $`\delta =0`$. Thus, the dependencies $`\lambda (\delta )`$ are different in the range $`3\delta 3`$ with different values of $`s_z(0)`$ in spite of the fact that the maximal Lyapunov exponent is computed over the rather long trajectory.

The model Hamiltonian (1) can be easily generalized to an ensemble of indistinguishable two-level atoms. In the semiclassical approximation, we have not observed any pronounced differences in the strength of chaos (that is characterized by the values of $`\lambda `$) with different initial internal atomic states. Interatomic quantum correlations, which occur through the mediation of the field generated by the atomic ensemble, have been shown in PK00  (where the model with hot moving atoms but without recoil has been considered) to play a significant role in the atom-field dynamics. Much more strong chaos has been numerically found PK00  in the vacuum Rabi oscillations with atoms initially prepared in the so-called superfluorescent state (with all the atoms to be uncorrelated initially and occupying their excited states) than with atoms initially prepared in the superradiant state (with initially strongly correlated atoms having a macroscopic electric dipole moment).

We numerically construct single-trajectory Poincaré sections of motion in the system (3) with three degrees of freedom and project them on the plane of the atomic external variables $`(x,p)`$. FIG. 2 presents these sections with the atom initially prepared close to the ground state and with two different values of the atom-field detuning $`\delta =1.2`$ (see (a) and (b) fragments) and $`\delta =1.92`$ (see (c) and (d) fragments). In the latter case chaos is not as strong as in the first case (see FIG. 1a). A fairly regular web structure, that is seen in the fragment (b) computed over a comparatively short integration time, breaks down with increasing integration time (the fragment (a)). For comparison, we present in FIG. 2c the Poincaré section computed under the same conditions as in FIG. 2a but with $`\delta =1.92`$ (an additional trajectory with $`p(0)=0.2`$ is plotted in the fragment (c)). FIG. 2d demonstrates the Poincaré section at $`\delta =1.92`$ with an increased initial momentum $`p(0)=20`$.

## V Statistical properties of the atom-photon interaction

In this section a detailed analysis of the atom and photon chaotic dynamics will be considered. A sensitive control parameter is $`\delta =(\omega _f\omega _a)/\mathrm{\Omega }_0`$, detuning of the field and atom frequencies. For the sake of convenience we specify two values of $`\delta `$: 1.2 and 1.92. It follows from calculating maximal Lyapunov exponents in the previous section that the smaller is $`\delta `$ (in the range $`0.5\delta 3`$), the stronger is mixing and chaos, and one can expect that the case with $`\delta =1.92`$ is, being chaotic (but not with the atom prepared in the upper state or close to it), more intermittent than the case $`\delta =1.2`$. This property of the atom-photon dynamics will be quantitatively characterized below.

The difference of a trajectory projection on the $`xp`$ plane $`(a_x=s_x=0`$) is evident from FIG. 3 where the density modulation has been used: a change of each density appears after $`\mathrm{\Delta }\tau `$ points of the mapping the trajectories ($`\mathrm{\Delta }\tau =75,000`$ for (a) and $`\mathrm{\Delta }\tau =49,000`$ for (b)). The narrow strips of the same density in FIG. 3b indicate a long stay of atom in the corresponding part of the $`xp`$ plane with oscillations in the potential well and a small change of the amplitude of the oscillations. In contrast to this pattern, the distribution of densities in FIG. 3a is more uniform manifesting much better mixing, although some traces of the intermittency persist.

The difference between $`\delta =1.2`$ and $`\delta =1.92`$ is also evident from FIG. 4 where a dependence $`x=x(\tau )`$ is shown. The intermittent case ($`\delta =1.92`$) has very long “flights” known also as Lévy flights SZK ; Z98 . There are two types of flights in FIG. 4b. One category of flights corresponds to the almost linear dependence of $`x=x(\tau )`$, while the other one corresponds to the stagnation of the trajectory near some value of $`x`$. FIG. 4c shows the flight in the $`xp`$ plane where the ballistic dynamics coexists or alternate stagnations. Both categories of flights are well understood from FIG. 2c,d: ballistic dynamics along $`x`$ in FIG. 2c is responsible for the linear dependence of $`x=x(\tau )`$, while the trajectory can stay very long near the saddle points as in FIG. 2d (the dark area near a saddle point). The case of $`\delta =1.2`$ in FIG. 2a,b is very different and flights of both categories are rare, if ever.

Just the presence of flights and intermittent behavior of the physical variables strong fluences the statistical properties of atoms and photons. We will use two important characteristics of the atom-photon variables: distribution of Poincaré recurrences and moments of the atom coordinate $`x`$. Consider a small phase volume $`\mathrm{\Delta }\mathrm{\Gamma }`$ and $`P(\mathrm{\Delta }\mathrm{\Gamma };\tau )`$ as a probability density of a trajectory to return first time back to $`\mathrm{\Delta }\mathrm{\Gamma }`$ at time instant $`\tau (\tau +d\tau )`$ if initially started at $`\mathrm{\Delta }\mathrm{\Gamma }`$ at $`\tau =0`$. Then the density probability to return first time to $`\mathrm{\Delta }\mathrm{\Gamma }`$ is

$$P(\tau )=\underset{\mathrm{\Delta }\mathrm{\Gamma }0}{lim}\frac{1}{\mathrm{\Delta }\mathrm{\Gamma }}P(\mathrm{\Delta }\mathrm{\Gamma };\tau )$$
(32)
with a normalization condition

$$_0^{\mathrm{}}P(\tau )𝑑\tau =\underset{\mathrm{\Delta }\mathrm{\Gamma }0}{lim}\frac{1}{\mathrm{\Delta }\mathrm{\Gamma }}_0^{\mathrm{}}P(\mathrm{\Delta }\mathrm{\Gamma };\tau )𝑑\tau =1.$$
(33)
The probability $`P(\tau )`$ does not depend on the choice of $`\mathrm{\Delta }\mathrm{\Gamma }`$ and for “good” chaotic mixing decays exponentially ZEN

$$P(\tau )=(1/h)e^{h\tau }$$
(34)
with the mean recurrence time

$$\tau _{\mathrm{rec}}=1/h=_0^{\mathrm{}}\tau P(\tau )𝑑\tau ,$$
(35)
and $`h`$ as Kolmogorov-Sinai entropy.

The general situation is more complicated since an algebraic behavior

$$P(\tau )1/\tau ^\gamma ,\tau \mathrm{},$$
(36)
is possible for large $`\tau `$ due to intermittent chaos. For bounded Hamiltonian dynamics, $`\tau _{\mathrm{rec}}<\mathrm{}`$ (Kac lemma), and the condition $`\gamma >2`$ should exist. Nevertheless, strongly intermittent dynamics sometimes does not permit us to achieve the limit at $`\tau \mathrm{}`$ and many different intermediate asymptotics can appear. FIG. 5 shows the distribution of recurrences that is close to the exponential one as in (34) for $`\delta =1.2`$ and to the algebraic one as in (36) with $`\gamma 2`$ $`(\tau >10^5`$) for $`\delta =1.92.`$

The difference due to intermittency also occurs for the moments

$$x^{2m}\tau ^{\mu (m)},$$
(37)
where the so-called transport exponent $`\mu (m)`$ varies for different time windows. The behavior of $`x^{2m}`$ is shown in FIG. 5. For $`\delta =1.2`$, the value $`\mu (1)`$ is close to 2 for $`\tau <10^3`$ and corresponds to the ballistic dynamics. For $`\tau >10^3`$ and $`\delta =1.2`$, we have $`\mu (1)1.13`$ that corresponds to a weak superdiffusion that is fairly close to the normal diffusion with $`\mu (1)=1`$. A very different behavior for moments appears for $`\delta =1.92`$ where there are many long-lasting flights. For $`\tau >10^3`$, $`\mu (1)2.2`$ that corresponds to a superballistic transport with an acceleration. This behavior can be explained as a result of long flights when atoms move in the photon’s field acquiring acceleration. This type of transport is self-similar and $`\mu (4)8.8=4\mu (1)`$.

## VI Manipulation of Atoms

In this section we would like to make a few comments related to the manipulation of atoms by changing different control parameters. As it was shown in Sec. V, a change of $`\delta `$ leads to a possibility of a sensitive control of the Lévy flights and, as a result, to cool the atoms which have the lower chaotic dispersion the longer the flight is. At the same time, simulations show fast mixing on the $`a_xa_y`$ and $`s_xs_y`$ planes. More precisely, spectral properties of the atomic dynamics are sensitively controlled by the parameter $`\delta `$. Let us demonstrate it using a simplified analysis.

Consider $`x=x(\tau )`$ as the only variable that describes the dynamics or the most essential part of the dynamics, and introduce a generation function

$$G(x,\tau ;\nu )=e^{i\nu [x(\tau )x_0]},x_0x(0).$$
(38)
Then

$$I(\tau ;x_0)_{\mathrm{}}^{\mathrm{}}𝑑\nu G(x,\tau ;\nu )=2\pi \delta [x(\tau )x_0].$$
(39)
The expression $`I(\tau ;x_0)`$ can be “coarse-grained” over $`x_0`$, i.e.

$$\begin{array}{c}I(\tau ;x_0)_{x_0}\frac{1}{\mathrm{\Delta }x_0}𝑑x_0I(\tau ,x_0)=\hfill \\ \hfill =\frac{2\pi }{\mathrm{\Delta }x_0}𝑑x_0\delta [x(\tau )x_0].\end{array}$$
(40)
The presence of $`\delta `$-function indicates recurrences to $`x_0`$ within an interval $`\mathrm{\Delta }x_0`$ at time instant $`\tau `$ within an interval $`\mathrm{\Delta }\tau _0`$. For $`\tau \mathrm{}`$ we can neglect $`n`$-triple recurrences from $`n2`$ and leave only the first recurrences. Then

$$\begin{array}{c}I(\tau ;x_0)_{x_0}=P(\tau )=\hfill \\ \hfill =\underset{x_0\mathrm{\Delta }x_0}{}\mathrm{const}/\left|\frac{dx(\tau )}{dx_0}\right|_{x(\tau )=x_0}.\end{array}$$
(41)
The expression (41) shows that for “good” chaotic systems

$$|dx(\tau )/dx_0|\mathrm{exp}(h\tau ),$$
(42)
and we arrive at (34). For the intermittent dynamics, the sum in (41) consists of two types of terms, the same as (42) and the algebraic growth

$$|dx/dx_0|\tau ^\gamma ,$$
(43)
with the value of $`\gamma `$ depending on the type of intermittency. For a fairly large $`\tau `$, the term (43) survives and we arrive at the case (36).

From another side,

$$\frac{^{2m}}{\nu ^{2m}}G(x,\tau ;\nu _{\nu =0}=(1)^m|x(\tau )x_0|^{2m},$$
(44)
and we obtain the moments of $`x(\tau )`$. This shows that the moments and their spectral properties are coupled to the recurrences distribution through the generating function $`G(x,\tau ;\nu )`$ which one would expect to obtain from experiments. When the moments are infinite, the expression (44) can be replaced by the following one

$$\frac{^\beta }{\nu ^\beta }G(x^\beta ,\tau ;\nu ^\beta _{\nu =0}=\mathrm{const}|x(\tau )x_0|^\beta ,$$
(45)
with an appropriate value of $`\beta `$ (see more discussion in SZ97 ).

The main way of controlling the properties of $`G(x,\tau ;\nu )`$ is to change the system’s topology in phase space. Speaking about the topology, we have in mind the phase pattern (see also the end of the Sec. II) that includes the singular points, curves, and partitioning of the domains of chaos and islands. To illustrate how the system is sensitive to small variations of the initial conditions that change the full energy, we show in FIGs. 6a and b the Poincaré sections with $`s_z(0)=0.863`$ and $`s_z(0)=0.8660254`$, respectively, at $`\delta =0.4`$ and under the other equal conditions. Very small difference in the values of the initial tipping angle between the direction of the Bloch vector and the axis $`z`$ gives rise to cardinally different motion with $`p(0)=2`$, chaotic oscillations in the wide range of the atomic momenta with $`s_z(0)=0.863`$ and small regular translational oscillations nearly the bottom of a potential well with $`s_z(0)=0.8660254`$. Transition from order to chaos takes place with the latter value of the initial atomic population inversion only with much more large values of the initial momentum, $`p(0)40`$.

The original system (3) has three degrees of freedom and it is not studied yet. Nevertheless, we were able to demonstrate by simulation a bifurcation of a hyperbolic point into the elliptic one, although we are not able to provide an analytical description at the moment since the FIG. 7 is just a projection of a trajectory in the 4-dimensional hyperspace onto the plane $`(p,x)`$. By a change of $`\delta `$ near $`\delta ^{}0.222`$, the saddle on the $`xp`$ plane with $`a_x=s_x=0`$ transforms into the elliptic point. The trapping potential well of the finite size on $`xp`$ plane occurs as a result of the bifurcation: this bifurcation will be studied in detail in another paper.

## VII Conclusion

A system with one or more cold atoms strongly coupled to a single mode of the cavity field is ideal for testing fundamentals of quantum mechanics and its corresponding to classical mechanics. Based on our understanding of the nonlinear dynamics of the atom-photon interaction in a standing-wave high-finesse cavity, new ways to manipulate and control atomic motion can be opened. We have shown that the motion is very sensitive to the atom-cavity detuning $`\delta `$. Varying $`\delta `$, one can design topology of the underlying phase space creating zones of trapping, quasi-trapping or acceleration, quasi-regular and stochastic webs, etc. It may provide new schemes for cooling, trapping, and accelerating atoms.

To give an idea about the values of the magnitudes we have used in numerical simulations, we need to estimate the range of values of the normalized recoil frequency $`\alpha =\mathrm{}k_f^2/m\mathrm{\Omega }_0`$ with real atoms and cavities. We will use the parameters of the real experiments with single atoms in the strong-coupling regime Y99 ; MF99 , for which the maximal atom-field coupling strength $`\mathrm{\Omega }_0`$ exceeds the decay rates of the cavity field and of the atomic dipole. Atoms were collected in a magneto-optical trap and cooled down to microKelvin temperatures, before entering a microscopic high-finesse Fabry-Perot cavity with $`Q10^6`$, $`\mathrm{\Omega }_02\pi (10^710^8)`$ Hz and $`k_f2\pi 10^6m^1`$. With these values of the parameters, one can estimate $`\alpha `$ to be in the range $`10^510^2`$ depending on the atomic mass and $`\mathrm{\Omega }_0`$.

## VIII Acknowledgment

We thank L.E. Kon’kov for his help in preparing FIG. 1. This work was supported by the Russian Foundation for Basic Research under Grant No. 02-02-17796; the Office of Naval Research, Grants No. N00014-96-1-0055 and No. N00014-97-1-0426, and U.S. Department of Energy Grant No. DE-FG02-92ER54184. The simulations were supported by the National Science Foundation cooperative agreement No. AC1-9619020 through computing resources provided by the National Partnership for Advanced Computational Infrastructure at the San Diego Supercomputer Center. S.P. thanks the Courant Institute of Mathematical Sciences, New York University, for the hospitality during his stay.


text:
# The inner structure of Zeldovich’ pancakes

## 1 Introduction

Structure formation in the Universe is a rich and challenging problem, which touches many sides of Physics. It is generally believed that the presently observed large scale structures have been generated by the process of the gravitational instability, acting on initially small density perturbations. The temperature fluctuations on the uniform Cosmic Microwave Background (CMB) radiation provide an image of these fluctuations, which were the seeds of the present large scale distribution of matter . The observed CMB angular power spectrum is dominated by a peak at $`1`$ degree of arc, and also shows structures at smaller scales.
According to present consensus, the constituents of our Universe are ordinary matter (5%), dark matter (25%) and dark energy (70%). The three forms can be summarily classified by how their energy densities change with a cosmic scale factor $`a`$: ordinary and dark matter behave as $`a^3`$, radiation, which interacts with ordinary matter, behaves as $`a^4`$, while dark energy, at least in the simplest models thereof, is independent of $`a`$. Hence, even if dark energy dominates today, dark matter and ordinary matter dominated earlier. Further, the Universe is flat, in agreement with an early epoch of inflation. Within that scenario, the primordial density perturbations were random Gaussian with power spectrum $`P(k)Ak`$, which is also consistent with the data. The ordinary matter content can be deduced from observations at e.g. visible or radio frequencies. The presence of dark matter can be inferred from the observed dynamics of cosmic objects, particularly from fast rotation of hydrogen clouds far outside the luminous disc of spiral galaxies, as well as high-velocity dispersion of galaxies in clusters . Dark energy, or quintessence , is equivalent to a non-zero cosmological constant, $`\mathrm{\Lambda }`$, in Einstein’s equations , and there is recent support for a non-zero $`\mathrm{\Lambda }`$ also from redshift observations. All of these direct measurements can be compared with theoretical cosmology and the observed angular structure of the CMB. A sequence of peaks should indeed arise from coherent acoustic oscillations in the baryon-photon fluid during an early epoch. Their amplitudes and relative positions provide another series of tests of cosmological models, and put a different series of constraints on the parameters of such models, c.f.  for a recent review.
The growth of small initial perturbations to large-scale structures is essentially a problem in classical gravitational physics, where all the above mainly enters in the initial conditions, and the how the Universe expands as a whole. The problem is nevertheless far from trivial. A principal difficulty is that a self-gravitational medium has no ground-state, and the dynamics is not ergodic. Contrary to e.g. molecular dynamics, one cannot appeal to a shadowing lemma, and argue that even if a simulation only solves the equations of motion in a rough approximation, the simulated system nevertheless stays close to the real system, with different initial conditions. The single most influential result in structure formation, Zeldovich’ pancake theory, was indeed developed without any recourse to simulation at all, but purely by extending linear analysis into the non-linear domain. The fact that this model and the adhesion model are still the benchmarks today shows the difficulty of firmly establishing the detailed characteristics of the mass agglomeration process.
The goal of this paper is to investigate the non-linear regime of the gravitational instability in the special case that the initial perturbations are planar. We also assume an Einstein-de Sitter Universe, and hence disregard the dark energy component. With these limitations, one can introduce a Lagrangian integration scheme, called the Quintic (Q) model , which is both fast, and exact from collision to collision. It is worth stressing that the Q model is just one of the possible representation of the dynamics . A conceptually equivalent formulation was previously proposed by . These two approaches share the characteristic that they allow to avoid sources of additional errors, like truncation, involved in any approximate numerical scheme (see ).

The focus is here on the inner structure of a pancake. The novelty with respect previous investigations  is that, by making use of a discrete map to approximate the exact Q dynamics, we derive a close analytical expression for the density profile. These results are validated by direct numerical simulation, and compared with analogous results for the adhesion model. We further investigate the pair-wise motion of two pancakes under their mutual attraction.
The paper is organized as follows: in Section 2 we present the general background, and in Section 3 we re-derive, for completeness, the Q model. In Section 4 the evolution of a two particle system is considered: an approximate discrete system is derived and validated by numerical investigations. The analysis is extended in Section 5 where the inner structure of a single cluster is investigated. Improving previous results in  we establish the density profile of a collapsing cluster. In Section 6 we study the dynamics of two massive clusters, and in Section 7 we sum up and discuss our results.

## 2 From Vlasov-Poisson equations to the adhesion model

The evolution of collision-less dark matter in a three dimensional expanding Universe is described by the kinetic Vlasov-Poisson equations. Assume an inertial reference frame, and label the position with $`𝐫`$. Then it is customary to introduce the co-moving coordinate $`𝐱`$ by the following transformation :

$$𝐫=a(t)𝐱,$$
(1)
where the scale factor $`a(t)`$ is function of proper world time. For an Einstein-de Sitter Universe

$$a=\left(\frac{t}{t_0}\right)^{2/3},$$
(2)
where $`t_0^2=6\pi G\rho (t_0)`$, $`\rho (t_0)`$ is the homogeneous density at time $`t_0`$  and $`G`$ is the gravitational constant. The Vlasov-Poisson equations then read:

$$\{\begin{array}{c}_tf+\frac{𝐩}{ma^2}_𝐱f_𝐱\psi _𝐩f=0\hfill \\ \\ ^2\psi =4\pi Ga^2\left(\rho \rho _b\right),\hfill \end{array}$$
(3)
where $`𝐩`$ is the variable conjugated to $`𝐱`$; $`f(𝐱,𝐩,t)`$ is the distribution function in the six-dimensional phase space $`(𝐱,𝐩)`$; $`\psi `$ is the gravitational potential and $`\rho _b`$ is the mean mass density. The particle density $`\rho (𝐱,t)`$ and velocities $`𝐮(𝐱,t)`$ are given, in term of $`f(𝐱,𝐩,t)`$, as:

$$\rho (𝐱,t)=\frac{m}{a^3}f(𝐱,𝐩,t)𝑑𝐩,$$
(4)
$$\rho (𝐱,t)𝐮(𝐱,t)=\frac{1}{a^4}𝐩f(𝐱,𝐩,t)𝑑𝐩.$$
(5)
It is well known, see e.g. , that (3) admit special solutions of the form

$$f(𝐱,𝐩,t)=\frac{a^3\rho (𝐱,t)}{m}\delta ^d(𝐩ma𝐮(𝐱,t)),$$
(6)
where $`d`$ is the dimension of space and $`\delta ^d(.)`$ the $`d`$-dimensional delta function. We will refer to this class as to single-speed solutions, because to each given ($`𝐱,t`$) corresponds a well defined velocity $`𝐮`$. Assuming (6), a closed system on the hydrodynamic level can be derived using (4) and (5):

$$\{\begin{array}{c}_t\rho +3\frac{\dot{a}}{a}\rho +\frac{1}{a}(\rho 𝐮)=0\hfill \\ \\ _t𝐮+\frac{\dot{a}}{a}𝐮+\frac{1}{a}(𝐮)𝐮=𝐠\hfill \\ \\ 𝐠=4\pi Ga(\rho \rho _b),\hfill \end{array}$$
(7)
where we have introduced $`𝐠=\psi /a`$ such that $`\times 𝐠=0`$. It should be stressed that system (7) is valid as long as the distribution function $`f(𝐱,𝐩,t)`$ is in the form (6). Beyond the time of caustic formation, when fast particles cross slow ones, the solution becomes multi-stream. Hence, the pressure-less and dissipation-less hydrodynamical equations of (7) are incomplete.
An ansatz that permits further progress is the so-called condition of parallelism which requires that the peculiar velocity is a potential field, that remains parallel to the gravitational peculiar acceleration field :

$$𝐠=F(t)𝐮$$
(8)
The important assumption, in the present discussion, is that the proportionality $`F(t)`$ is space-independent. That can only strictly be true in the linear regime, where

$$F(t)=4\pi G\rho _bb/\dot{b},$$
(9)
and $`b`$ is the amplitude of the perturbation. A natural further assumption is that only the linearly growing mode is excited , which gives $`bt^{\frac{2}{3}}`$ in our case. Defining the new velocity field $`𝐯=𝐮/(a\dot{b})`$, system (7) reduces to free motion in Eulerian coordinates

$$_b𝐯+(𝐯)𝐯=0,$$
(10)
This is Zeldovich’ pancake model in the original formulation . It is also to referred as the Zeldovich approximation, because the contribution of the decaying mode of the density field has been neglected. The ansatz that the condition of parallelism holds also in the multi-stream region could be approximatively true, but is still a bold assumption to make.
After caustic formation one may think that the resulting change in the gravitational force could be modeled by an effective diffusive term. In the adhesion model, , one hence introduces a term of the form $`\nu ^2𝐯`$ in the right hand side of the equation (10):

$$\{\begin{array}{c}_b𝐯+(𝐯)𝐯=\nu ^2𝐯\hfill \\ \\ 𝐯=\stackrel{~}{\psi }\hfill \\ \\ _b\rho +(\rho 𝐯)=0,\hfill \end{array}$$
(11)
where $`\stackrel{~}{\psi }=\psi /(\dot{b}F(t))`$. In order for the diffusion term to have a smoothing effect only in those regions where the particles crossing takes place, the phenomenological viscosity parameter $`\nu `$ should be small.
Although numerical experiments suggest qualitative agreement, the exact relationship between (3) and (11) is still an open problem. In a recent paper  a derivation of (11) was outlined in one dimension. Comparisons with numerical results suggested however that at least the diffusion term should be modified, to give instead transport equation in the class recently studied by Buchert and co-workers . We will here study the system (3) directly, and establish the inner structure of a mass agglomeration without the ansatz of parallelism.

## 3 The Q model, and its relation to the Zeldovich approximation and the adhesion model

We will consider the evolution of the one-dimensional perturbation in the expanding three dimension Universe, using the discrete approximation of initial condition. Let us assume, that in some interval $`L`$ along the axis $`x`$ all the matter is concentrated in the $`N`$ sheets, which we will from now on call particles.
The Newtonian equations of motion for such $`N`$ particles, interacting via gravity, follow from the Lagrangian :

$$=\underset{i}{}\frac{1}{2}m_i\dot{r_i}^2m_i\varphi (r_i,t),$$
(12)
where $`_r^2\varphi =4\pi G\rho `$. In the point particle picture the density profile reads:

$$\rho (x_i,t)=\underset{x_j}{}m_ja^3\delta (x_ix_j),$$
(13)
where $`x_i`$ is the co-moving coordinate of the $`i`$’th particle, in the direction of which the density and velocities vary.
Expressing (12) as a function of the proper coordinate, $`x_i`$, and assuming (13), the equation of motion of the $`i`$’th particle reads:

$$\frac{d^2x_i}{dt^2}+2\frac{\dot{a}}{a}\frac{dx_i}{dt}4\pi G\rho _b(t)x_i=a^3E_{grav}(x_i,t),$$
(14)
where $`\rho _b(t)`$ is the mean mass density at time $`t`$ and

$$E_{grav}(x_i,t)=2\pi G\underset{j}{}m_j\text{sign}(x_ix_j).$$
(15)
In an Einstein-de Sitter Universe we can make a nonlinear change of variable, $`\tau =t_0\mathrm{log}t/t_0`$, such that (14) takes the form:

$$\frac{d^2x_i}{d\tau ^2}+\frac{1}{3t_0}\frac{dx_i}{d\tau }\frac{2}{3t_0^2}x_i=E_{grav}(x_i,\tau )\text{Q model},$$
(16)
where $`t_0^2=6\pi G\rho (t_0)`$. This is the representation of the Vlasov-Poisson equations which we call the Quintic (Q) model.
The interest of this formulation is that, as for the classical static self-gravitating systems in one dimension, $`E_{grav}`$ is a Lagrangian invariant, proportional to the net mass difference to the right and to the left of a given particle, at a given time. Hence, in between collisions, (16) has an explicit solution :

$$x_i(\tau )=c_1^i\mathrm{exp}(\frac{2(\tau \tau ^n)}{3t_0})+c_2^i\mathrm{exp}(\frac{(\tau \tau ^n)}{t_0})+K_i^n,$$
(17)
where $`K_i^n=(3t_0^2/2)E_{grav}(x_i,\tau )`$, is constant. The coefficients $`c_1^i`$ and $`c_2^i`$ are determined by $`x_i^n=x_i(\tau ^n)`$ and $`w_i^n=\dot{x_i}(\tau ^n)`$, i.e. by the states of the particle at the time of the last crossing, and read:

$$\{\begin{array}{c}c_1^i=\frac{3}{5}\left[x_i^n+t_0w_i^nK_i\right]\hfill \\ \\ c_2^i=\frac{2}{5}\left[x_i^n\frac{3}{2}t_0w_i^nK_i\right].\hfill \end{array}$$
(18)
The form of equation (17) suggests introducing an auxiliary variable $`z=\mathrm{exp}((\tau \tau ^n)/3t_0)`$. The crossing times between neighboring particles (i.e. $`i`$,$`i+1`$) can, hence, be computed by solving numerically the following quintic equation:

$$f(z)=A_{i,i+1}^nz^5B_{i,i+1}^nz^3+C_{i,i+1}^n=0,$$
(19)
where:

$$\{\begin{array}{c}A_{i,i+1}^n=\frac{3}{5}\left[\mathrm{\Delta }x_i^n+t_0\mathrm{\Delta }w_i^n(K_{i+1}^nK_i^n)\right]\hfill \\ \\ \\ B_{i,i+1}^n=(K_{i+1}^nK_i^n)=const\hfill \\ \\ C_{i,i+1}^n=\frac{2}{5}\left[\mathrm{\Delta }x_i^n\frac{3}{2}t_0\mathrm{\Delta }w_i^n(K_{i+1}^nK_i^n)\right],\hfill \end{array}$$
(20)
and $`\mathrm{\Delta }x_i^n=x_{i+1}^nx_i^n`$, $`\mathrm{\Delta }w_i^n=w_{i+1}^nw_i^n`$. Thus, the evolution of the system is recovered by using a version of the event-driven scheme discussed in  and . Details of the numerical implementation are given in .
Let as assume that in the interval $`L`$ all the particles have equal mass $`m_j=L\rho _0/N`$. Then, the gravitational force in equation (16) transforms in

$$E_{grav}(x_i,\tau )=\frac{L}{3Nt_0^2}\underset{j}{}\text{sign}(x_ix_j)=\frac{L}{3Nt_0^2}N_i.$$
(21)
where $`N_i=N_{i,right}N_{i,left}`$ is the difference between the number of particles to the right and to the left of a given particle. At the initial time $`\tau =\tau ^0=0`$, the positions of all the particles are ordered ($`x_{i+1}>x_i`$), and we have $`N_i=N(2i1)`$. With $`K_i^0`$ in the solution (17) constant we get

$$K_i^0=\frac{L}{2N}N_i^0=\frac{L}{2}+\frac{L}{N}(i\frac{1}{2}).$$
(22)
If we now pick the specific initial conditions such that all the $`c_2^i`$’s in (17) are zero, we recover, in the discrete setting, the special case of correlated velocity and density perturbations, such that only the linearly increasing mode is excited. Note that this separation is valid until the first particle crossing, i.e. well into the non-linear regime. In this time interval we have

$$x_i(\tau )=\frac{3}{2}w_i^0t_0\mathrm{exp}(\frac{2\tau }{3t_0})\frac{L}{2}+\frac{L}{N}(i\frac{1}{2}),$$
(23)
Hence, defining the new “time” $`b=b_0\mathrm{exp}(2\tau /3t_0)`$ (with density dimension), the solution of (23) reduces to free motion of the particles. In other words, we have re-derived Zeldovich’ approximation in the framework of the Q model. In the Zeldovich approximation, this solution is assumed to hold for all time, even when the solution becomes multi-stream, while in the adhesion model we have merging of the particles after their crossing. The real system must necessarily lie somewhere in between.

## 4 A two particles system: an iterative map

Consider two particles of equal mass $`m`$ interacting by the Q dynamics, confined in a box of size $`L`$. We choose units of length and mass such that $`L=1`$ and $`m=1/2`$. To streamline the following, we first first introduce the dimensionless variables:

$`\theta `$ $`=`$ $`{\displaystyle \frac{\tau }{t_0}},`$ (24)
$`q_i(\theta )`$ $`=`$ $`{\displaystyle \frac{x_i}{L}},`$ (25)
$`\beta _i(\theta )`$ $`=`$ $`{\displaystyle \frac{6\dot{x}_it_0}{L}}.`$ (26)
Thus, before the time of first crossing, equation (17) takes the form:

$$q_i(\theta )=c_1^i\mathrm{exp}(\frac{2\theta }{3})+c_2^i\mathrm{exp}(\theta )\frac{1}{4}.$$
(27)
where $`c_1^i`$ and $`c_2^i`$ are given by (18). We select a special class of initial condition: the velocities of the particles are arbitrarily assigned, while the positions are determined by setting $`c_2^i`$ to zero (23). This simply means that there is no decaying mode initially, as discussed in the end of the previous section, and we can therefore immediately integrate the equations of motion of every particle up to its first crossing. In addition we chose a system of reference such that the center of mass is at rest at the origin. That means that the initial velocities are $`w_1^0=w_2^0>0`$. For simplicity we will suppress in the following the index $`i`$ and consider the evolution of the leftmost of the two particles.
The time of first crossing, $`\theta _{cross}`$, is then deduced from equation (27) and reads

$$\theta _{cross}=\frac{3}{2}\mathrm{ln}\beta (0).$$
(28)
The rescaled particle’s velocity is $`\beta _{cross,0}=\beta (\theta _{cross})=1`$. Remark that $`\theta _{cross}`$ is solely determined by the initial particle velocity, while $`\beta _{cross,0}`$ is a dimension-less constant, not dependent on the states of the particles at $`\theta =0`$. This means that while the time at which the particles cross ($`\theta _{cross}`$) depends on the initial velocity, their further evolution does not. On the other hand, after crossing the initial balance between positions and velocities no longer holds, and the condition $`c_2^i=0`$ is no longer true. From that time on both the growing and decaying modes in (27) have to be considered in the analysis.
For times larger than $`\theta _{cross}`$, the separation between the particles increases. Conversely, as an effect of the gravitational attraction, their velocities are progressively reduced. Hence, a turning point is reached where the elongation is maximal and the velocities are zero. The inversion time, $`\theta _{turn}`$, follows by differentiating equation (17) with respect to time, and imposing the condition $`\dot{q}(\theta _{turn})=0`$. Thus:

$$\theta _{turn}=\theta _{cross}+\stackrel{~}{\theta }$$
(29)
where $`\stackrel{~}{\theta }=\frac{3}{5}\mathrm{ln}6`$. The corresponding maximum elongation, $`q_{max,0}=q(\theta _{turn})`$, reads:

$$q_{max,0}=\frac{1}{4}\frac{1}{5}6^{\frac{3}{5}}\frac{1}{20}6^{\frac{2}{5}}=\alpha 0.08.$$
(30)
It does not depend on the initial velocity.
Note that the particle inversion occurs before the initial position (i.e. $`|q_i(0)|=\frac{1}{4}(1\beta (0)`$) is recovered: this is the effect of the friction-like term in equation (16), which is responsible for the progressive localization of the system around the center of mass. An analogous damping is shown to occur for the particles velocities. In fact, for $`\theta >\theta _{turn}`$, the particles move inward approaching the time of successive crossing. Due to the inner symmetry, the latter will take place in the origin. The particle rescaled velocity corresponding to the first return reads:

$$\beta _{ret}=\beta _{cross,1}=0.83\beta _{cross,0}=0.83.$$
(31)
By a simple generalization of the previous discussion, it can be shown that, for larger times, the particle executes damped oscillations. The particles move back and forth, displaying a progressive reduction of their mutual distance: the crossings become more frequent as time is increased. Take an index $`n`$ to label particle crossings in the system. The variable $`\theta _n`$ represents the time elapsed between the $`n1`$’th and $`n`$’th crossings, and $`\beta _{cross,n}`$ is the dimensionless velocity of the particle at the $`n`$’th crossing. Then we make the further assumptions of smallness:

$`|\beta _{cross,n}|1.`$ (32)
For large enough $`n`$ these are reasonable approximation. Further, we assume (32) to hold even for an early stage of the evolution: we will return on this point, and provide numerical support to the validity of this ansatz. Within these approximation the dynamics of (17) is reduced to the following simple map:

$$\begin{array}{ccc}\beta _{cross,n+1}\hfill & =& \beta _{cross,n}\left(1\frac{4}{3}\beta _{cross,n}\right)\hfill \\ \theta _{n+1}\hfill & =& 2|\beta _{cross,n}|\hfill \\ \mathrm{\Delta }q_{max,n+1}\hfill & =& \frac{1}{12}\beta _{cross,n}^2\hfill \end{array}$$
(33)
where $`\mathrm{\Delta }q_{max,n}`$ stands for the maximal elongation reached by the particles after $`n`$ encounters. It is interesting that system (33) is actually equivalent to the motion of an impurity advected by an unsteady Burgers flow with a particular numerical relation between frictional and inertial forces, see .
Consider then the finite difference formula:

$$\frac{\mathrm{\Delta }|\beta _{cross}|}{\mathrm{\Delta }\theta }=\frac{|\beta _{cross,n+1}||\beta _{cross,n}|}{\theta _{n+1}},$$
(34)
and rewrite the right hand-side using the first two relations of (33). Approximating the finite differences with differentials and integrating, using as initial condition for $`|\beta _{cross}(\theta =\theta _{cross})|`$ the velocity of at first crossing $`|\beta _{cross,0}|`$ and for $`|\mathrm{\Delta }q_{max}(\theta =\theta _{turn})|`$ the first maximum elongation (30), one ends up with:

$$|\beta _{cross}(\theta )|=|\beta _{cross,0}|\mathrm{exp}\left(\frac{1}{9}(\theta \theta _{cross})\right),\beta _{cross,0}=1,$$
(35)
$$|\mathrm{\Delta }q_{max,0}(\theta )|=2q_{max}\mathrm{exp}\left(\frac{2}{9}(\theta \theta _{turn})\right),q_{max,0}=\alpha .$$
(36)
Here $`|\beta _{cross}(\theta )|`$ is the dimension-less velocity of the particles at the moment of their crossing and $`\mathrm{\Delta }q_{max}(\theta )`$ the maximum separation between the two particles, at rescaled time $`\theta `$. Both the velocity and the inter-particle distance (and also the time elapsed between two successive crossing $`\theta _n(\theta )|\beta _{cross}(\theta )|)`$ decrease exponentially, as function of time $`\theta `$, or, equivalently, a power-law decays is displayed in term of cosmological time $`t`$: $`|\beta _{cross}(t)|t^{1/9}`$ and $`|\mathrm{\Delta }q_{max}(t)|t^{2/9}`$.
Let us turn to numerical simulations to test the reliability of our results. The initial condition is displayed in the small inset of Figure 1: here $`c_2=0`$, thus before the time of first crossing, only the contribution of the growing mode has to be considered. In the main plot of Figure 1 the absolute value of the particle velocity, $`\beta _{cross}`$, is represented as function of rescaled time $`\theta `$, in scale lin-log. The solid line refers to the analytical prediction (35). In Figure 2 the particle maximum separation $`\mathrm{\Delta }q_{max}(\theta )`$ is plotted versus $`\theta `$, in linear-logarithmic scale (36). The solid line shows the theoretical result (36). In both case good agreement is displayed, even in the initial stage of the dynamics, thus confirming a posteriori the validity of our assumptions.

## 5 Density profile: Q model vs. adhesion model

The objective of this section is to investigate the density profile of a collapsing cluster in the Q model, and compare that to the adhesion model. We show that collapse in the Q model, which we recall is nothing but one-dimensional self-gravitating dynamics on the background of an Einstein-de Sitter Universe, is less pronounced. We derive a general expression for the density profile. To this end, we begin by analysing the motion a particle in an external halo, moving in the potential from a much more massive and very localized agglomeration. Then, we reconstruct the whole density distribution by combining together the contribution from each particle. These theoretical results are then compared first with numerical simulations of the Q model, and then with predictions based on the Zeldovich approximation and the adhesion model.
Consider a system of $`N`$ particles of equal mass, confined in box of finite size $`L=1`$, symmetric with respect to a reference point located at the origin. First, assume $`N2l`$ particles to form an inner bulk of width $`\mathrm{\Delta }`$, localized in the center of the box, with $`l`$ much smaller than $`N`$, and $`\mathrm{\Delta }`$ much smaller than $`L`$. The remaining $`2l`$ particles are assumed uniformly distributed in the outer regions at the borders of the box.
To first approximation we can neglect the fine structure of the cluster, and replace it with an heavy, structure-less, macro-particle, carrying the same amount of mass. In addition, the gravitational interaction between the particles in the halo is also neglected. As a consequence, these particles are then not affected by the discreteness of the inner distribution: a sudden change in the acceleration is experienced when the origin is crossed.
As in the previous section, we set initial conditions such that only the growing mode is active, before the first particles crossing occurs, and use the same dimension-less variable. Consider the motion of particle $`i`$, situated to the left of the agglomeration and define:

$$\gamma _i=\frac{1}{2}\frac{1}{N}(i\frac{1}{2}).$$
(37)
The positive parameter $`\gamma _i`$ represents the absolute value of the initial position of the $`i`$’th particle ($`w_i^0=0`$ in (23)). In analogy with (28)-(31), one finds:

$$\theta _{i,cross}=\frac{3}{2}\mathrm{ln}\left(\frac{\beta _i}{4\gamma _i}\right),$$
(38)
$$\beta _{i,cross}=\beta _i(\theta _{cross})=4\gamma _i,$$
(39)
$$\theta _{i,turn}=\theta _{i,cross}+\stackrel{~}{\theta }=\theta _{i,cross}+\frac{3}{5}\mathrm{ln}6,$$
(40)
$$q_{i,max,0}=4\alpha \gamma _i.$$
(41)
In this approximation, the $`i`$’th particle crosses for the first time the origin with rescaled velocity $`\beta _{i,cross,0}`$. Note that the latter is solely determined by the particle initial position. The corresponding maximum elongation $`q_{i,max,0}`$ is $`4\alpha 0.32`$ times smaller then the unperturbed distance $`\gamma _i`$.
By repeating the same procedure as in Section 4, one ends up with the following relations for, respectively, $`|\beta _{i,cross}|`$ and $`|\mathrm{\Delta }q_{i,max}|`$:

$$|\beta _{i,cross}(\theta )|=4\gamma _i\mathrm{exp}\left(\frac{1}{9}(\theta \theta _{i,cross})\right),$$
(42)
$$|\mathrm{\Delta }q_{i,max}(\theta )|=4\alpha \gamma _i\mathrm{exp}\left(\frac{2}{9}(\theta \theta _{i,turn})\right).$$
(43)
To validate our analysis we performed numerical simulations and studied the single particle behavior. In our experiments we considered the following initial condition: a finite, but large, fraction of all the particles, say $`N2\mathrm{}`$, is uniformly distributed in a narrow region, say $`\mathrm{\Delta }`$, centered around the reference origin. The remaining particles, $`\mathrm{}`$ of them, are symmetrically placed in the two lateral regions, with spatially uniform distribution. As a consequence, the inter-particle distance in the central bulk is smaller than in the external halos. In the limit when $`\mathrm{\Delta }0`$ and $`\mathrm{}=1`$, we reach the conditions assumed in the preceding derivation. Velocities are given as a smooth function of positions.
In the upper panel of Figure 3 a phase space portrait of a single particle is represented: the damping, both in $`q`$ and $`\beta `$, is clearly displayed. In the lower panel, the rescaled particle position $`q`$ is plotted vs. time $`\theta `$: the thin solid line refers to the simulation, while the thick curve is plotted from relation (43). Numerical simulations show good agreement with analytical predictions (42) and (43).
The previous results can also be extended to the case when particles are initially uniformly distributed in the finite size box. Given particle j, we neglect the interaction with the external masses (i.e. $`|q_j|>|\gamma _i|`$) and mimic the attraction of the inner particles by considering the forces exerted by a single massive object, localized in the origin. In this approximation, according to the lines of the preceding discussion, the behaviors of the $`|\beta _{i,cross}(\theta )|`$ and $`|\mathrm{\Delta }q_{i,max}(\theta )|`$ are also described by equations (42) and (43). Despite the drastic assumptions involved, numerical simulations show relatively good agreement with theoretical predictions, provided the velocity distribution are smooth enough.

The goal of the remaining discussion is to derive an analytical estimate of the density profile by making use of these results. In the continuum limit ($`N\mathrm{}`$), the density $`\rho (q,\theta )`$ reads:

$$\rho (q,\theta )=\rho _0/|\mathrm{\Delta }q/\mathrm{\Delta }\gamma |,$$
(44)
where we introduced the Jacobian of the transformation from Lagrangian to Eulerian coordinates. The density field has a rather complex structure including singularities $`\rho (q,\theta )(qq_m)^{1/2}`$, where $`q_m`$ stands for the particle turning point. A rough characteristics of the cluster size, $`q_s`$, can be deduced from equation (41). In fact, before the time the outmost particle reaches, as the last one, the massive bulk, the size of the agglomeration, $`q_s`$, is given by:

$$q_s(\theta )=4\alpha \gamma _i(\theta ),$$
(45)
where $`\gamma _i(\theta )`$ is the initial coordinate of the particle which falls into the origin at the time $`\theta \stackrel{~}{\theta }`$. From (38), (40) it follows:

$$\gamma _i(\theta )=\frac{\beta _i}{4}\mathrm{exp}\left(\frac{2}{3}(\theta \stackrel{~}{\theta })\right)=\frac{\beta _i}{6^{2/5}4}\mathrm{exp}\left(\frac{2}{3}\theta \right)$$
(46)
When the continuous limit is recovered, $`\beta _i=\beta (\gamma _{})`$, $`\gamma _{}`$ being the Lagrangian coordinate of the most periferic particle of the cluster, at time $`\theta `$. Then assume that $`\beta (\gamma )0`$, as the edge of the box is approached ($`\gamma 1/2`$). Then from (46), $`\gamma _{}`$ tends to $`1/2`$, at large times, and consequently, from (45), $`q_s(\theta )2\alpha 0.16`$. This observation on is in good agreement with the results of numerical simulations presented in a in previous work .
Inside the pancake a multistream flow has developed: the density at coordinate $`q`$ follows by summing the contribution of each stream. Consider the mean density distribution assuming the length scale $`\mathrm{\Delta }q`$ much larger than the distance between two successive peaks i.e. larger than $`|q_{m+1}q_m|`$. The sum can be then approximated by an integral over the Lagrangian coordinate $`\gamma `$. By introducing the density distribution function $`f_\rho (\gamma ,\theta )`$ associated to the Lagrangian coordinate $`\gamma `$, we can write:

$$\rho (q,\theta )=_{\gamma _{min}}^{\gamma _{max}}f_\rho (\gamma ,\theta )𝑑\gamma .$$
(47)
Assume that all the particles belonging to a small Lagrangian interval, $`\mathrm{\Delta }\gamma `$, are uniformly distributed inside the corresponding Eulerian segment $`\mathrm{\Delta }q_{max}(\theta ,\gamma )`$ (43):

$$\mathrm{\Delta }q_{max}(\theta ,\gamma )=|\mathrm{\Delta }q_{i,max}(\theta )|=\frac{\chi \gamma ^{\frac{4}{3}}}{\beta ^{\frac{1}{3}}(\gamma )}\mathrm{exp}\left(\frac{2}{9}\theta \right),\chi =4^{4/3}6^{2/15}\alpha $$
(48)
and consider the mass conservation:

$$\rho _0\mathrm{\Delta }\gamma =\left(f(\gamma ,\theta )\mathrm{\Delta }\gamma \right)\mathrm{\Delta }q_{max}(\theta ,\gamma ),$$
(49)
The interval $`[\gamma _{min}(q,\theta ),\gamma _{max}(\theta )]`$ selects the particles that contributes to the density in $`q`$, at time $`\theta `$. In particular they should have reached the origin, thus:

$$\frac{\beta (\gamma _{max})}{4}\mathrm{exp}\left(\frac{2}{3}\theta \right)=\gamma _{max}.$$
(50)
On the other hand the amplitudes of the particles’ oscillations, $`\mathrm{\Delta }q_{max}(\theta ,\gamma )`$ have to be larger than $`q`$, thus $`\gamma _{min}`$ is found by solving the following implicit equation.

$$q=\chi \mathrm{exp}\left(\frac{2\theta }{9}\right)\frac{\gamma _{min}^{\frac{4}{3}}}{\beta ^{\frac{1}{3}}(\gamma _{min})}$$
(51)
Finally, the density reads:

$$\rho (q,\theta )=\rho _0_{\gamma _{min}}^{\gamma _{max}}\frac{d\gamma }{\mathrm{\Delta }q_{max}(\theta ,\gamma )}=\frac{\rho _0}{\chi }\mathrm{exp}\left(\frac{2}{9}\theta \right)_{\gamma _{min}}^{\gamma _{max}}\frac{\beta ^{\frac{1}{3}}(\gamma )}{\gamma ^{\frac{4}{3}}}𝑑\gamma ,$$
(52)
and $`\rho (q,\theta )=0`$ when:

$$|q|>q_w(\theta )=\mathrm{\Delta }q_{max}(\theta ,\gamma _{max}).$$
(53)
Here $`q_w(\theta )`$ is the width of the mean density distribution. The density profile (52) depends on the initial particle velocity, namely $`\beta (\gamma )`$. In the following, we will solve analytically the integral in (52) for special class of initial conditions, and compare with numerics. In addition, we will discuss the analogous case in the adhesion model.

### 5.1 Step Profile of Initial Velocities

Assume the initial velocities profile to be given by:

$$\beta (q)=\beta _0\text{sign}(q),$$
(54)
where $`\beta _0`$ is a constant coefficient and, without loss of generality, we focus on the region $`q>0`$.
First consider the evolution of the density in the Zeldovich approximation, with initial positions as in equation (23). Since the velocity profile at $`\theta =0`$ is a step function, the particles with initial negative (resp. positive) spatial coordinate will move as a whole towards the origin, keeping their mutual distances constant. Thus before the last left particle reaches the origin, a two-stream flow develops over a finite region of width $`2|q_{w,zeld}|`$, where $`|q_{w,zeld}|=\beta _0\mathrm{exp}(2\theta /3)/4`$. Note that the the density in each flow is equal to the initial density $`\rho _0`$.
Then follow the development under the Q dynamics later in time. By inserting (54) in (52) and performing the integral, one obtains:

$$\rho (q,\theta )|_{|q|<q_w}=\rho _0\frac{3}{2q_p}\left[\left(\frac{q_p}{q}\right)^{1/4}\left(\frac{1}{2\gamma _{max}}\right)^{1/3}\right],$$
(55)
where $`\gamma _{max}`$ is solution of equation (50) and $`q_p(\theta )`$ is the spatial scale of inner structure

$$q_p=\frac{\chi }{2^{\frac{4}{3}}\beta _0^{\frac{1}{3}}}\mathrm{exp}\left(\frac{2}{9}\theta \right).$$
(56)
We recall that $`q_w(\theta )`$ stands for the width of the pancake, see eq. (53). Consider also the mass function $`M(q,\theta )`$:

$$M(q,\theta )=_0^q\rho (\xi ,\theta )𝑑\xi ,$$
(57)
Inserting eq. (55) and performing the integral one ends up with:

$$M(q,\theta )=\rho _0\left[2\left(\frac{q}{q_p}\right)^{3/4}\frac{3}{2}\left(\frac{1}{2\gamma _{max}}\right)^{1/3}\left(\frac{q}{q_p}\right)\right].$$
(58)
In the initial stage of the evolution $`\beta (\gamma _{max})=\beta _0`$, thus, from equation (50):

$$\gamma _{max}=\frac{\beta _0}{4}\mathrm{exp}\left(\frac{2}{3}\theta \right).$$
(59)
The pancake’s density distribution is therefore characterized by two length scales: the spatial scale of inner structure, $`q_p(\theta )`$, which decays exponentially, see (56), and the width of the agglomeration which grows exponentially according to:

$$q_w(\theta )=\frac{\chi \beta _0}{4^{\frac{1}{3}}}\mathrm{exp}\left(\frac{2}{3}\theta \right)0.1\beta _0\mathrm{exp}\left(\frac{2}{3}\theta \right).$$
(60)
It is worth stressing that during the initial stage of the pancake’s formation, the width $`q_w`$ increases with time slower than as predicted in the Zeldovich picture. From equation (55) it follows that the density in the bulk increases exponentially in time-like coordinate $`\theta `$, according to $`\rho (q,\theta )q^{(1/4)}\mathrm{exp}(\theta /6)`$. In other words, the system tends to display a progressively denser core, localized near the origin. The mass contained in the cluster is also increasing exponentially, $`m_p(\theta )=2M(q_w,\theta )=\rho _0\beta \mathrm{exp}(2\theta /3)/2`$, in agreement with the predictions of the Zeldovich model. The discussion above applies to the co-moving frame, but physical density is measured in the inertial frame; the two are related by (1). We will return to this point in Sect. 7. Let us however in brief state that the collapse in co-moving coordinates is slower than the expansion of the Universe. The physical density therefore on the contrary decreases at around almost all mass points, and in an inertial coordinate system the agglomerations spread out over time, albeit more slowly than the expansion of the Universe as a whole.
Equations (58)-(60) suggest that, initially, both the mean density distribution and mean mass function are self-similar. In particular:

$$M(q,\theta )=m_p\left[2\left(\frac{q}{q_w}\right)^{3/4}\frac{3}{2}\left(\frac{q}{q_w}\right)\right]=m_p(\theta )\overline{M}(q/q_w(\theta )).$$
(61)
In the very late time evolution of an isolated cluster all particles have reached the inner bulk (i.e. $`\theta >3/2\mathrm{ln}(\beta _0/2)`$ ). The box length, $`L`$, is normalized to one so that, $`m_p(\theta )=\rho _0`$. Then $`\gamma _{max}=1/2`$ and equations (55),(57) indicate the a self-similar collapse:

$$M(q,\theta )=\overline{M}(q/q_p(\theta )).$$
(62)
Now both the density and the mass function are solely characterized by a unique length scale, which is both the inner length scale, $`q_p(\theta )`$, and the outer length scale, $`q_w`$. In this case $`q_p`$ is the coordinate of the out-most particle, and the density for $`q>q_p`$ is identically equal to zero.
Summing up, we have shown that for an initial condition in the form of equation (54), the mean mass function $`M(q,\theta )`$, displays a self-similar collapse, in different stages of the evolution. Introducing the rescaled variable $`x=q/q_w(\theta )`$, the function $`\overline{M}(q/q_w(\theta ))`$ takes the universal form:

$$\overline{M}(x)=\left[2x^{3/4}\frac{3}{2}x\right].$$
(63)
for $`x=[0,1]`$, where $`\overline{M}(1)`$ is equal to one half. From (63) it follows that the matter is mainly concentrated around the pancake center. In particular, one-half of the whole mass (i.e. $`2\overline{M}(x_{0.5})=1/2`$) is found to lay in a symmetric interval defined by $`x_{0.5}=0.14`$, while $`90\%`$ of it, is distributed in a segment equal to half the size of entire cluster ($`x_{0.9}=0.54`$).

To validate our theoretical analysis we now turn to numerical simulations and consider the late stage of evolution of an initial perturbation in the form (54). The particles are initially uniformly distributed in space. In Figure 4 the normalized density profile is plotted at late stage of evolution (thin solid line) and superposed to the theoretical prediction (55) (thick solid line). The same curves are represented in scale log-log in the left inset: here the circles refer to the simulation. No parameter need to be adjusted by numerical fitting and the agreement has to be considered satisfying. The corresponding phase space portrait is reported in the right inset of Figure 4, where the characteristic spiral behavior is clearly displayed. In figure 5 the cumulative mass function is plotted as a function of the variable $`\xi =q/q_p`$.

In a previous work of two of the authors, , we studied the late time evolution of an isolated perturbation within the framework of the Q dynamics. In particular we measured the progressive contraction of the inner region of the agglomeration, when compared to the overall Universe expansion. This effect was computed by the width of region, $`\mathrm{\Delta }q_\mu `$, that contains a fraction $`\mu `$ of the whole mass of the system, centered around the position of maximum density. It turned out that the interval $`\mathrm{\Delta }q_\mu `$ shrinks in time according to a power-law of cosmological time $`t`$, with exponent $`2/9`$. An heuristic interpretation was also provided in . To make a bridge between our present investigations and these earlier findings we define $`q_\mu `$ (63) such that $`2M(q_\mu ,\theta )=\rho _0\mu =const`$ . From equation (62) it follows that $`\mathrm{\Delta }q_\mu `$ scales proportionally to $`q_p`$; thus we are led to assume $`q_\mu =x_\mu q_p`$, where $`x_\mu `$ is a positive coefficient. Equation(56) implies:

$$q_\mu =\frac{\chi x_\mu }{2^{\frac{4}{3}}\beta _0^{\frac{1}{3}}}\mathrm{exp}\left(\frac{2}{9}\theta \right)$$
(64)
in qualitative agreement with the results in . As an example, let us focus on the case $`\mu =1/2`$. In this case, $`x_\mu =0.14`$. We performed numerical simulations starting with the same class of initial condition as in Figure 4. In Figure 6, $`q_{1/2}`$ is plotted vs. $`\theta `$: circles refer to the numerical simulation and the solid line represents equation (64). The slopes of the two curves in figure 6 agree, but the amplitudes do not. We can make the difference smaller by adjusting some of the hypotheses made in the derivation above. In particular, we considered the particles belonging to the cluster to be uniformly distributed in the finite interval $`\mathrm{\Delta }q_{max}`$ given by (48). In reality, the particles spend more time on the border. This effect is explicitly modeled by assuming:

$$q_{p,eff}=cq_p$$
(65)
where $`c`$ is an “ad hoc” numerical factor, larger than one.

To complete this section, let us compare our theoretical expression (55) with the density profile derived in the framework of the adhesion model. These results, which we recall in the following, are detailed in . Let as consider the stationary solution of Burgers equation :

$$v(x,b)=v_{st}(x)=U\mathrm{tanh}\left(\frac{x}{\delta }\right),$$
(66)
where $`\delta =U/2\nu `$ is the width of the shock. This solution is also the asymptotic solution of Burgers equation for the initial the step-function profile previously considered. The trajectories of individual particles $`x(b)`$ satisfy the following equation

$$\frac{dx(b)}{db}=v(x,b),$$
(67)
where the field $`v(x,b)`$ is determined by the solution of the Burgers equation. Assuming that $`v(x,b)`$ is a stationary solution of Burgers equation from (66),(67) we have the following expression for the coordinates of the particles

$$x(b,y)=\delta \text{Arc}\mathrm{sinh}\left[\mathrm{sinh}(y/\delta )\mathrm{exp}(Ub/\delta )\right]$$
(68)
where $`y`$ stands for the Lagrangian coordinate. From the conservation of mass, the following expression for the density is derived:

$$\rho (x,b)=\frac{\rho _0\mathrm{cosh}\left(\frac{x}{\delta }\right)}{\sqrt{\mathrm{sinh}^2\left(\frac{x}{\delta }\right)+\mathrm{exp}\left(\frac{2Ub}{\delta }\right)}},$$
(69)
Solutions (68) do not apply to the intermediate stage of evolution. Thus we will focus on the asymptotic behavior, for time $`b\delta /U`$ and consider the regions $`|x|>\delta _p`$, where $`\delta _p=\delta \mathrm{exp}(Ub/\delta )`$. The particles in the agglomeration behave according to

$$x(b,y)=y\mathrm{exp}(Ub/\delta ).$$
(70)
Thus in the adhesion model a faster collapse is displayed, when compared to the exact Q dynamics: in the latter a power-law decay in “Burgers time” $`b`$, is in fact produced, namely $`q_p(b,y)=y(b/b_0)^{1/3}`$. As concern the density profile, for $`x<\delta _p`$, one gets:

$$\rho (x,b)=\frac{\rho _0\mathrm{exp}\left(\frac{Ub}{\delta }\right)}{\sqrt{1+\left(\frac{x^2}{\delta ^2}\right)\mathrm{exp}\left(\frac{2Ub}{\delta }\right)}},$$
(71)
hence maximum value of the density increases exponentially in “Burgers time”, $`\rho _{max}(0,b)=\rho _0\mathrm{exp}(Ub/\delta )`$. In the interval $`\delta _px\delta `$ the density transforms into a time-independent power-law distribution

$$\rho (x,b)=\rho _0\frac{\delta }{x}=\rho _0\frac{2\nu }{Ux}.$$
(72)
Thus in the adhesion model we have that the asymptotic density distribution for the initial step profile is localized in the regions $`|x|<\delta =const`$, has non-integrable time independent power law tails (72) with cutoff scale $`\delta _p=\delta \mathrm{exp}(Ub/\delta )`$. In the Q model we have instead an integrable power-law distribution, see (55).

When considering a perturbation in a finite box, in co-moving coordinates, both the Quintic and the adhesion models show the formation of a dense structure. However, in the adhesion model, the width $`q_p`$ of the pancake decreases faster ( $`g_p\mathrm{exp}b`$) then in Q-model ( $`g_pb^{1/3}`$ ).
These discrepancies suggest that, at least in one dimension, the adhesion approach is valid only as an approximate model of structure formation. This observation agrees with the conclusion in  where a new transport equation was proposed.

## 6 Interaction of two massive clusters

In this section we consider the dynamics of two interacting clusters. As a first approximation, one can neglect the inner discrete structure of a cluster, by introducing a single macro-particle to mimic its dynamics. Thus, the problem is first reduced to the study of the evolution of a two particles system, investigated in Section 4. According to this picture, the agglomerations move back and forth displaying damped oscillation. This observation agrees with previous numerical investigations reported in . Following the discussion of the previous section, each massive agglomeration experiences a contraction of the inner core: the competition of these two effects determines the late-time asymptotics of the system.

Numerical simulation are performed starting with the initial condition displayed in Figure 7 a). The distance between the two centers of mass, $`\delta `$, is plotted vs. time $`\theta `$, in Figure 7 c). The damping tendency is evident (thin solid line) and the theoretical prediction (36) (thick solid lines) agrees with the numerical findings. At late time the two clusters still present a finite separation, see 7 b). This is in a good agreement with the theoretical predictions of Sections 4 and 5.1, both the widths of the agglomerations and amplitude of their mutual oscillations decay as the same exponential function of time, hence their ratio is a constant number. Therefore, the system composed by two interacting clusters, with finite inner structure, evolves approximately in a self-similar way.

## 7 Concluding remarks

In this paper we considered the problem of structure formation in a Einstein-de Sitter Universe, focusing on planar perturbations. This was done by using the Q model, or Quintic model, a Lagrangian representation previously derived in , in which one of the steps is to solve a large number of quintic equations. We underline however that the Q model is nothing but the motion of a system of gravitationally interacting particles, on the background of an expanding Universe, with convenient choices of time and space coordinates. We showed that the Q model can be approximated by an even simpler iterative map. This observation permitted to establish a connection with  where the motion of turbulent flows with strong friction was analyzed, and to derive theoretical predictions for both the density profile and the mass function in a collapsing cluster. These predictions were cross-validated with direct numerical simulation.
The main conclusion of this work is that the inner structure of a Zeldovich pancake is different from the adhesion model, with finite viscosity in the Burgers equation. That there is a discrepancy is what one would expect, but this had nevertheless to our knowledge never been shown before. As a mathematical problem, a qualitative difference is that in co-moving coordinates, the size of collapsing cluster goes down as power-law with cosmological time in the self-gravitating system, but as a stretched exponential in the adhesion model <sup>1</sup><sup>1</sup>1In “Burgers time” the collapse is simply exponential, see subsection 5.1.. As a physical problem, the difference is even more dramatic: in the self-gravitating system the characteristic scale of the density distribution is given in co-moving coordinates by equation (56), which decreases, with cosmological time $`t`$, as $`\left(t/t_0\right)^{\frac{2}{9}}`$. The co-moving system of coordinates is related to an inertial system by equation (1), where the proportionality is $`\left(t/t_0\right)^{\frac{2}{3}}`$. Combining the two, the characteristic scale of a collapsing cluster hence increases in time, as $`d(t)\left(t/t_0\right)^{\frac{4}{9}}`$. Since the density distribution of equation (55) is self-similar, the maximum density attained depends on the graininess of the initial conditions. We have not investigated the law of the increase of the maximum density, but it is at least possible that with continuous initial density distributions arbitrarily high densities can be produced, also if measured in an inertial system of reference. A more interesting characteristics is however the fraction of total mass $`M`$ where (in an inertial system) the density is higher than some cut-off $`\rho _{cr}`$. That quantity scales as $`\left(\frac{M}{\rho _{cr}d(t)}\right)^3`$. Since $`d(t)`$ is an increasing function of time in the inertial system, that fraction decreases as $`\left(t/t_0\right)^{\frac{4}{3}}`$. In this sense one-dimensional self-gravitating motion, on the background of an expanding Universe, is not effective in bringing about mass collapse.
In the adhesion model the characteristic scale decreases, both in a co-moving and an inertial system. Here both the maximum real-space density and the fraction of mass at a density larger than some cut-off increases. The adhesion model is therefore, in contrast, effective in bringing about mass collapse in one dimension. In other words, the adhesion model over-predicts the clustering tendency, and overestimates the peak density. A Universe in which most mass lies in regions where density decreases can never give rise to galaxies, galaxy clusters and other significant objects. The analysis given here therefore suggests, quantitatively, why a one-dimensional description of mass agglomeration is insufficient.
We have also investigated the mutual motion of two pancakes, and shown them to be closely similar to the decaying oscillations of two mutually attracting “heavy particles”.

## 8 Acknowledgements

This work was mainly supported by a grant from the Swedish Royal Academy of Sciences (KVA), further support from the Swedish Research Council through grants NFR F 650-19981250 (D.F) and NFR I 510-930 (E.A.), and RFBR (02-02-17374, 00-15-96619) and ”Universities of Russia” (S.N.G. and A.Yu.M) is also gratefully acknowledged. E.A. thanks the Edwin Schrödinger Institute (Vienna) for hospitality, and for opportunities to discuss and write up this work. D.F. thanks Ph. Choquard for enlighting discussions. We thank A. Noullez for fruitful discussions.