text:
# On the equivalence of pairing correlations and intrinsic vortical currents in rotating nuclei

## I Introduction and methods

In a recent paper , some of the authors of the present article conjectured that the well-known dynamical effects of pairing correlations in rotating nuclei, explicited in the significant decrease of the moments of inertia from rigid body values , could be pictured as a vortical intrinsic collective motion coupled to the global rotation of the so-called intrinsic (body-fixed) frame relatively to the laboratory system. Namely, one intuitively would think of such a mode as aligned and counter-rotating with respect to the global angular velocity so as to produce the above quoted reduction of the moments of inertia. Furthermore, if one wants to preserve the nuclear shape in the presence of this intrinsic vortical mode the latter should be tangential to the grossly-defined nuclear surface.

As often suggested after R. Y. Cusson  (see also Ref.  for some relevant references), a simple ansatz for such a coupling is provided by the S-ellipsoid fluid dynamics as studied classically in great details by Chandrasekhar in the context of celestial self-gravitating objects . It makes use of a velocity field which is linear in the coordinates and produces a coupling of a global rotation with an intrinsic mode corresponding to a motion tangential to the body surface which is supposed to be ellipsoidal. Furthermore, the vorticities (i.e., the curl) of the two velocity fields are assumed to be aligned or anti-aligned along a principal axis of the ellipsoid. As shown, e.g., in Ref.  (see Fig. 1 there) such a simple parameterization is capable of yielding a variety of modes in the laboratory frame from rigid-body rotation to irrotational modes as well as various shear modes.

The connection of this S-ellipsoid fluid dynamics with pairing correlations may be discussed at various levels. First, one may note that a tangential intrinsic mode which is not changing the body shape (in the $`\stackrel{}{r}`$-sector of the phase-space) is only able to redistribute the density function in the $`\stackrel{}{p}`$-sector of the phase-space. Clearly, up to self-consistent couplings between the momentum redistribution and the body shape, the BCS pairing correlations behave exactly in the same way (redistributing the single-particle momenta). Second, in the semiclassical (Thomas-Fermi) calculations of Ref. , upon widely varying the pairing strength in some deformed rare-earth rotating nuclei, the authors have unambiguously exhibited current patterns very similar indeed to the classical current patterns of S-ellipsoids. Such currents had been shown before to exist in paired solutions at finite spin . They were however embedded in the midst of shell-effect—generated intrinsic currents whose existence had been found long ago (see the non-paired rotating solutions of ) and which will be further discussed in Section II A below.

However, the arguments above sketched are purely qualitative. It is the aim of this paper to investigate how quantitatively correct they might be.

To achieve this goal, it is first necessary to find a tool to quantify the intrinsic vortical content of a given current distribution. Within the context of S-ellipsoid fluids, it is rather natural to use the Kelvin circulation whose components are defined by (see, e.g., )

$$\widehat{K}_k=\underset{i,j}{}ϵ_{ijk}\left(\frac{a_j}{a_i}x_ip_j\frac{a_i}{a_j}x_jp_i\right)$$
(1)
where $`ϵ_{ijk}`$ is the totally antisymmetrical third rank tensor and the $`a_i`$’s are characteristic lengths associated with the ellipsoidal shape of the nuclei (e.g., the semi-axes). These lengths are proportional to $`\sqrt{x_i^2}`$. These operators clearly appear as the components of a doubly-stretched (in both $`\stackrel{}{r}`$ and $`\stackrel{}{p}`$) orbital angular momentum. Physically, they correspond to a rotation in the intrinsic frame of a sphere obtained by stretching the ellipsoidal intrinsic distribution. Indeed as shown by Lebowitz (see, e.g., Ref. ), the S-ellipsoid coupling mode may be understood as a global rotation along some principal axis $`\stackrel{}{u}`$, followed by the stretching to a sphere and the rotation (according to the same axis $`\stackrel{}{u}`$ of the considered ellipsoid) as above described, and finally by the inverse stretching back to the original ellipsoidal shape.

To check the validity of the above mentioned conjecture, we compare the results of two types of microscopic calculations, one involving pairing correlations in a purely rotating formalism, the other with no pairing but imposing the S-ellipsoid dynamics for the global rotation and intrinsic vortical mode velocity fields as well as for their coupling.

The pairing correlations are treated à la Hartree-Fock-Bogoliubov (HFB) using a Skyrme effective interaction in the particle-hole channel and a seniority force in the pairing channel within a code developed by H. Laftchiev . The triaxial character of the solutions arising from the time-reversal symmetry breaking inherent to rotating nuclei is treated within the Fourier decomposition formalism described in Ref. . Namely, the wavefunctions and relevant densities are defined in cylindrical coordinates on a mesh in the two spatial directions $`(r,z)`$ and as Fourier series in the azimuthal angle $`\theta `$. An approximate description of particle-number projection within the Lipkin-Nogami (LN) scheme  is also available in the code.

The S-ellipsoid coupling case is described within the Hartree-Fock-Skyrme generalized Routhian formalism of Refs. . This formalism makes it possible to describe the S-ellipsoid dynamics upon using a double constraint on both angular momentum and Kelvin circulation via two angular velocities $`\mathrm{\Omega }`$ and $`\omega `$ (relevant to each collective mode) along the first-axis.

These two microscopic approaches have been used in the following scheme. In a first step, we have used the HFB formalism with a constraint on angular momentum using realistic pairing strengths. We thus generate vortical currents in the intrinsic frame which can be quantified for instance by the mean value of the Kelvin circulation operator measured in the intrinsic frame and defined as

$$\widehat{K}_1=(\stackrel{~}{\stackrel{}{r}}\times \stackrel{~}{\stackrel{}{j}})_1\mathrm{d}^3r.$$
(2)
There $`\stackrel{~}{\stackrel{}{r}}`$ is the stretched position vector and $`\stackrel{~}{\stackrel{}{j}}`$ is the stretched current in the intrinsic frame whose components are obtained after subtracting the rotating-part to the usual current-density $`\stackrel{}{j}`$ as

$$\stackrel{~}{j}_i=a_i\left[\stackrel{}{j}\frac{m}{\mathrm{}}(\stackrel{}{\mathrm{\Omega }}\times \stackrel{}{r})\rho \right]_i.$$
(3)
where $`\stackrel{}{j}`$ is the usual quantum mechanical current density defined as

$$\stackrel{}{j}=\frac{1}{2i}(\stackrel{}{}_\stackrel{}{r}^{}\stackrel{}{}_\stackrel{}{r})\rho (\stackrel{}{r},\stackrel{}{r}^{})|_{\stackrel{}{r}=\stackrel{}{r}^{}}$$
(4)
The mean value $`I`$ of the angular momentum operator is defined through

$$\sqrt{I(I+1)}=\widehat{I}=(\stackrel{}{r}\times \stackrel{}{j}+\stackrel{}{\rho })_1\mathrm{d}^3r$$
(5)
where $`\stackrel{}{\rho }`$ is the spin-vector density.

In a second step, we perform pure Hartree-Fock calculations (i.e., without pairing) with a double constraint on both angular momentum and Kelvin circulation to obtain the same values $`I`$ and $`\widehat{K}_1`$ as in the preceding step. In the following, this type of calculations will be referred to as HF+V calculations.

The paper will be organized as follows. In the next Section we will present and compare the results of the numerical calculations yielded by the two formalisms, namely HFB and HF+V in several mass regions. Section III will be devoted to a summary of the salient features exhibited in Section II and to some perspectives for further studies.

## II Results

Numerical calculations have been performed for three heavy nuclei: two in the $`A150`$ and $`A190`$ superdeformation regions (namely for the 150Gd and the 192Hg nuclei) including LN approximate particle-number projection scheme and one very heavy nucleus (254No) without LN projection. The SkM\* parameterization  of the Skyrme effective force has been used. The pairing strengths $`G_q`$ for the charge state $`q`$ are given for each calculated nuclei in table I together with the corresponding single-particle configuration space in use. For all three nuclei, kinetic ($`𝔍^{(1)}`$) and dynamical ($`𝔍^{(2)}`$) moments of inertia are calculated using the formulae

$$𝔍^{(1)}=\frac{\widehat{I}}{\mathrm{\Omega }},$$
(6)
and

$$𝔍^{(2)}=\frac{\widehat{I}}{\mathrm{\Omega }}.$$
(7)
As was already mentioned in a previous paper , using these expressions for the moments of inertia, we assume that they can be computed as well from the energy-corrected LN wavefunctions. We also assume that the angular velocity dependence of the estimate of LN $`\lambda _2`$ parameter is small, otherwise $`𝔍^{(2)}`$ would have to be evaluated through second derivatives of the energy with respect to $`\mathrm{\Omega }`$.

### A The yrast superdeformed band of 150Gd

Up to now, fourteen superdeformed bands have been experimentally found in 150Gd . We will deal here with the yrast superdeformed band only. Using the Cranked-HF approach, it was impossible  to reproduce the trends of the experimental kinetic and dynamical moments of inertia of this band as can be seen in Fig. 1 (dotted line).

The results presented in the upper panel of Fig. 1 for the dynamical moment of inertia of this nucleus show a much better agreement with experimental data within the HFB+LN formalism as compared to the HF one. Indeed, we now reproduce the decrease of this moment of inertia along the band with theoretical values which are rather close to the experimental ones as was previously shown within pure HFB calculations by Bonche et al. . In the lower panel, the agreement of HFB+LN kinetic moment of inertia with the experimental data of Ref. , though much better than the HF one, is seen to be rather poor. However, the experimental value of this moment of inertia relies on an assumption of the lowest spin value of the observed band \[see eq. (6)\] which has been assumed in  to be $`I_{\mathrm{init}}=32\mathrm{}`$. Taking now for the initial spin of the band $`I_{\mathrm{init}}=34\mathrm{}`$ shifts up the experimental data and then yields a very good agreement with HFB+LN values of the moment of inertia.

Let us now focus on the intrinsic vortical content of the current density. In the upper part of Fig. 2 the Kelvin circulation mean value is plotted as a function of the angular velocity for 150Gd. One may first notice that the Kelvin circulation mean value in the HF case is not zero, even though the intrinsic vorticity constraint is absent. It is negative (with a rather small absolute value of about 1 $`\mathrm{}`$) for the protons while it is positive and larger (about 2 $`\mathrm{}`$) for the neutrons. These contributions should most likely be attributed to small inhomogeneities in the density due to shell effects. It is almost constant for the neutrons and slightly decreasing with $`\mathrm{\Omega }`$ for the protons. Dubbing this current patterns as shell effects may also be understood in rather general terms. In Ref.  within purely semiclassical microscopic calculations, it has been shown that for finite systems small surface-peaked counter-rotating intrinsic currents do appear. This effect is analogous to the diamagnetism which is observed in an electronic system submitted to a constant magnetic field. Thence the origin of all intrinsic currents (in Routhian HF calculations) which are not surface-peaked and counter-rotating (as obtained years ago in, e.g., ) are clearly and properly to be attributed to shell effects.

Comparing now the HF curves with the HFB ones, we see that the introduction of pairing correlations shifts down the Kelvin circulation mean value by about 2 $`\mathrm{}`$ for each charge state at low angular velocities and still 0.5 $`\mathrm{}`$ units at the end of the band. This decrease of the Kelvin circulation in the HFB case is related, as it has been already demonstrated in , to the fact that the collective effect of pairing correlations is equivalent to a counter-rotating intrinsic vortical motion (see also Ref. ). The combination of the above mentioned HF variation of $`\widehat{K}_1`$ with $`\mathrm{\Omega }`$ for protons and neutrons together with the regular decrease of the pairing correction as described previously yields an increase of the neutron contribution and an almost constant pattern for the protons. As a result the total Kelvin circulation exhibits an $`\mathrm{\Omega }`$ variation similar to the neutron one.

This decreasing pairing effect on the 150Gd Kelvin circulation is nicely correlated, as can be seen on the lower panel of Fig. 2, with the decrease of the pairing energy (defined as a trace of the product of the abnormal density with the gap potential). The latter is a well-known behavior, known as the Mottelson-Valatin effect , similar to the effect of a magnetic field on a superconductor below the critical point . Indeed as the pairing energy almost vanishes at high angular velocity, the Kelvin circulation in the HFB formalism almost reaches its HF value.

As said in the Introduction, the Kelvin circulation is a powerful tool to measure the intrinsic vortical content of the current distribution. In Fig. 3 (150Gd bottom part) we present the current patterns in the intrinsic frame of two rotational states of the HFB superdeformed yrast band of 150Gd, namely for $`I=32\mathrm{}`$ and $`I=46\mathrm{}`$ (corresponding to angular velocities $`\mathrm{\Omega }=381`$ keV and $`\mathrm{\Omega }=516`$ keV). In both cases, the neutron current patterns are rather disordered (due to shell structure effects). This is related to positive values of the neutron contribution to the Kelvin circulations. On the contrary, proton current patterns are rather well oriented, at least around the nuclear surface, along elliptic lines parallel to the surface, the more so in the $`I=32\mathrm{}`$ case where total Kelvin circulation and pairing energies are greater (in absolute value) than for $`I=46\mathrm{}`$. This trend is due to the fact (observed in Fig. 2) that the proton HF contribution to $`|\widehat{K}_1|`$ is very small. In the HF case (upper part), the current patterns show no particular order for both spins, and for protons as well as for neutrons. This clearly shows that the inclusion of pairing correlations favors the appearance of S-ellipsoid–like currents in rotating nuclei.

We then have performed HF+V calculations, namely constraining HF solutions to have $`I`$ and $`\widehat{K}_1`$ values \[see eqs. (2) and (5)\] identical to their HFB counterparts. The constraint in use is of an isoscalar type, that is constraining on the total (neutron + proton) angular momentum and Kelvin circulation. However, it can be seen in Fig. 2 that the proton and neutron contributions to the Kelvin circulation in HF+V calculations are very close to what is obtained in HFB calculations even though these contributions have not been constrained separately. Moreover, the current patterns for the HFB and HF+V approaches (see Fig. 3) are similar. Hence we have successfully grafted the superfluid part of the HFB currents into the HF calculations constraining merely the expectation values of the $`\widehat{K}_1`$ isoscalar operator.

The angular velocities within the three studied approaches, namely HF, HFB and HF+V, are plotted in Fig. 4 against the angular momentum. It is seen that the HFB and HF+V values fully overlap:

$$\mathrm{\Omega }_{HFB}\mathrm{\Omega }_{HF+V}.$$
(8)
It is worth noting that these two angular velocities are greater than $`\mathrm{\Omega }_{HF}`$. This is not surprising, since the counter-rotating collective effect of pairing has to be compensated in HFB calculations at a given angular momentum by an increase of the angular velocity. However, it is striking to see that the identity of the two momenta of the current-distribution, namely $`I`$ and $`\widehat{K}_1`$ in HFB and HF+V calculations, leads to identical rotational angular velocities $`\mathrm{\Omega }`$. It proves indeed, more quantitatively than the plots of Fig. 3, that the current patterns in the two approaches are similar everywhere, which could hardly be demonstrated by merely looking at integrated quantities.

The identity of the HFB and HF+V angular velocities seen in Fig. 4 yields an identity of the kinetic moments of inertia clearly demonstrated in eq. (6). As for the HF+V dynamical moments of inertia, even though they do not fully coincide with their HFB counterparts, they are indeed very close and even reproduce better the experimental data. In this nuclei, the HF+V and HFB formalisms yield solutions whose axial quadrupole moments remain constant over the entire band (differing by less than 2 percent around $`Q_{20}3900`$ fm2). Hence the equivalence of the dynamical moments of inertia in the two formalisms is clearly free from any deformation effect and reveals similar superfluid properties.

### B The yrast superdeformed band of 192Hg

We have used also the HF, HFB+LN and HF+V formalisms in order to describe the yrast superdeformed band of 192Hg and compare their results with the available experimental data .

As already noted by Gall et al.  the very low amount of pairing correlations attained above $`I30\mathrm{}`$ (as exemplified by the pairing energies for both neutrons and protons in the corresponding lower part of Fig. 2) makes it necessary to correct for the particle number symmetry breaking. This has been performed upon using the approximate Lipkin-Nogami particle number projection on top of the HFB calculation which has been found necessary to approach a proper description of the rotational band of this nuclei and avoid spurious drops in the moments of inertia due to sudden pairing-energy disappearance both for neutrons and protons. Whereas the LN prescription is known to be efficient in restoring some pairing correlations in weak-pairing regime as would be the case in HFB calculation for these nuclear states, the question of knowing whether or not it provides the correct amount of pairing correlations is still controversial (see, e.g., the discussion of Peru ). A more exact treatment of pairing correlations could be obtained for instance by generalizing to time-reversal symmetry breaking systems the method of Pillet et al.  which explicitly conserves the particle number.

The 192Hg dynamical moments of inertia are plotted in the upper panel of Fig. 1 within the three mentioned formalisms together with their experimental counterparts. As was already the case for 150Gd, the HF calculations fail to reproduce the experimental values as well as the behavior of this moment as a function of the angular velocity since it remains almost constant over the whole rotational band. By adding pairing correlations within the HFB+LN formalism, it is possible to reproduce the experimental values at low spins. However, using a simple seniority force as done here, one can only reproduce the increasing trend of the experimental moment of inertia and our theoretical value becomes to high at the upper bound of the rotational band. In view of some residual instabilities of the solution due to the presence of single-particle states in the vicinity of the cutoff energy defining the configuration space where pairing correlations are treated, we have introduced for this nucleus a cutoff factor of the Fermi-function type.

Similarly to what we have done for the 150Gd nucleus, we have attempted to model the collective effects of pairing correlations in this nucleus by performing HF+V calculations. As seen in Fig. 1, the calculated HF+V dynamical moments of inertia nicely match the HFB+LN ones, even better than for 150Gd. The axial deformation of the 192Hg nuclear states in the rotational band calculated within these two formalisms differ by less than 3 percent with an axial quadrupole moment mean value of $`Q_{20}44`$ barns.

In view of the rather low expectation values of the Kelvin circulation in the HF regime (see Fig. 2) its behavior as a function of the angular momentum in the HFB approach reflects mostly the interplay between the global rotation and the pairing correlations. One observes indeed a very interesting parabolic pattern. At low spins its low value is due to the fact that the intensity of the intrinsic currents is obviously related to the intensity of the global rotational currents which they are counteracting. On the other hand these intrinsic currents are, as we have assumed and demonstrated here, related with the intensity of pairing correlations. This balance explains the appearance of a maximum for $`|\widehat{K}_1|`$ at some “intermediate” spin value (around $`I20\mathrm{}`$). From the plotted value of $`\widehat{K}_1`$ for both charge states in Fig. 2, one may infer that as it was the case for 150Gd, the pairing generated intrinsic currents are almost entirely of an isoscalar character.

The agreement between HFB and HF+V dynamical moments of inertia clearly points into the direction of identical current patterns for both theoretical descriptions. This is indeed so, as shown in Fig. 3. Finally, one observes in Fig. 4 that the similarity of the kinetic moments of inertia in the two theoretical descriptions is reflected in the corresponding curves giving the spin-dependence of the angular velocity.

### C The ground-state deformed band of 254No

Recently, the existence of the rotational ground-state band of the heavy nucleus 254No has been assessed  by an observation of the $`\gamma `$-transition energies from spins $`I=4\mathrm{}`$ up to $`I=20\mathrm{}`$ in coincidence with the $`\alpha `$-decay chain of this nucleus. This high spin structure is the first one studied in $`Z>100`$ nuclei.

It has been shown by several authors  with different effective forces and pairing parameterizations that the properties of this isotope could be reproduced within HFB-type formalisms. The calculations presented here have been performed in the spirit of Ref. , that is fixing the pairing-force parameterization in order to reproduce the first transition energies. As a result, it happens that higher transition energies are fairly well reproduced. This procedure has been found in  to be more successful within a pure HFB than a HFB+LN formalism. Thus we will discuss here in terms of intrinsic vortical currents the results of our HFB calculations.

The dynamical moments of inertia obtained within the three formalisms HF, HFB and HF+V are plotted in the corresponding upper panel of Fig. 1 as functions of the angular velocity. As was the case for the two previously studied nuclei, it is still impossible here to reproduce the experimental data with pure HF calculations. On the contrary, the results of our HFB calculations fit rather well in the experimental patterns, and the HF+V results fully overlap the HFB ones. A slightly less good agreement between both theoretical estimates of $`𝔍^{(1)}`$ among themselves and with their experimental counterpart is obtained as exhibited in Fig. 1.

In Fig. 2, the HF and HFB estimates of $`\widehat{K}_1`$ for the 254No neutron current distributions are seen to be very close. This is clearly related as exhibited in the lower part of this figure to the very small amount of pairing correlations for this charge state at spin less than $`6\mathrm{}`$ (whereas for higher spin values, there are no correlated solutions). As for the protons now, the shell effect contribution to $`\widehat{K}_1`$ (i.e., the HF value) is quite significant, actually of the same order of magnitude as the pairing contribution (i.e., the difference between the HFB and the HF value). As a very specific feature of the HF plus vorticity calculations in this nucleus (versus those above reported in 150Gd and 192Hg) one notices that a purely isoscalar constraint on $`\widehat{K}_1`$ is not effective in reproducing the HFB results. Actually, whereas it is tempting to attribute such an isovector behavior for a $`Z100`$ nucleus (as opposed to $`Z80`$ or less nuclei) to the increasing relative importance of the Coulomb interaction (scaled as $`Z^2`$), the observed behavior of $`\delta |\widehat{K}_1|`$ (difference between HFB and HF+V values of $`|\widehat{K}_1|`$) which is positive for protons and negative for neutrons is yet to be understood (and should be assessed by systematic calculations).

In Fig. 3 intrinsic currents corresponding to two very different rotational regimes ($`I=2`$ and $`20\mathrm{}`$) are displayed for 254No. At low spins the neutron currents are small yet roughly consistent with a tangential pattern at the surface. This is due to the existence of a low pairing correlations regime. At high spin of course the HF and HFB neutron current patterns are coinciding and hardly tangential at the surface. This situation is quite at variance with what is obtained for the protons at large spins, where the prevalence of large pairing correlations entails the existence of rather well marked tangential counter-rotating currents.

Finally, the behavior of the angular velocity of both HFB and HF+V calculations on the one hand and of the HF calculations on the other hand, as a function of the total angular momentum is found quite similar here (see Fig. 4) to what was found before for 150Gd and 192Hg.

## III Conclusions

In the above, we have examined in a Routhian approach, the equivalence of the dynamical effects yielded by pairing correlations à la Bogoliubov with those obtained when a specific ansatz for the collective flow is imposed onto a purely Hartree-Fock type of microscopic calculations. This has been achieved in rather different contexts as far as the nucleon number ($`A`$ ranging from 150 to 250), the deformation (considering normally deformed bands as well as superdeformed bands), the range of spin values as well as the pairing correlation contents. In all cases however, we have found rather convincingly that constraining Hartree-Fock Routhian solutions to have the currents obtained by HFB Routhian calculations at the same spin values, was a sound substitute to these HFB calculations themselves.

For that purpose, the so-called Kelvin circulation operator has appeared as an efficient measure of the intrinsic currents. Since this operator is well suited to the collective dynamics dubbed as S-type ellipsoids by Chandrasekhar (implying a collective field which is linear in the coordinates), the actual HFB dynamical pattern should not deviate very much from this simple model pattern. That has been actually assessed by a mere inspection of current lines obtained in our HFB calculations in cases where at large enough spin values strong pairing correlations where present. Pairing correlations appear indeed as generating tangential counter-rotating currents in the intrinsic frame. This equivalence has been actually confirmed by a convergent array of evidences:

1. the similarity of dynamical moments of inertia,
2. the congruence of kinetic moment of inertia, or equivalently (at the same value of the angular momentum) of the associated angular velocity,
3. the correspondence between currents obtained in both approaches.

As for the currents, it has appeared that shell effects where visible in situations where either the collective intrinsic flow was too small (at low angular velocity for the global rotation) or when pairing correlations where vanishing. These shell effects generating density inhomogeneities are responsible for erratic current patterns (sometimes roughly similar to localized vortices) which are in general neither counter-rotating nor tangential.

As above discussed, upon increasing the spin value one experiences a competition between two processes. The first one is an increase of the intrinsic currents, a reaction effect which should be (and actually is, all things kept equal) increasing with the global rotation angular velocity. This type of currents decreases, as well known, the kinetic energy at a given angular velocity. The second one is the Mottelson-Valatin rotation-generated reduction of the pairing correlations on. The competition arises from the fact that this intrinsic flow should be decreasing with the decrease of the pairing correlations. One may thus easily understand the maximum obtained for the absolute value of the expectation value of the Kelvin circulation obtained for the 192Hg nucleus at intermediate spin values.

In our calculations an isoscalar constraint on the Kelvin circulation operator has been used. Whereas it has been deemed to be adequate for the two lighter nuclei under study, this has been hardly the case for the heavier one. Confirming and, if so, understanding this $`A`$-dependence or more likely $`Z`$-dependence, would require a rather systematic study which has not been performed here.

The present study is now extended in two directions.

At high spins the decrease of pairing correlations has often been extrapolated towards the appearance of a phase transition (associated thus with a complete disappearance of these correlations). It is clearly not very appropriate to quantitatively assess this feature in traditional approaches because one reaches there a low pairing regime where the Bogoliubov particle number symmetry-breaking approach is not very adequate there as well known. We are currently undertaking  Routhian calculations which explicitly conserve the particle number according to the general lines proposed in Ref. .

Even though we are representing here in a rather simple way the dynamical effects of pairing correlations at finite angular momentum as obtained in a Bogoliubov approach, we are still needing such lengthy calculations to yield the correct amount of intrinsic currents (through the expectation value of the Kelvin circulation operator). The challenge is to try to correlate for a given nucleus (i.e., for given nucleonic numbers, given deformations and given single particle level densities near the chemical potential at zero spin or equivalently the amount of pairing correlations at zero spin) the Kelvin circulation expectation value as a function of the angular velocity of the global rotation. This is currently achieved through a simple model approach inspired from the Mottelson-Valatin rotation anti-pairing effect and will be presented soon in an other paper .

## IV ACKNOWLEDGMENTS

Part of this work has been funded through an agreement (# 12533) between the BAS (Bulgaria) and the CNRS (France) and another (# 97-30) between the JINR (Russia) and the IN2P3/CNRS (France) which are gratefully acknowledged. One author (P. Q.) would like to thank the Theoretical Division of the LANL for the hospitality extended to him during his stay at Los Alamos. Two of the authors (H. L. and P. Q.) also wish to thank the INT (Seatle) for their kind invitation and the opportunity provided there to discuss some parts of this paper.


format: text

text:
# Charge and Salt Driven Reentrant Order-Disorder and Gas-Solid Transitions in Charged Colloids

## I INTRODUCTION

Monodisperse charge stabilized colloidal suspensions are studied with great interest because they exhibit a variety of structural orderings analogous to that occur in atomic systemsta01 ; ak98 . Unlike in atomic systems, the range and strength of the effective interparticle interaction among the constituent particles are tunable over a wide range in colloidal system. Hence, different structural orderings can be realized with ease in these systems at ambient conditions . Though the most dominant interparticle interaction among like-charged particles in a charge stabilized colloidal suspension is screened Coulomb repulsion, experimental observations in bulk suspensions such as reentrant transition ak88 ; ta96 , vapor-liquid transitionta96 ; ta92 , stable voids coexisting with ordered or disordered regions ise88 ; ito94 ; ta97 ; ta98  have served as evidences for the existence of a long-range attractive component in the effective interparticle interaction of like-charged colloidsta99 . Purely a repulsive pair-potential given by Derjaguin-Landau-Verwey-Overbeek(DLVO) theory der48  can not explain the inhomogeneous nature or the phase separation observed in like-charged colloidal suspensions. The DLVO theory is advanced by Sogami-Ise by considering the counterion mediated attraction and regarding the macroionic part as a one-component systemsog84 . The effective pair-potential $`U(r)`$ sog84  obtained from the Gibb’s free energy is found to have a long-range attractive term in addition to the usual screened Coulomb repulsive term given by DLVO theory. Tata et al’s computer simulations using this $`U(r)`$ could explain the above mentioned experimental observations satisfactorily ta01 ; ta96 ; ta90 ; ta93 ; ta96a ; ta98a .

Recently there have been several alternative theories van99 ; war00 ; chan01  to understand the observed phase separation in charged colloids. These theories, which are also known as ”volume term theories” van99 ; war00  have strong similarity in approach but differ in technical details and are the extensions of Debye -Huckel theory for asymmetric electrolytes. These theories consider the total Helmholtz free energy of the system that arises from Coulomb interaction between all species in the system. The outcome of all these theories is that pair-potential between colloidal particles can remain repulsive as a function of interparticle distance $`r`$, yet a gas-liquid or gas-solid phase separation can occur due to volume terms which are responsible for the spinodal instability. Recently Schmitzsch00  has brought out the connection between volume term theories and the long-range attractive term that is present in Sogami-Ise potential $`U(r)`$. It is shown that Sogami-Ise theory includes volume terms and the attractive term in $`U(r)`$ arises from the electrostatic interactions that are present in volume term theories. Further, the criticism on Sogami-Ise theory by Overbeek ob87  and Woodward’swood88  has been countered by Smalleysma90  and Schmitzsch93 . Hence, we use $`U(r)`$ as the model potential in our MC simulations.

The objective of this paper is two fold. (1) Carry out simulations as a function of $`\sigma `$ for volume fractions and salt concentrations which are the same as that used in experiments of Yamanaka et al yam98  and look for the reentrant solid-liquid phase behaviour observed by them . The reentered disordered state observed at high values of $`\sigma `$ is known to be inhomogeneous but the nature of the underlying phase transition has not been identified. (2) To identify conditions under which the reentered inhomogeneous disordered state will transform to a homogeneous state. Towards this end we perform a detailed MC study over a wide range of salt concentrations as $`C_s`$ is known to alter the position $`R_m`$ of the minimum as well as the depth $`U_m`$ of $`U(r)`$. Our simulations as a function of $`\sigma `$ show the reentrant order-disorder transition in agreement with experimental observations and we go on to predict a inhomogeneous to homogeneous transition as a function of $`C_s`$, that is yet to be observed in charged colloidal suspensions.

The details of MC simulations are given in the Sec. II. The effect of charge density on the phase behaviour of suspensions with different volume fractions at low salt concentrations are discussed in Sec. III, where the suspensions exhibit reentrant solid-liquid transition. Section IV includes the results on the influence of salt concentration on the reentered disordered state. Comparison of our simulation results with the results of reported experiments is presented in Section V. A brief summary with conclusions are given in section VI.

## II DETAILS OF SIMULATION

Monte Carlo (MC) simulations are carried out using Metropolis algorithm with periodic boundary conditions for a canonical ensemble (constant $`N`$, $`V`$, $`T`$ where $`N`$, $`V`$ and $`T`$ are respectively the number of particles, volume and temperature). Particles of diameter $`d`$(=120nm) are assumed to interact via a pair potential $`U(r)`$ having the functional form

$$U(r)=2\frac{(Ze)^2}{ϵ}\left(\frac{\mathrm{sinh}(\kappa d/2)}{\kappa d}\right)^2\left(\frac{A}{r}\kappa \right)\mathrm{exp}(\kappa r)$$
(1)
where $`A=2+\kappa d`$ coth($`\kappa d/2`$) and the inverse Debye screening length $`\kappa `$ is given as

$$\kappa ^2=4\pi e^2(n_pZ+C_s)/(ϵk_BT).$$
(2)
$`Ze`$ is effective charge on the particle (related to the surface charge density by $`\sigma =Ze/\pi d^2`$), $`C_s`$ the salt concentration, $`T`$ the temperature(298K), $`ϵ`$ the dielectric constant of water and $`k_B`$ is the Boltzmann constant. The position of the potential minimum $`R_m`$ is given as $`R_m=\{A+[A(A+4)]^{1/2}\}/2\kappa `$ and its depth by $`U_m=U(R_m)`$. Both $`R_m`$ and $`U_m`$ depend on $`\sigma `$ and $`C_s`$. For the required volume fraction $`\varphi `$ (= $`n_p\pi d^3/6`$), the length $`l`$ of the MC cell is fixed from the relation $`l^3=N/n_p`$. By performing simulation with different number of particles, N, we found in our earlier simulations ta96 ; ta93 ; ta98a  that the results are the same within the statistical error for $`N`$ 432. Hence N is fixed at 432.

The results presented in this paper correspond to $`\varphi `$=0.005 and 0.03. At these volume fractions the deionised suspensions are known to crystallize into a body centered cubic (bcc) ordering sir89 , hence we have chosen particles placed on a bcc lattice as the initial configuration in all our simulations. We monitor the total interaction energy $`U_T`$ and the first peak height of the structure factor, $`S_{max}`$ for identifying thermal equilibrium unambiguously bin79 ; ta92a . Most of the simulations away from the transition (e.g., freezing) took approximately 8 $`\times `$ 105 configurations to reach equilibrium, while those close to the transition needed nearly 6 $`\times `$ 106 configurations. A MC Step (MCS) is defined as $`N`$ attempted moves during which, on an average, each particle gets a chance to move. The step size to move the particles during MC evolution process is chosen in such a way that the trial acceptance ratio is always around 50$`\%`$. After reaching equilibrium pair-correlation function $`g(r)`$, coordinate averaged pair-correlation function $`g_c(r)`$ and $`S_{max}`$ are calculated using procedures reported earlier ta92a ; ta91 . $`g_c(r)`$ is obtained by averaging the coordinates of the particles over a sufficiently large number of configurationsta91 . Since $`g_c(r)`$ is free from thermal broadening and it helps in identifying the crystal structure unambiguously. Further, for suspensions which exhibit solid-like behaviour, the $`g_c(r)`$ shows sharper peaks as compared to the corresponding $`g(r)`$. The mean square displacement $`<r^2(t)>`$ for a chosen particle is defined as $`<r^2(m)>`$ = $`<𝐫_i(m+n)𝐫_i(n)^2>`$, where $`𝐫_i(m)`$ is the position of $`i`$ th particle after $`m`$ MCS and $`<\mathrm{}>`$ denotes the configurational average over the initial configurations $`n`$. Cluster size and its distribution are calculated as described in our earlier workta93 . The total fraction of particles participating in the clustering $`F_c`$, defined as the ratio of the total number of particles that participate in the clustering to the total number of particles $`N`$, has also been obtained for suspensions which exhibited inhomogeneous nature (phase separation). We identify the nature of suspension whether it is homogeneous or inhomogeneous by calculating the ratio $`d_s/d_0`$, where $`d_0`$ (= $`\sqrt{3}/2[1/n_p]^{1/3}`$) is the average interparticle separation calculated from the particle concentration $`n_p`$ and assuming a bcc-like coordination in the liquid-like order and $`d_s`$ is the average interparticle distance estimated from the first peak position in $`g(r)`$. For a homogeneous suspension $`d_s/d_0`$=1 and for an inhomogeneous suspension $`d_s/d_0`$ $`<`$ 1.

## III EFFECT OF CHARGE DENSITY

We have carried out MC simulations for two volume fractions $`\varphi `$ = 0.005 and $`\varphi `$ = 0.03 as a function of charge density on the particles keeping the salt concentration at 2$`\mu `$M. These parameters are same as that used by Yamanka et al.yam98  in their experiments. Fig. 1 shows the pair-correlation functions and projections of the corresponding MC cell onto a $`xy`$-plane for suspensions with $`\varphi `$=0.005 and $`\sigma `$=0.1, 0.23 and 0.5 $`\mu `$C/cm2. For low values of $`\sigma `$ (=0.1 $`\mu `$C/cm2), $`g(r)`$ and $`g_c(r)`$ showed a decay as a function of $`r`$, the corresponding particle positions(Fig. 1(B)) showed disorder and the ratio $`d_s/d_0`$ is found to be one (see table I). These observations suggest the suspension at this value of $`\sigma `$ is homogeneous and liquid-like ordered. On the other hand for suspension with $`\sigma `$= 0.27 $`\mu `$C/cm2 we find $`g(r)`$ does not decay with $`r`$ and the positions of sharp peaks in $`g_c(r)`$ clearly indicate that the ordering is bcc. The projection of particles in the MC cell shows perfect ordering and the calculated $`d_s/d_0`$ is found to be one confirming the bcc ordered crystalline structure is homogeneous. At high values of $`\sigma `$ ( = 0.5 $`\mu `$C/cm2) $`g(r)`$ and $`g_c(r)`$ are found to decay a function of $`r`$ once again . Further, we observe the first peak in $`g(r)`$ shifted to smaller $`r`$ implying $`d_s/d_0`$ $`<`$ 1. These observations suggest that suspension is inhomogeneous and disordered. The corresponding projected particles in the MC cell also confirm this and showed particle free regions (voids) coexisting with disordered dense region. Since the $`g_c(r)`$ shows sharper features than the corresponding $`g(r)`$, the dense disordered region is identified to be solid-like (amorphous).

In order to identify the values of $`\sigma `$ at which the homogeneous liquid (HL) freezes into a homogeneous crystalline (HC) order and then to an inhomogeneous disordered (also referred as a phase separated (PS)) state, the structural parameter $`S_{max}`$ is calculated as a function of $`\sigma `$ and is shown in Fig. 2. The sudden increase in $`S_{max}`$ at $`\sigma `$=0.23 $`\mu `$C/cm2 (Fig. 2(a)) corresponds to the freezing transition and the sudden drop in $`S_{max}`$ at $`\sigma `$=0.33 $`\mu `$C/cm2 corresponds to the transition from homogeneous crystalline state to a phase separated (gas-solid coexistence) state. Simulations for suspension with $`\varphi `$=0.03 showed similar phase behaviour except that the corresponding transitions have been found to occur at lower values of $`\sigma `$ (see Table I).

These observations are understood from the dependence of $`U(r)`$ on $`\sigma `$. At low charge density the suspensions remain homogeneous because the particles which are at a distance $`d_0`$ will always experience screened Coulomb repulsive interaction as $`d_0<R_m`$ (see curve $`a`$ of Fig. 3). The structural ordering in the homogeneous state depends upon the strength of $`U(r)`$ at $`r`$ = $`d_0`$. As the charge density on the particles is increased the strength of repulsive interaction increases. This results in freezing into a homogeneous crystalline state. Since $`R_m`$ decreases monotonically with increase in $`\sigma `$, particles which are at separation of $`d_0`$ will experience attraction when $`R_md_0`$ (see curve $`b`$). Due to strong attraction the particles condense into a solid-like dense phase leaving some particles in the rare phase (gas-like) which appear as voids when the fraction of the volume occupied by the rare phase is smaller than that of the dense phase. Thus the occurrence of freezing of a homogeneous liquid into a homogeneous crystalline state at lower values of $`\sigma `$ and the system reentering into a disordered inhomogeneous phase (gas-solid coexistence) at higher values of $`\sigma `$ are understandable.

## IV EFFECT OF SALT ON INHOMOGENEOUS DISORDERED PHASE

There have been reports of theoretical ta98a ; ta91 ; rob88 ; sen91  and experimental investigationsak98 ; sir89 ; mon89  on the effect of $`C_s`$ on melting/freezing of a homogeneous colloidal suspension. However, the conditions under which an inhomogeneous disordered state, which occurs in highly charged colloids, becomes homogeneous are not known. When the well depth $`U_m`$ $`<`$ $`k_BT`$, the attraction between particles becomes weak, hence the suspension can remain homogeneous. It can be seen from Fig. 3 (curve $`d`$) that the well depth $`U_m`$ decreases to values smaller than $`k_BT`$ at high salt concentrations. So, we performed MC simulations over a wide range of salt concentrations keeping the other suspension parameters fixed and results are summarized in Table II. Fig. 4 shows the effect of increasing salt concentration on the disordered dense phase which coexisted with voids at $`C_s`$ = 2$`\mu `$M (Fig. 1(c)). For low values of $`C_s`$, $`g(r)`$ shows a decay, but the structural correlations persist to the full length of the MC cell (Fig. 1(c)) suggesting that the disordered dense phase is connected. The fraction of volume occupied by the voids is found to be smaller than that occupied by the dense phase. Hence the voids constitute the minority phase at low values of $`C_s`$. The particle positions in the MC cell(see Fig. 1(c)) confirm this. When the salt concentration is increased we observe structural correlations in $`g(r)`$ (Fig. 4(a)) to extend up to a few times the diameter of the particle. The first peak height is also found to be very high and the first peak position matches exactly with $`R_m`$. These observations imply the formation of dense phase clusters with $`R_m`$ as the average distance between particles within the cluster. The other sharp peaks appearing at distances beyond the first peak corresponds to short-range ordering of particles within the clusters (intracluster ordering).

Since the clustering occurs at large values of $`C_s`$, the ordering of particles within a cluster could be liquid-like or solid-like (amorphous). In order to distinguish between the two types of ordering we have calculated $`g_c(r)`$, g(r) and the mean square displacement $`<r^2(t)>`$ of particle within the cluster. We compare the behaviour of these quantities with that corresponding to a vapor-liquid coexistence state. This vapor-liquid coexistence state obtained by performing simulations for suspension parameters at which a vapor-liquid coexistence has been reported ta01 ; ta96 ; ta92 . It can be seen from Fig. 5 that the peaks in $`g_c(r)`$ (curve $`a`$) corresponding to $`C_s`$=550$`\mu `$M are much sharper with increased peak height as compared to the corresponding $`g(r)`$ (Fig. 4(a)). Whereas the $`g_c(r)`$ corresponding to the vapor-liquid coexistence (dotted curve in Fig. 5) is much lower in peak height and also poorer in features as compared to the corresponding g(r) (inset of Fig. 5). The mean square displacement of particles inside the cluster corresponding to $`C_s`$=550$`\mu `$M shows a saturation behaviour as a function of MC time and undergo much smaller displacement as compared to that of particle within a liquid-like ordered cluster. These observations unambiguously suggest that the clusters formed at higher values of $`C_s`$ for highly charged colloids are solid-like. These clusters can be observed either in the scattering or microscopy experiments provided the density of solvent is matched with the density of colloidal particles. Tata et al have observed the coexistence of voids with dense phase disordered regions ta97 ; ta98  and existence of clusters under density matched conditions ta98  in highly charged polycholorostyrene sulfonate particles (PCSS) dispersed in an aqueous medium. Density matching for PCSS particles has been achieved by redispersing them in 60% glycerol aqueous solution.

Simulations at a salt concentration of 1000 $`\mu `$M showed only a single peak in g(r) (see Fig. 4(b)) at $`r=R_m`$. The corresponding projection of particles in the MC cell (see inset of Fig. 4 (b)) do not show clustering and are distributed uniformly. These observations suggest that the suspension with $`C_s`$= 1000 $`\mu `$M is homogeneous and noninteracting (gas-like). Thus we observe for $`C_s`$ $``$ 550 $`\mu `$M an inhomogeneous state in the form of dense phase clusters having solid-like ordering and a homogeneous gas-like disordered state at $`C_s`$=1000$`\mu `$M. Hence the inhomogeneous (PS) to homogenous disordered transition is expected to occur in between these salt concentrations. We identify this transition by monitoring the total fraction of particles participating in the clustering $`F_c`$ (Fig. 6(a)) and the first peak height $`g_{max}`$ of $`g(r)`$ (Fig. 6 (b)). Note the sudden change in both of these parameters around $`C_s`$ = 570$`\mu `$M as $`C_s`$ is increased. This change is associated with the inhomogeneous (PS) to homogenous (HG) transition or in other words transition from a gas-solid coexistence to a gas phase. The large values of $`g_{max}`$ in the inhomogeneous state arises due to scaling of $`g(r)`$ with respect to particle concentration $`n_p`$. The dense phase concentration $`n_d`$ is estimated from the first peak position ( $``$ $`R_m`$) of $`g(r)`$ and is found to be several times higher than $`n_p`$ (see Table II). In the case of suspension with $`\varphi `$=0.03 the transition is found to occur at $`C_s`$=270$`\mu `$M which is lower than $`C_s`$ = 570$`\mu `$M for $`\varphi `$ = 0.005. Suspensions of higher volume fraction undergo PS to HG transition at lower values of $`C_s`$. This is due to increased screening of particles in suspensions of higher volume fraction as compared to that in suspensions with lower $`\varphi `$.

## V COMPARISON WITH EXPERIMENTS

It is important to compare our simulations results with that observed experimentally. Our simulations as a function of $`\sigma `$ revealed the existence of the three phases viz., a homogenous liquid (HL) at low $`\sigma `$, a homogenous crystalline (HC) state at intermediate values of $`\sigma `$ and an inhomogeneous disordered state (PS) at high values of $`\sigma `$. Further, we showed that the homogenous liquid freezes into a homogeneous crystalline state at about $`\sigma `$ $``$ 0.22 $`\mu `$C/cm2 and disorders once again at about $`\sigma `$ $`ge`$ 0.33 $`\mu `$C/cm2. Yamanaka et al.yam98  have reported freezing of a homogenous liquid-like ordered aqueous suspension of charged particles with $`\varphi `$=0.005 into a homogenous crystalline state for charge densities beyond 0.24 $`\mu `$C/cm2. The suspension is found to disorder once again when $`\sigma `$ is increased beyond 0.40 $`\mu `$C/cm2. This reentered disordered state is identified to be inhomogeneous from the scattering measurements. The measured average interparticle separation $`d_{exp}`$ is found to be less than $`d_0`$yam98 ; yam96 , which implies the nonspacefilling nature of the suspension at high charge densities. Thus present simulations are in good agreement with the experimental observations. Further, our simulations at higher volume fraction have revealed that the values of $`\sigma `$ at which HL to HC transition and HC to PS transition occur are lower as compared to that in suspensions of lower $`\varphi `$. This result is also in agreement with experimental observationsyam98 .

Present simulations on the reentrant inhomogeneous disordered state as a function of salt concentration revealed that the inhomogeneous state undergoes a transition to homogeneous gas phase beyond a critical salt concentration. This gas-solid transition which can be observed in highly charged colloids by varying salt concentration, has not yet investigated experimentally. Further, we notice from these simulations the absence of vapor-liquid coexistence for highly charged colloids. In this context, it is worth mentioning the result from MC simulations on a system of hard particles with an attractive Yukawa interaction by Hagen and Frankel ha94 . They showed the disappearance of vapor-liquid coexistence when the range of attractive part of the Yukawa potential is less than approximately one sixth of hard-core diameter. In present simulations the range of attractive and repulsive parts of $`U(r)`$ have been varied by changing $`\sigma `$ and $`C_s`$. The range of attraction is found to be smaller for suspensions of high charge density particles as compared to the suspensions of low charge density particles. Hence, one can observe a vapor-liquid coexistence in suspension of low charge density particles ta96 ; ta92  but not in suspensions of highly charged particles.

## VI CONCLUSIONS

MC simulations using a pair-potential with long-range attractive term show a homogenous liquid, homogenous crystal and a disordered inhomogeneous phase as the charge density is varied, suggesting a reentrant order-disorder transition, which is in agreement with the observed experimental results. The reentered inhomogeneous phase is identified as the gas-solid coexistence state and occurs due to strong attraction experienced by the highly charged particles. Suspensions exhibit homogeneous liquid-like and homogenous crystalline orders at lower values of $`\sigma `$ due to the dominant screened Coulomb repulsion at average interparticle separation. The inhomogeneous disordered state undergoes a transition to a homogenous gas state upon increasing the salt concentration. This transition is analogous to the gas-solid transition observed in atomic and colloid-polymer mixtures as a function of temperature and free polymer, respectivelyile95 . We believe that our prediction of gas-solid transition driven by the salt concentration would motivate experimentalists to investigate highly charged colloids in search of this transition as well as for detailed phase diagrams with a gas-liquid or gas-solid coexisting phases.

Authors wish to thank M.C. Valsakumar for valuable discussions and helpful suggestions.

Author to whom correspondence should be addressed.


text:
# Exponential velocity tails in a driven inelastic Maxwell model

## Abstract

The problem of the steady-state velocity distribution in a driven inelastic Maxwell model of shaken granular material is revisited. Numerical solution of the master equation and analytical arguments show that the model has bilateral exponential velocity tails ($`P(v)e^{|v|/\sqrt{D}}`$), where $`D`$ is the amplitude of the noise. Previous study of this model predicted Gaussian tails ($`P(v)e^{av^2}`$).

Recently, granular systems have been intensively studied SWINNEY ; MUZZIO . One reason of such an interest is the fact that in many respects these systems are highly unconventional and even exotic KADANOFF . As a result, very often even standard laws of statistical physics require certain modification when applied to granular systems.

As an example, let us consider the Maxwell-Boltzmann law for velocity distribution $`P(v)`$ of atoms or molecules in a gaseous state, which states that $`P(v)\mathrm{e}^{av^2}`$. Under certain experimental conditions, a granular system can be considered as a gas and a natural question is what is its velocity distribution. Numerous theoretical and experimental studies do not provide a simple answer to this question. On the contrary, they show that $`P(v)`$ depends on certain details of the system as for example how energy is transferred by a thermostat into the system, in order to balance the energy loss during inelastic collisions. Theoretical works have shown that $`P(v)`$ might be of the form $`\mathrm{e}^{av^2}`$ SWINNEY1 ; BIBEN , $`\mathrm{e}^{av^{3/2}}`$ NOIJE ; BARRAT  or $`\mathrm{e}^{av}`$ MONTANERO . Under certain conditions experiments show clear deviations from the Maxwell-Boltzmann law, but it is still rather difficult to decide what is the form of $`P(v)`$ in real granular systems OLAFSEN .

Very often granular systems are described using the so-called inelastic hard-sphere models, which could be then analysed using corresponding Boltzmann equations. These off-lattice, two- or three-dimensional systems are however very difficult to study, especially when we want to explore large-velocity regions of the phase space. A possible alternative is to construct simplified models for which more accurate calculations are possible. One class of such models are Maxwell models, for which the collision term in the corresponding Boltzmann equation is velocity independent. Recent studies show that for Maxwell models $`P(v)`$ might take either exponential ($`\mathrm{e}^{av}`$ERNST  or Gaussian ($`\mathrm{e}^{av^2}`$) form BEN ; CARILLO .

A particularly interesting Maxwell model was proposed by Ben-Naim and Krapivsky (BN-K) BEN . These authors presented an elegant solution of the master equation of their model using a Fourier transform method. However, their conclusion that $`P(v)`$ has a Gaussian decay is based on certain approximation whose validity is difficult to assess. Namely, they infer the large-velocity behaviour of $`P(v)`$ from the low-$`k`$ behaviour of the Fourier transform $`\widehat{P}(k)`$. Since BN-K model is one of the very few models for which exact or numerical but very precise calculations can be made, it would be desirable to clarify the validity of this approach.

In the present paper we reexamine the BN-K model. Analysing numerically the solution of the master equation of the model we obtain the velocity distribution $`P(v)`$. Asymptotically ($`v\mathrm{}`$), this quantity shows bilateral exponential tails ($`P(v)\mathrm{e}^{a|v|}`$) and such a behaviour is seen over more than 10 decades. We also present analytical arguments which support existence of bilateral exponential tails in this model. From our analysis it follows that $`a=\frac{1}{\sqrt{D}}`$, where $`D`$ is the amplitude of noise which simulates the input of energy into the system.

To introduce the model, we consider a collection of particles that are characterized by a single parameter, their velocity $`v`$. In this model positions of particles are not specified hence the model neglects any spatial correlations. Particles undergo two-body inelastic collisions that change their velocities according to $`(v_1,v_2)(v_1^{},v_2^{})`$, where

$$\left(\begin{array}{c}v_1^{}\\ v_2^{}\end{array}\right)=\left(\begin{array}{cc}\gamma & 1\gamma \\ 1\gamma & \gamma \end{array}\right)\left(\begin{array}{c}v_1\\ v_2\end{array}\right),$$
(1)
and $`\gamma `$ is the inelasticity parameter ($`0<\gamma <1`$). Particles which participate in a collision are chosen at random. In addition to the collision rule (1), particles are subjected to the uncorrelated white noise of strength $`D`$. Existence of noise ensures that the model has a well-defined nontrivial steady state. In the steady-state the velocity distribution $`P(v)`$ satisfies the following equation BEN

$$DP^{\prime \prime }(v)=P(v)+\frac{1}{1\gamma }_{\mathrm{}}^+\mathrm{}𝑑uP(u)P\left(\frac{v\gamma u}{1\gamma }\right),$$
(2)
where double prime <sup>′′</sup> denotes the second derivative with respect to velocity. To solve eq. (2) BN-K introduced the Fourier transform of the velocity distribution $`\widehat{P}(k)=𝑑ve^{ikv}P(v)`$, which in the steady state satisfies

$$(1+Dk^2)\widehat{P}(k)=\widehat{P}[\gamma k]\widehat{P}[(1\gamma )k].$$
(3)
This equation admits the following solution

$$\widehat{P}(k)=\underset{i=0}{\overset{\mathrm{}}{}}\underset{j=0}{\overset{i}{}}[1+\gamma ^{2j}(1\gamma )^{2(ij)}Dk^2].$$
(4)
After some transformations eq. (4) can be written as

$$\widehat{P}(k)=\mathrm{exp}\left[\underset{n=1}{\overset{\mathrm{}}{}}\frac{(Dk^2)^n}{na_{2n}(\gamma )}\right],$$
(5)
where $`a_n(\gamma )=1(1\gamma )^n\gamma ^n`$. To obtain a large velocity behaviour of $`P(v)`$ BN-K truncate the series in eq. (5) keeping only the first, quadratic in $`k`$, term

$$\widehat{P}(k)\mathrm{exp}\left[\frac{Dk^2}{a_2(\gamma )}\right].$$
(6)
Subsequently, taking the inverse Fourier transform of this Gaussian function they obtain $`P(v)\mathrm{exp}(\frac{v^2}{2v_0^2})`$, where $`v_0^2=2D/a_2(\gamma )`$.

In our opinion, the above procedure of truncating the series and inverting the resulting function is not well justified. The approximation (6) agrees with the exact expression (5) only up to the $`k^2`$ term of the Taylor expansion. Equally consistent approximation can be given as

$$\widehat{P}(k)\frac{c^2}{c^2+k^2},$$
(7)
where $`c=\sqrt{\frac{a_2(\gamma )}{D}}`$. But the inverse Fourier transform of (7) is $`\frac{c}{2}\mathrm{exp}(c|v|)`$, which is qualitatively different than the Gaussian decay obtained by BN-K. Both (6) and (7) are consistent with the exact solution (5) up to the $`k^2`$ order but differ at the higher orders. Apparently, these higher order terms qualitatively affect the large-velocity limit of $`P(v)`$. Consequently, without more detailed arguments such approximations are not justified. Indeed, as we will show below, $`P(v)`$ has a bilateral exponential decay but with a different coefficient $`c`$.

To check the validity of the BN-K approach we calculate $`P(v)`$ numerically as an inverse Fourier transform of $`\widehat{P}(k)`$. One way to compute $`\widehat{P}(k)`$ is to evaluate the infinite sums in the logarithm of eq. (4). However, it turns out that using eq. (5) leads to a better precision. Let us notice, however, that the series in eq. (5) converges but only for $`|k|<k_c=\frac{1}{\sqrt{D}}`$. To calculate $`\widehat{P}(k)`$ for $`|k|>k_c`$ we can then use eq. (3), provided that $`\gamma k`$ and $`(1\gamma )k`$ fall within the range of convergence. If not, we have to calculate $`\widehat{P}(\gamma k)`$ and $`\widehat{P}((1\gamma )k)`$ referring again to eq. (3). Implementing this recursive procedure, we calculated $`\widehat{P}(k)`$ and then using the Fast Fourier Transformation algorithm we obtained $`P(v)`$. Our results for $`D=1`$ are shown in Fig. 1. This figure clearly shows that $`P(v)`$ decays exponentially. For $`\gamma =0.5`$ and 0.8 such a behaviour is seen for more than 10 decades. For $`\gamma `$ close to unity there is a $`\mathrm{exp}(av^2)`$ decay for small velocities which asymptotically is replaced by the bilateral exponential decay ($`\mathrm{exp}(c|v|)`$). Moreover, as shown in Fig. 1 the asymptotic slope $`c`$ is independent of $`\gamma `$ and is approximately one. As we will show below this slope only depends on the amplitude of noise $`D`$.

Before doing that let us notice that the $`D`$-dependence of the velocity distribution $`P(v)`$ can be easily inferred from the fact that $`D`$ and $`k`$ enter its Fourier transform only through $`Dk^2`$ terms (see eqs. (4) and (5)). From this property one can easily obtain that

$$P(v,D)=\frac{1}{\sqrt{D}}P\left(\frac{v}{\sqrt{D}},D=1\right),$$
(8)
where we explicitly indicated the dependence on the noise $`D`$.

In the following we provide some analytical arguments showing that $`P(v)`$ decays exponentially. Let us assume that in eq. (2) the second term (gain) can be neglected. Then eq. (2) simplifies to $`DP^{\prime \prime }(v)=P(v)`$ and the normalizable solution reads

$$P(v)e^{\frac{|v|}{\sqrt{D}}}.$$
(9)
Let us notice that the exponential part is $`\gamma `$-independent (since $`\gamma `$ enters the master equation only through the neglected gain term). Moreover, such a solution is in a very good agreement with numerical calculations (see Fig. 1).

For more general models it is known that neglecting the gain term is justified for $`v\mathrm{}`$ when the resulting solution decays faster than exponentially NOIJE ; BARRAT . Our solution (9) is thus a marginal case. However, we can show that for our model the gain term in the limit $`v\mathrm{}`$ indeed can be neglected. First, let us evaluate the gain term for the solution (9). Elementary integration for $`D=1`$ and $`v>0`$ gives

$$\begin{array}{c}\frac{1}{1\gamma }_{\mathrm{}}^+\mathrm{}𝑑u\mathrm{exp}(|u|)\mathrm{exp}\left[\left|\frac{v\gamma u}{1\gamma }\right|\right]\hfill \\ \hfill =\frac{2}{2\gamma 1}\left[\gamma \mathrm{exp}\left(\frac{v}{\gamma }\right)+(\gamma 1)\mathrm{exp}\left(\frac{v}{1\gamma }\right)\right]\end{array}$$
(10)
One can see that since $`0<\gamma <1`$, the gain term decays exponentially with $`v`$ but faster than the solution (9) (i.e., a non-neglected loss term in (2)). We expect that eq. (9) is only an asymptotic ($`v\mathrm{}`$) solution of the master equation (2). Thus, for velocities $`uO(1)`$ the distribution $`P(u)`$ deviates from the asymptotic form (9). This will modify the integral (10) but only in the vicinity of $`u=0`$ and $`u=\frac{v}{\gamma }`$. As we argue below, such a modification of $`P(u)`$ still leads to the gain term decaying faster than the loss term.

Indeed, the contribution around $`u=0`$ is a product of $`P(u)`$ and of the exponential term $`\mathrm{exp}\left(\frac{v\gamma u}{1\gamma }\right)`$. Thus, for large $`v`$ a modification of $`P(u)`$ around $`u=0`$ will change our estimation (10) but only at the order of $`\mathrm{exp}(\frac{v}{1\gamma })`$. Similarly, one can show that a modification of $`P(u)`$ around $`u=\frac{v}{\gamma }`$ will change (10) by a factor of the order of $`\mathrm{exp}(\frac{v}{\gamma })`$. Consequently, the gain term again decays faster than the solution $`P(v)`$, which justifies its neglect. Let us also notice that when $`\gamma `$ approaches 0 or 1 the model becomes energy conserving and the bilateral exponential distribution (9) is no longer expected to hold. But it is easy to notice that in such cases the decay of the gain term matches the decay of the solution $`P(v)`$ and it cannot be neglected.

As a comment let us notice that knowing the second moment of $`P(v)`$ BEN  enables us to calculate the effective temperature $`T`$ of our systems defined as an averaged square velocity

$$T=v^2=\frac{D}{\gamma (1\gamma )}.$$
(11)
To check the validity of our calculations we also made Monte Carlo simulations of this model. Results, which are shown in Fig. 2 confirm the bilateral exponential decay of $`P(v)`$ although the accuracy is this time much lower. Our Monte Carlo data are rescaled in such a way that their variance matches that obtained using the Fourier transform method and a very good agreement is seen even on the logarithmic scale.

In conclusion, we have shown that the Maxwell model proposed by Ben-Naim and Krapivsky has velocity distribution decaying as a bilateral exponential. Together with the recent results by Ernst and Brito ERNST  it indicates that such a decay might be of more generic nature for this class of models. As a possible extension, it would be desirable to examine some other Maxwell models where velocities are not scalars but rather $`d`$-dimensional vectors. Actually, such models were already studied and the analysis indicates that, for increasing $`d`$, correlations between velocities and deviation from the pure Gaussian distribution decrease BEN1 . One possibility is that there might be a critical dimension $`d_c`$ and such that for $`d<d_c`$ the velocity distribution has a bilateral exponential decay (as in the present model) while it it has Gaussian decay for $`dd_c`$. Analysis of such models is however left for the future.

###### Acknowledgements.

We acknowledge interesting discussion with François Coppex. This work was partially supported by the Swiss National Science Foundation and the project OFES 00-0578 ”COSYC OF SENS” and the Hungarian Academy of Sciences (Grant No. OTKA T029792).