Nonlinear dynamics of electron-positron clusters

Electron-positron clusters are studied using a quantum hydrodynamic model that includes Coulomb and exchange interactions. A variational Lagrangian method is used to determine their stationary and dynamical properties. The cluster static features are validated against existing Hartree-Fock calculations. In the linear response regime, we investigate both dipole and monopole (breathing) modes. The dipole mode is reminiscent of the surface plasmon mode usually observed in metal clusters. The nonlinear regime is explored by means of numerical simulations. We show that, by exciting the cluster with a chirped laser pulse with slowly varying frequency (autoresonance), it is possible to efficiently separate the electron and positron populations on a timescale of a few tens of femtoseconds.


Introduction
The positron was the first antiparticle to be discovered experimentally (in 1932) and it appears naturally as a negative energy solution of the Dirac equation [1]. Positron physics is of great fundamental and practical interest, ranging from condensed matter physics to astrophysics and biological physics. For instance, positron techniques are useful for the investigation of defects in solids and solid surfaces [2]. In medicine, the widespread use of positron-emission tomography (PET) for diagnostics and treatment monitoring requires a sound understanding of the physical and biological effects of positrons on living organisms. In astrophysics and cosmology, understanding the imbalance of matter versus antimatter is one of the major challenges of today's theoretical physics. Finally, recent projects aiming to elucidate the gravitational behavior of antimatter require the careful manipulation of positrons in order to produce anti-hydrogen atoms [3].
Positrons can be easily obtained from the β + decay of radioactive isotopes, e.g. from 22 Na. The positrons generated in this reaction exhibit a broad energy spectrum that extends up to 540 keV. For practical use in antimatter studies, positrons need to be cooled down to a few electron-volts by means of a moderator and are subsequently stored in a trap [3]. Slow positrons implanted into a porous silica film may efficiently form positronium atoms (Ps) [4,5]. Positronium is a bound state constituted of an electron and a positron. It is the lightest particle-antiparticle "atom", with a relatively long lifetime of 125 ps for the singlet state (para-positronium) and 142 ns for the triplet state (ortho-positronium).
Positronium may also be viewed as the simplest "many-body" electron-positron (ep) system. For larger numbers of particles, various other states are possible, depending on the density and temperature of the system [6]. At high temperatures and low densities, classical e-p plasmas can be formed, both relativistic and nonrelativistic [7,8,9,10]. Lowering the temperature below the Ps ground-state binding energy (E Ps = 6.8 eV), the electrons and the positrons recombine to form a positronium gas. At still lower temperatures, Ps atoms may even form a Bose-Einstein condensate (BEC) [11]. Because of its light mass, the critical temperature of a Ps gas is much higher than, for instance, that of an alkali atom gas with the same number density, which is obviously an interesting feature from an experimental point of view. For instance, for a Ps density n = 10 24 m −3 , the critical temperature would be as high as T c ≈ 15 K [12]. The realization of a BEC of Ps atoms is a very exciting and challenging project, as such a system could lead to the simultaneous coherent decay of all Ps atoms, thus acting as a powerful gamma ray source.
The critical temperature T C ∼ 2 n 2/3 /(mk B ) (where m is the electron mass) is attained when the de Broglie thermal wavelength of the Ps atoms, λ B = / √ mk B T , becomes comparable to the mean interparticle distance, measured for instance by the Wigner-Seitz radius r s = (4πn/3) −1/3 . At even higher densities, when the interparticle distance is of the order of the Ps ground-state radius 2a 0 (where a 0 is the Bohr radius), the e-p system effectively behaves as a two-component degenerate Fermi liquid [6]. This is the "metallic" phase of electron-positron matter. For a wide range of values of r s , the ground-state properties of e-p infinite matter were studied by Boronski and Nieminem [13] using two-component density-functional theory (DFT) within the local density approximation (LDA).
However, it is well known that finite-size metallic systems can also exist. These systems -known as metal clusters [14,15,16] -are usually composed of 10 ≤ N ≤ 10 4 ions and an equal number of electrons (although charged clusters have also been studied). In many respects, they behave as giant "atoms" where the positive charge is not localized in the nucleus but is distributed more or less uniformly within the cluster. For metal clusters, the positive charges are ions, whose mass is thousands of times that of the electrons. Thus, it is appropriate, when studying effects occurring on a short timescale (< 100 fs), to assume that the ions are effectively immobile. Further, for large clusters the ions can be modeled by a uniform positive charge density (jellium approximation). The electronic ground state, obtained by DFT methods, reveals a shell structure with discrete energy levels, akin to those of ordinary atoms.
Recently, it has been suggested [17,18] that clusters made of electrons and positrons could also exist. Using either two-component DFT or Hartree-Fock methods, it was shown that such e-p clusters have ground-state densities similar to that of metals (r s /a 0 ≈ 3.5) and also display an electronic shell structure. Of course, for e-p clusters it is not possible to use the jellium approximation, as both species have the same mass and are thus equally mobile. For the same reason, e-p clusters are locally neutral in their ground state, whereas metal clusters display a "spill-out" effect [15], whereby the electron density extends slightly further than the ion density.
The experimental realization of stable e-p clusters is still ahead of us, mainly because extremely high densities are required, n ≈ 10 28 m −3 (at such densities, the Fermi temperature is very high, T F ≈ 10 4 K, so that the e-p gas is degenerate even at room temperature). However, remarkable advances have been made in the confinement and cooling of positron plasmas. Temperatures smaller than 5 K and densities larger than 10 16 m −3 have been achieved in recent years [19,20]. Further, recent technological progress is opening up the possibility of employing intense laser radiation to trigger physical processes beyond atomic-physics energy scales, such as electron-positron pair production at high densities [21,22]. Indeed, very recently, in realistic simulations of a 10 PW laser striking a solid target, it was demonstrated that a maximum positron density of 10 26 m −3 can be obtained, seven orders of magnitude greater than achieved in previous experiments [23].
On the other hand, dense gases of interacting Ps atoms have been created by irradiating a thin film of nanoporous silica with intense positron bursts [12], reaching densities of the order of n ≈ 10 21 Ps/m 3 , just three orders of magnitude lower than the density needed to form a BEC of Ps atoms with T C ≈ 15 K. All in all, both Ps BECs and e-p clusters represent admissible states of electron-positron matter under extreme conditions of temperature and density and in that respect they deserve to be properly investigated theoretically.
As of today, no results exist on the dynamics of e-p clusters, either in the linear or nonlinear regimes. In contrast, the linear response of metal clusters has been the subject of intense investigations in the last few decades [24,25]. A strong dipole resonance is observed near the Mie frequency, which for spherical clusters in the jellium approximation is given by the bulk plasma frequency divided by √ 3. Using more sophisticated approaches, it can be shown that the resonant frequency actually depends on the cluster size. The nonlinear electronic response was investigated more recently by means of phase-space or hydrodynamic methods [26].
In the present work, we aim at characterizing the linear and nonlinear response of e-p clusters. We use a variational approach based on a two-component quantum hydrodynamic method that incorporates the kinetic energy, the Coulomb interaction, and the exchange energy, but neglects higher-order correlations (Sec. 2). Thus, our method can be viewed as an approximation of the two-component Hartree-Fock equations. Only one adjustable parameter (related to the gradient correction of the exchange energy) appears in our model and is determined by matching our solution for the ground state (Sec. 3) with that obtained through Hartree-Fock calculations. Subsequently, we study the linear response of the e-p clusters, which reveals both dipole and monopole resonances (Sec. 4). Finally, we investigate numerically the nonlinear dynamics and show that the electron and positron populations can be effectively separated on a timescale much shorter than that of mutual annihilation (Sec. 5). Conclusions are presented in Sec. 6.

Quantum hydrodynamic model and Lagrangian method
In our approach, the electron-positron system is governed by a set of quantum hydrodynamic equations for the densities n i and the mean velocities u i , where the subscript i = e, p denotes each species [27]. Hydrodynamic methods have been used successfully in the past to model the electron dynamics in molecular systems [28], metal clusters and nanoparticles [29,30], thin metal films [31], and quantum plasmas [32,33]. In the following, we will always use atomic units (a.u.) such that space is normalized to the Bohr radius a 0 = 4πε 0 2 /(me 2 ), energy to the Hartree energy E H = me 4 /(4πε 0 ) 2 , and time to τ H = /E H . In a.u., the hydrodynamic equations read as follows: where q i = ±1 for positrons and electrons. The four terms on the right-hand side of Eq. (2) represent respectively the pressure, the Hartree (Coulomb) potential, the exchange potential, and the von Weizsäcker correction (sometimes referred to as the Bohm potential).
We neglect correlations altogether in our model. This is done for two main reasons: first, the correlation functional for finite-size systems of same-mass particles is extremely complicated, particularly the cross terms (electron correlations due to the positron cloud, and vice-versa) [13]; second, results that only include the exchange term can be directly compared to exact calculations performed with the Hartree-Fock equations [17].
The Hartree potential V H satisfies Poisson's equation The exchange potential is derived from the exchange energy functional where the first term is the usual LDA expression and the second term is a gradient correction [34]. The latter is on the same level of approximation as the von Weizsäcker correction to the kinetic energy. The parameter β will be determined numerically by comparing our results to exact Hartree-Fock calculations (see Sec. 3). The obtained value β = 0.0135 is larger than the best-fit value usually employed in atomic-structure calculations, which is β ≈ 0.005 [34]. The exchange potential is then obtained as the functional derivative of the exchange energy with respect to the density: Finally, for the pressure we use the expression of the Fermi pressure for a zerotemperature electron (or positron) gas: where E F [n i ] is the Fermi energy. Since the ground state density of e-p clusters is similar to that of metals (see next section), their Fermi temperature is of the order T F ∼ 10 4 K, which justifies the zero-temperature assumption. It can be shown that the hydrodynamic equations (1)-(2) can be derived from the following Lagrangian density L [35]: where the independent fields are n i , S i , and V H . The velocity fields u i follow from the auxiliary functions Our purpose now is to derive -using the variational approach detailed in Ref. [35] -a set of evolution equations for a small number of macroscopic quantities that characterize the electron and positron density profiles. With this aim in mind, we assume that the density profiles are Gaussian functions where N is the total number of electrons and positrons (assuming overall charge neutrality), whereas σ i (t) and d i (t) are time-dependent functions defining respectively the size of the electron and positron clouds and the displacement in the z direction. We allow for a displacement along the z axis because we will later suggest using a laser pulse to excite the dipole mode (see Secs. [4][5].
Of course, such a Gaussian Ansatz is not exact and may even differ significantly from the electron density obtained, for instance, from a HF calculation. Nevertheless, it is a useful and relatively safe procedure to obtain a mathematically treatable set of equations that can be solved either exactly or with minimal numerical effort. For instance, it was noticed in Ref. [35] that, even when the actual density is not well approximated by a Gaussian profile, the resonant frequencies computed with our technique are still very close to the exact ones.
For the above density profiles, the exact solution of Poisson's equation (2) is where Erf(s) is the error function and s 2 = x 2 +y 2 +(z −d i ) 2 . In addition, the continuity equation (1) is exactly solved by the following velocity field with where the dot denotes derivation with respect to time andx,ŷ,ẑ are unit vectors along each direction. An irrelevant additive function of time was ignored in Eq. (9). We can now compute the Lagrangian L = −N −1 L dr, where the multiplicative factor was introduced for convenience of notation. In atomic units, the Lagrangian reads as: where the coefficients correspond respectively to the kinetic energy (Fermi pressure), the Hartree energy, and the exchange energy. The corresponding equations of motion are obtained from the standard Euler-Lagrange equations d dt where ζ = {d e , d p , σ e , σ p }. The result is: where we have defined Σ 2 = σ 2 e + σ 2 p . It is preferable to use center-of-mass and relative coordinates, defined as D = (d e + d p )/2 and d = d e − d p . We obtain thatD = 0 and In the next sections, the above equations will be used to study the steady state properties, as well as the linear and nonlinear responses of e-p clusters.

Steady states
Yatsyshin et al. [17] and Solovyov et al. [18] have investigated the stationary properties of electron-positron clusters using respectively Hartree-Fock (HF) and DFT calculations. In all cases, they obtain neutral states where the electron and positron density profiles are locally identical. Within the framework of our model, we also look for neutral equilibria for which σ e = σ p = σ 0 and d = 0. Substituting into Eq. (15), the Hartree terms cancel out, because the equilibrium is locally neutral. We obtain: The steady state is then given by a Gaussian density with a width equal to The corresponding Wigner-Seitz radius r s can be defined using the peak value of the density: So far, our model still contains a free parameter, namely the coefficient β appearing in the gradient correction of the exchange energy functional, Eq. (3). For atomic systems, this parameter has usually been determined by comparison with exact HF calculations, yielding a value β ≈ 0.005 [34]. Here, we follow the same procedure and compare our analytical result, Eq. (17), with the HF calculations published in Ref. [17]. It turns out that the best fit on the Wigner-Seitz radius is obtained for β = 0.0135 and we will retain this value for all forthcoming calculations. Of course, oscillations of r s with the system size N -which are due to shell effects and are present in the HF calculations -cannot be recovered with our simple model.
The system size σ 0 and the Wigner-Seitz radius are shown in Fig. 1 as a function of the number of particles. As expected, σ 0 goes like N 1/3 for large N, whereas r s /a 0 varies between 3.45 in the bulk (N ≫ 1) and 3.9 for small values of N (although it should be pointed out that the validity of our approach becomes questionable for very small systems). In addition, the bulk value of r s does not depend on the chosen value of β.
These findings are in accordance with standard results obtained for metal clusters [26]. This is important, as one can envision to excite such e-p systems with the same optical means that are currently used to study the dynamical response of metal clusters.

Linear response
We now investigate the linear response of the e-p cluster under external excitations. Two types of mode can be studied in the framework of our model: (i) a dipole mode, where the electron and positron clouds oscillate with respect to each other and (ii) a breathing or monopole mode [35,36], where the radii σ i of both clouds oscillate, either in phase or in antiphase. It is important to stress that, at the level of the linear response, the dipole and breathing modes are completely decoupled, although of course some coupling will occur in the nonlinear regime.

Dipole mode
Let us rewrite the equation for the dipole d(t), assuming that σ e = σ p = σ 0 : The right-hand side of Eq. (19) can be written as −∂V /∂d, which implicitly defines the effective potential V (d). On Fig. 2 we plot the effective potential together with the standard Coulomb potential in vacuum: V Coul (d) = −2N/d. The factor 2N is the total number of charges (electrons and positrons). We can see that the effective potential is a "regularized" version of the Coulomb potential: the divergence at d = 0 has disappeared and we are left with a smooth potential well. We now linearize Eq. (19) around the equilibrium d = 0, using the expansion Erf(x) = (2/ √ π)(x − x 3 /3 + . . .). We obtaind + Ω 2 d d = 0, with the linear frequency: or in SI units: We now define the average density as [35] n ≡ n 2 (r)dr n(r)dr = N (2π) 3/2 σ 3 0 (22) and the plasma frequency computed with the reduced massm = m/2: The linear dipole frequency then becomes: which is the same as the Mie frequency for spherical metal clusters [15].

Breathing modes
We linearize Eqs. (14)-(15) around the equilibrium σ e = σ p = σ 0 and d = 0, by writing σ i = σ 0 +σ i and d =d, withσ i ,d ≪ σ 0 . We obtain that the equations ford and σ i completely decouple. The former yields the dipole mode described in the preceding subsection, while the equation forσ i reads as: By Fourier transforming in the time variable, i.e., replacingσ i with −Ω 2σ i , the linearized equations can be written as follows: The relevant dispersion relation can be written as A(Ω) = ±B, which yields the two resonant frequencies: In the large N limit, σ 0 → (2C K /C X ) N 1/3 , and we obtain We note that the solution Ω − corresponds to A = −B and thereforeσ e =σ p : this is a neutral mode where the electron and positron densities fluctuate in phase. In contrast, the solution Ω + corresponds toσ e = −σ p : it is a nonneutral mode, with the density fluctuations oscillating in antiphase. We further note that for the neutral mode Ω − only the Bohm, exchange and Fermi pressure terms play a role: these are all local terms, so that the mode disappears for very large clusters. In contrast, the nonneutral mode Ω + depends on the Hartree potential, which is nonlocal, hence this mode persists for all cluster sizes.

Nonlinear response and autoresonant excitation
We now turn our attention to the excitation of the electron and positron dynamics by means of electromagnetic waves (laser pulses). First, it should be noted that the relevant linear frequencies computed in the preceding section are of the order of a few electronvolts. For instance, for N = 100, one finds Ω d = 3.50 eV, Ω + = 4.31 eV, and Ω − = 0.45 eV. These frequencies fall within the visible or near ultraviolet (UV) spectrum, which is good news, as visible and near-UV lasers are commonly employed in ultrafast optics experiments. For such lasers, the wavelength is several hundred nanometers long, i.e., much larger than the size of a typical e-p cluster (see Fig. 1). This means that only the dipole mode can be excited directly, just like for ordinary metal clusters. The only hope to excite the breathing modes is via nonlinear coupling to the dipole mode.
In order to model the interaction of an electromagnetic wave with our e-p system, we introduce an external homogeneous electric field parallel to the z direction, E = E(t)ẑ. This amounts to adding a term −2E(t) on the right-hand side of Eq. (14) (the prefactor −2 appears because the force is −E for the electrons and +E the for positrons). Since the system is globally neutral, a homogeneous field does not affect the center-of-mass equation of motion, which remainsD = 0. Notice that this procedure is completely consistent with our Lagrangian approach and could have been obtained rigourously from the start by adding an external energy term − i q i zE(t)n i to the Lagrangian density in Eq. (5).
The idea here is to excite the dipole mode in order to separate the electron and positron populations using an oscillating dipolar electric field. Of course, this can always be achieved by using a sufficiently strong field, comparable to the electric field that binds the electrons and positrons together, which is of the order of 1 a.u. = 5.14 × 10 11 V/m.
One could hope to lower the required field by exciting the system at the resonant dipole frequency, i.e. E(t) = E 0 cos(ω 0 t) with ω 0 = Ω d . However, the effective potential between the electron and positron clouds is not harmonic, as is apparent from Fig. 2. As the distance d(t) grows under the influence of the resonant field, the system increases its energy and eventually reaches the anharmonic region of the confining potential. At this point, the laser frequency will no longer match the energy-dependent frequency of the confining potential, so that the resonance condition is lost and absorption of the laser light becomes inefficient. This is apparent from Fig. 3, where we show a numerical solution of the fully nonlinear equations (14)-(15) obtained with a second-order leap-frog method, for a cluster with N = 100. The distance d between the electron and positron clouds always remains much smaller than the size σ i of each cloud, so that no separation is achieved. We also note that the neutral breathing frequency Ω − is excited nonlinearly; indeed, both σ e and σ p oscillate in phase (they are indistinguishable on the figure) with a period very close to 2π/Ω − = 377 a.u. The above limitation can be overcome by resorting to autoresonant excitation [37]. Basically, autoresonance occurs when a classical nonlinear oscillator is externally excited by an oscillating field with slowly varying frequency. In our notation where E 0 is the excitation amplitude, g(t) is a Gaussian envelope function with peak value equal to unity, and ω 0 is equal to Ω d in our case. For instance, when α < 0, the time-dependent frequency ω(t) = ω 0 + α(t − t 0 ) is initially larger than the linear frequency, reaches ω 0 at t = t 0 , and then goes on slowly decreasing with a rate equal to α. It can be shown that, for |α| ≪ ω 2 0 and E 0 above a certain threshold, the instantaneous oscillator frequency becomes "locked" to the instantaneous excitation frequency, so that the resonance condition is always satisfied. In that case, the amplitude of the oscillations grows indefinitely and without saturation, until of course some other effect kicks in. It was previously shown [37] that the threshold behaves as E th 0 ∼ ω 1/2 0 |α| 3/4 , implying that the amplitude can be arbitrarily small, provided that the external frequency varies slowly enough. Autoresonant excitation has been been fruitfully applied to several systems, including charged antiparticles [38], the quantum pendulum [39,40], and semiconductor quantum dots [41].
We now apply an autoresonant excitation to an e-p cluster with N = 100, using the same amplitude as in Fig. 3, E 0 = 0.005 a.u., and α = −10 −4 (α must be negative because here the frequency is a decreasing function of the energy). The envelope function g(t) = exp[−(t−t 0 ) 2 /2∆ 2 ] reaches its maximum at t 0 , i.e., the instant when the external excitation frequency coincides with the linear resonant frequency Ω d . The width of the pulse is ∆ = 740 a.u. ≈ 18 fs, which is a realistic duration for current femtosecond laser pulses.
The results of fully nonlinear numerical calculations for the autoresonant case are shown in Fig. 4. With the same field amplitude as before, we observe that both the dipole d(t) (distance between the centers of the two clouds) and the widths σ i grow to very large values. This is a clear sign that the autoresonant technique allows us to go well beyond the resonant excitation at the linear frequency.
Nevertheless, since both d and σ i keep growing, the overlap between the electron and positron densities is not necessarily decreasing with time. This overlap is the truly interesting quantity, because it tells us whether the two populations are well separated or not. The overlap is also proportional to the probability of e-p annihilation, which we have neglected so far but may play a role over longer times. We define the normalized overlap I(t) as: where n 0 (r) is the initial equilibrium density, so that I(0) = 1. If we assume that σ e (t) = σ p (t) (which was observed to be true for all cases that we simulated), then the overlap can be computed analytically and yields: I(t) = exp(−d 2 /2σ 2 e,p ). In Fig. 5 we plot the time history of the overlap and of the "velocities"ḋ andσ i . The overlap quickly drops to zero as soon as the autoresonant mechanism becomes effective, which confirms that the two species do become well separated. The evolution ofḋ andσ i shows that both the dipole and the widths grow linearly in time but with different velocities, with |ḋ| > |σ i |.

Conclusion
In this paper we studied the static and dynamical properties of electron-positron clusters. Our model takes into account quantum and finite-size effects and incorporates the Coulomb forces and exchange interactions. Using a Lagrangian approach and a Gaussian ansatz for the density profiles, we were able to derive a set of ordinary differential equations for the radii of the electron and positron clouds, σ e and σ p , and the distance d between their centers of mass. The only free parameter of the model, β, related to the gradient correction to the exchange interaction, was determined by matching our static results with those issued from exact Hartree-Fock calculations.
The above approach allowed us to investigate for the first time the dynamical properties of e-p clusters. We first concentrated on the linear response, which revealed three resonant frequencies: a dipole mode (oscillations of d) and two breathing modes (oscillations of the σ i ). For typical parameters, these resonant frequencies lie within or near the visible spectrum.
The dipole mode can be excited with an external oscillating electric field, such as that provided by a laser pulse. However, the resonant excitation rapidly becomes inefficient when the dipole d grows beyond the harmonic part of the confining potential and starts exploring the nonlinear region, where the frequency is energy-dependent. This drawback was overcome by resorting to autoresonance, whereby the excitation frequency varies slowly during the pulse. The autoresonant technique allowed us to efficiently separate the electron and positron populations using a laser pulse in the visible -or near visible -range, with a peak electric field E 0 = 0.005 a.u. = 2.57 × 10 9 V/m and pulse duration ≈ 20 fs. These values are largely independent on the cluster size N and may be achieved experimentally using current ultrafast spectroscopy techniques.
Finally, it is important to stress that the species separation could be achieved in very short times (≈ 20 fs). This is much shorter than the lifetime of electrons and positrons in a positronium atom, which is 125 ps for the singlet state (para-positronium) and 142 ns for the triplet state (ortho-positronium). On the other hand, for a nonrelativistic e-p plasma (free e-p pairs), the annihilation rate γ D can be estimated with the Dirac formula [42,43]: γ D = πr 2 0 cn, where r 0 is the classical electron radius and c is the speed of light in vacuum. Using the typical density of an e-p cluster, this yields γ D ≈ 3 × 10 8 s −1 , i.e., an annihilation event roughly every 3 ns. This is again much longer than the separation time computed above.
Of course, the e-p lifetime in bound clusters like those studied in this work may differ significantly from the values observed for positronium atoms or for free e-p plasmas. Nevertheless, since the difference between the separation time and the estimated annihilation time is so large (4-5 orders of magnitudes), one may reasonably expect that effective separation can be achieved before annihilation starts playing a significant role.