text:
# Impurity in a Maxwellian Unforced Granular Fluid

## I Introduction

Granular media are typically polydisperse. For example, sand and grains have a broad range of particle sizes and shapes. Such granular mixtures exhibit size segregation, a ubiquitous collective phenomena that underlies diverse processes including for example production and transport of powders in industry, sand dunes propagation and volcanic flows in geophysics . The “Brazil Nut” problem where an impurity is immersed in a uniform granular media is an extreme realization of a granular mixture as it corresponds to the vanishing volume fraction limit of a binary mixture. While this problem has been extensively studied, dynamics of such impurities are not fully understood .

Understanding the velocity statistics of granular mixtures is a necessary step in describing multiphase granular flows. Experimental and theoretical studies show that in general, different components of a granular mixture are characterized by different typical speeds, i.e., granular temperatures are usually distinct . However, analytical treatment of this case is difficult given the coupling between the different components of the mixture. An additional complication arises from the large number of parameters, including several restitution coefficients, governing the dynamics.

In contrast, the impurity problem is more amenable to analytical treatment. First, the impurity is directly enslaved to the fluid and second, there are fewer parameters . In this paper, we study the impurity problem in the framework of the Maxwell model that assumes that the collision rate is uniform, i.e., independent of the relative velocity of the colliding particles. This simplifies the Boltzmann collision operator and thence this model is widely used in kinetic theory. The Maxwell model is analytically tractable even when the collisions are inelastic as shown in a number of recent studies of uniform and polydisperse granular gases .

We obtain analytic results for the velocity distributions valid for arbitrary spatial dimension and collision parameters. We consider two versions of the Maxwell model. In the first, termed the Inelastic Maxwell Model (IMM), the collision rates are completely independent of the relative velocities of the colliding particles. In the second, termed the Improved Inelastic Maxwell Model (IIMM), the collision rates are proportional to the average relative velocity of the colliding particles.

In the IMM, there are two phases separated by a critical impurity mass $`m_{}(r_p,r_q)`$, determined by the restitution coefficients $`r_p`$ and $`r_q`$ characterizing fluid-fluid and fluid-impurity particle collisions, respectively. When the impurity mass $`m`$ is smaller than the critical mass, $`m<m_{}`$, different temperatures characterize the fluid and the impurity. The impurity velocity distribution exhibits similar characteristics as does the fluid; in particular, the same exponent governs the algebraic decay of the large velocity tail. Asymptotically, the impurity to fluid temperature ratio diverges as the critical mass is approached and is infinite when $`mm_{}`$. In this regime, the impurity essentially scatters of an ensemble of static fluid particles. The governing Lorentz-Boltzmann equation simplifies considerably. In general, moments of the velocity distribution exhibit multiscaling asymptotic behavior and the velocity distribution consists of replicas of the initial conditions. On the other hand, this phase transition is suppressed in the IIMM. The temperature ratio remains finite and the system is always in the light impurity phase. Nevertheless, secondary transitions affecting large velocity moments remain.

The rest of this paper is organized as follows. In Sec. II, we describe the IMM and present the Boltzmann equation for the fluid velocity distribution and the Lorentz-Boltzmann equation for the impurity velocity distribution. We then determine the temperature and show that a phase transition occurs. The light impurity phase is analyzed first and then the heavy impurity phase, emphasizing the behavior of the velocity distribution and its moments. In Sec. III, we discuss the IIMM and show that the system always remains in the light impurity phase. We conclude with a summary.

## II The Inelastic Maxwell Model

We study dynamics of a single impurity particle in a uniform background of identical inelastic spheres. Without loss of generality, we set the mass of the fluid particles to unity, while the mass of the impurity is denoted by $`m`$. Particles interact via binary collisions that lead to exchange of momentum along the impact direction. Collisions between two fluid particles are characterized by the collision parameter $`p`$, and collisions between the impurity and any fluid particle are characterized by the collision parameter $`q`$. When a particle of velocity $`𝐮_1`$ collides with a fluid particle of velocity $`𝐮_2`$ its post-collision velocity $`𝐯_1`$ is given by

$`𝐯_1`$ $`=`$ $`𝐮_1(1p)(𝐠𝐧)𝐧,`$ (1)
$`𝐯_1`$ $`=`$ $`𝐮_1(1q)(𝐠𝐧)𝐧,`$ (2)
with $`𝐠=𝐮_1𝐮_2`$ the relative velocity and $`𝐧`$ the unit vector parallel to the impact direction. The first equation gives the velocity of a fluid particle and the second that of an impurity. The collision rules (1)–(2) are derived by employing momentum conservation combined with the fact that in an inelastic collision, the component of the relative velocity parallel to the impact direction is reduced by a factor equal to the fluid-fluid (impurity-fluid) restitution coefficient $`r_p`$ ($`r_q`$). The restitution coefficients are related to the collision parameters via

$`r_p`$ $`=`$ $`12p,`$ (3)
$`r_q`$ $`=`$ $`m(m+1)q.`$ (4)
Since the restitution coefficients obey $`0r1`$, the collision parameters satisfy $`0p1/2`$ and $`\frac{m1}{m+1}q\frac{m}{m+1}`$. The energy dissipated in each collision is equal to $`A(𝐠𝐧)^2`$, with $`A=2p(1p)`$ and $`A=m(1q)[2(m+1)(1q)]`$ for fluid-fluid and impurity-fluid collisions, respectively.

We consider a collision process where random pairs of particles undergo inelastic collisions with a random impact direction. This inelastic Maxwell model (IMM) is described by a Boltzmann equation with a uniform collision rate. The model is a straightforward generalization of the classical Maxwell model. Specifically, Maxwell showed  that for elastic “Maxwell molecules” interacting via a repulsive $`r^{2(d1)}`$ potential in $`d`$ spatial dimensions, the collision rate does not depend on the magnitude of the relative velocity.

Let $`P(𝐯,t)`$ and $`Q(𝐯,t)`$ be the normalized velocity distributions of the background and the impurity, respectively. In this collision process, no correlations develop and the governing Boltzmann equation and the Lorentz-Boltzmann equation are

$`{\displaystyle \frac{P(𝐯,t)}{t}}`$ $`=`$ $`{\displaystyle 𝑑𝐧𝑑𝐮_1P(𝐮_1,t)𝑑𝐮_2P(𝐮_2,t)}`$ (5)
$`\times `$ $`\left\{\delta \left[𝐯𝐮_1+(1p)(𝐠𝐧)𝐧\right]\delta (𝐯𝐮_1)\right\},`$ (6)
$`{\displaystyle \frac{Q(𝐯,t)}{t}}`$ $`=`$ $`{\displaystyle 𝑑𝐧𝑑𝐮_1Q(𝐮_1,t)𝑑𝐮_2P(𝐮_2,t)}`$ (7)
$`\times `$ $`\left\{\delta \left[𝐯𝐮_1+(1q)(𝐠𝐧)𝐧\right]\delta (𝐯𝐮_1)\right\}.`$ (8)
The angular integration over the impact direction is normalized to unity, $`𝑑𝐧=1`$. The collision integrals directly reflect the collision rules.

In writing Eqs. (5)–(7), fluid-fluid and impurity-fluid collision rates were assumed to be the same and were set to unity for convenience. Hence the average number of collisions experienced by a particle equals time. The results detailed in this section are exact for a mean-field collision process where, irrespective of their type, random pairs of particles are chosen to undergo inelastic collisions according to (1)–(2). A more realistic collision rates model is treated in Sec. III.

The Boltzmann and Lorentz-Boltzmann equations simplify in Fourier space. Thence we use the Fourier transforms of the velocity distribution functions,

$`F(𝐤,t)`$ $`=`$ $`{\displaystyle 𝑑𝐯e^{i𝐤𝐯}P(𝐯,t)},`$ (9)
$`G(𝐤,t)`$ $`=`$ $`{\displaystyle 𝑑𝐯e^{i𝐤𝐯}Q(𝐯,t)}.`$ (10)
The convolution structure of Eqs. (5)–(7) implies that the collision terms factorize

$`{\displaystyle \frac{}{t}}F(𝐤,t)+F(𝐤,t)`$ $`=`$ $`{\displaystyle 𝑑𝐧F[𝐤𝐩,t]F[𝐩,t]},`$ (11)
$`{\displaystyle \frac{}{t}}G(𝐤,t)+G(𝐤,t)`$ $`=`$ $`{\displaystyle 𝑑𝐧G[𝐤𝐪,t]F[𝐪,t]},`$ (12)
with $`𝐩=(1p)(𝐤𝐧)𝐧`$ and $`𝐪=(1q)(𝐤𝐧)𝐧`$. These equations reflect the momentum transfer occurring during collisions. They also hint that the model is analytically tractable. For example, expansion in powers of the wave number shows that moments of the velocity distributions obey closed hierarchies of evolution equations.

### A The Temperature

We study the freely evolving case where in the absence of energy input the system “cools” indefinitely . The fluid temperature is defined as the average square of the velocity: $`T=\frac{1}{d}𝑑𝐯v^2P(𝐯,t)`$ with $`v|𝐯|`$. From the Boltzmann equation (5), the fluid temperature evolves according to $`\frac{d}{dt}T=\lambda T`$ with $`\lambda =2p(1p)𝑑𝐧n_1^2`$. The angular integration is readily performed, $`𝑑𝐧n_1^2=1/d`$, by using symmetry and the identity $`n_1^2+\mathrm{}+n_d^2=1`$. Hence the fluid temperature satisfies

$$\frac{d}{dt}T=\frac{2p(1p)}{d}T,$$
(13)
implying that the fluid temperature decays exponentially with time

$$T(t)=T_0e^{2p(1p)t/d}.$$
(14)
Similarly, from the Lorentz-Boltzmann equation (7) one can obtain the governing equation for the impurity temperature $`\mathrm{\Theta }(t)=\frac{1}{d}𝑑𝐯v^2Q(𝐯,t)`$. One finds that the impurity temperature is coupled to the fluid temperature via a simple linear rate equation

$$\frac{d}{dt}\mathrm{\Theta }=\frac{1q^2}{d}\mathrm{\Theta }+\frac{(1q)^2}{d}T.$$
(15)
The solution to this equation is a linear combination of two exponentials

$$\mathrm{\Theta }(t)=\left(\mathrm{\Theta }_0cT_0\right)e^{(1q^2)t/d}+cT_0e^{2p(1p)t/d},$$
(16)
with the constant

$$c=\frac{(1q)^2}{1q^22p(1p)}.$$
(17)
Therefore, there are two different regimes of behavior. When $`1q^2>2p(1p)`$, the impurity temperature is proportional to the fluid temperature asymptotically, $`\frac{\mathrm{\Theta }(t)}{T(t)}c`$ as $`t\mathrm{}`$. In the complementary region $`2p(1p)>1q^2`$ an extreme violation of equipartition occurs, as the ratio of the fluid temperature to the impurity temperature vanishes. The impurity is very energetic compared with the fluid and it practically sees a static fluid. From Eq. (17) we find that at the transition point $`q=\sqrt{12p(1p)}`$. Employing relations (3)–(4) between the restitution coefficients and the collision parameters we obtain the critical mass

$$m_{}=\frac{r_q+\sqrt{(1+r_p^2)/2}}{1\sqrt{(1+r_p^2)/2}}.$$
(18)
The heavy impurity phase arises when the impurity is a bit heavier than the fluid particles: Even when the fluid-fluid collisions are completely inelastic ($`r_p=0`$), the critical mass satisfies $`m_{}>1+\sqrt{2}`$. For weakly dissipative fluids ($`r1`$), the critical mass diverges, $`m_{}(1r_p)^1`$ (see Fig. 1).

Note now a few features of the light impurity phase. First, Eq. (17) generalizes the elastic fluid ($`p=0`$) result $`c=(1q)/(1+q)`$. That result was actually established for a hard sphere fluid, so at least asymptotically both the IMM and the improved model predict the same impurity temperature when the fluid is elastic. Further, the initial impurity temperature becomes irrelevant and the impurity is governed by the fluid background. The average energies of the impurity and fluid particles are asymptotically equal when $`m\mathrm{\Theta }/T1`$, or when $`mc=1`$. Thus, energy equipartition occurs on a particular surface in the three dimensional space $`(m,r_p,r_q)`$ where $`m(1q)^2=1q^22p(1p)`$. Using relations (3)–(4) between the collision parameters and the restitution coefficients, energy equipartition occurs when the impurity mass $`m_{\mathrm{eq}}`$ is given by

$$m_{\mathrm{eq}}=\frac{1+r_p^22r_q^2}{1r_q^2}.$$
(19)
As expected, this mass equals unity when $`r_p=r_q`$. Curiously, $`m_{\mathrm{eq}}`$ vanishes when $`r_q^2=(1+r_p^2)/2`$, indicating that for $`r_q>\sqrt{(1+r_p^2)/2}`$, energy equipartition does not occur, as shown in Fig. 2. However, in a generic point in the parameter space $`(m,r_p,r_q)`$ equipartition does break down . This is a signature of the dissipative and nonequilibrium nature of the system.

In summary, when the impurity mass is smaller than the critical mass $`m<m_{}`$, the fluid temperature governs the impurity temperature. Otherwise, the impurity is infinitely more energetic than the fluid asymptotically. As $`t\mathrm{}`$ one finds

$$\frac{\mathrm{\Theta }(t)}{T(t)}\{\begin{array}{cc}c\hfill & m<m_{}\text{,}\hfill \\ \mathrm{}\hfill & mm_{}\text{.}\hfill \end{array}$$
(20)
At the critical impurity mass, $`m=m_{}`$, the solution to Eq. (15) shows that the temperature ratio diverges linearly with time, $`\frac{\mathrm{\Theta }(t)}{T(t)}2p(1p)t`$. Interestingly, the dependence on the dimension is secondary as it only sets the overall time scale (the transformation $`tt/d`$ absorbs the dimension dependence). For example, both the critical mass, $`m_{}`$, and the temperature ratio $`c`$ are independent of $`d`$.

### B The light impurity phase

For a light impurity, $`m<m_{}`$, the fluid velocity statistics govern basic characteristics of the impurity statistics. We have seen this for the temperature, as the initial impurity temperature is irrelevant asymptotically, and the ratio of the two temperatures approaches a constant $`\frac{\mathrm{\Theta }(t)}{T(t)}c`$. The question arises how more detailed characteristics such as the impurity velocity distribution emerge from the fluid velocity distribution. This question can be answered in detail in one dimension where the explicit scaling solution for the fluid velocity distribution is known . Below we primarily focus on the one-dimensional case and derive an explicit exact solution. We then briefly comment on the higher-dimensional case, where the solution is very cumbersome.

The velocity distribution of the fluid approaches a scaling form asymptotically, $`P(v,t)T^{1/2}𝒫\left(vT^{1/2}\right)`$. Following the earlier treatment of the fluid velocity distribution, we employ the Fourier transform technique. The fluid Fourier transform is thus characterized by the scaling form $`F(k,t)=f\left(|k|T^{1/2}\right)`$. The governing equation for the corresponding scaling function is obtained from (11) to give

$$p(1p)zf^{}(z)+f(z)=f(pz)f(zpz),$$
(21)
and the solution is $`f(z)=(1+z)e^z`$ .

The behavior of the fluid suggests that the impurity velocity distribution also approaches a scaling form $`Q(v,t)T^{1/2}𝒬\left(vT^{1/2}\right)`$. Then the corresponding Fourier transform reads $`G(k,t)=g\left(|k|T^{1/2}\right)`$. From the Lorentz-Boltzmann equation (12), this scaling function satisfies the linear equation

$$p(1p)zg^{}(z)+g(z)=(1+az)e^{az}g(qz)$$
(22)
with $`a=1q`$. The fluid scaling function is a combination of $`z^ne^z`$ with $`n=0`$ and $`n=1`$. We therefore seek a similar solution of Eq. (22)

$$g(z)=\underset{n=0}{\overset{\mathrm{}}{}}A_nz^ne^z.$$
(23)
Inserting (23) into (22), the exponential factors cancel. Equating terms of the order $`z^n`$, yields the following recursion relation for the coefficients:

$$A_n=\frac{q^{n1}(1q)p(1p)}{1q^nnp(1p)}A_{n1}.$$
(24)
From the normalization $`g(0)=1`$ we get $`A_0=1`$. The next two coefficients are $`A_1=1`$, $`A_2=(1c)/2`$. Using these values, one computes the small $`z`$ expansion of the solution (23): $`g(z)=1cz^2/2+𝒪(z^3)`$, and the first non-trivial coefficient is indeed consistent the definition of the Fourier transform.

The Fourier transform can be inverted to obtain the impurity velocity distribution function explicitly. The inverse Fourier transform of $`e^{\kappa z}`$ is $`\frac{1}{\pi }\frac{\kappa }{\kappa ^2+w^2}`$; the inverse transforms of $`z^ne^z`$ can be obtained using successive differentiation with respect to $`\kappa `$. In general, the velocity scaling function is

$$𝒬(w)=\frac{1}{\pi }\underset{n=0}{\overset{\mathrm{}}{}}A_n(1)^n\frac{d^n}{d\kappa ^n}\frac{\kappa }{\kappa ^2+w^2}|_{\kappa =1}.$$
(25)
In other words, one can express the solution as a combination of powers of Lorentzians

$$𝒬(w)=\frac{2}{\pi }\underset{n=2}{\overset{\mathrm{}}{}}B_n\left(\frac{1}{1+w^2}\right)^n.$$
(26)
The coefficients $`B_n`$ are linear combinations of the coefficients $`A_k`$’s with $`kn+1`$, e.g., $`B_2=13A_2+3A_3`$ and $`B_3=4A_224A_3+60A_4`$. The first squared Lorentzian term dominates the tail of the velocity distribution $`𝒬(w)𝒫(w)w^4`$, as $`w\mathrm{}`$. This algebraic behavior prevails for $`w1`$. Therefore, the impurity has the same extremal velocity statistics as the fluid.

An interesting aspect of this solution is that for particular values of the collision parameter $`q`$, the infinite sum terminates at a finite order. Of course, when $`q_1=p`$ then

$$𝒬(w)=𝒫(w)=\frac{2}{\pi }\frac{1}{(1+w^2)^2}.$$
(27)
This is reflected by the vanishing coefficient $`A_2=0`$. The multiplicative recursion relation (24) then shows that $`A_n=0`$ for all $`n2`$. Furthermore, for $`q_2^2(1q_2)=p(1p)`$, one has $`A_3=0`$, and the velocity distribution (26) contains only two terms

$$𝒬(w)=\frac{2}{\pi }\left[\frac{13A_2}{\left(1+w^2\right)^2}+\frac{4A_2}{\left(1+w^2\right)^3}\right],$$
(28)
with $`A_2=q/(1+2q)`$. Similarly, the solution can be a finite sum with $`k`$ terms when $`p`$ and $`q_k`$ are related via $`q_k^k(1q_k)=p(1p)`$ thereby imposing $`A_n=0`$ for $`n>k`$. For given fluid collision parameter $`p`$, the above relation has the solution $`q_k`$ only for sufficiently small $`k`$. Thus, there are a few special values of the impurity collision parameter $`q_k`$ for which the scaled impurity velocity distribution is a linear combination of simple rational functions $`\left(1+w^2\right)^n`$ with $`n=2,\mathrm{},k+1`$.

While the scaling functions underlying the impurity and the fluid are similar, more subtle features may differ. In particular, the full time dependent behavior, as characterized by the moments of the impurity distribution, exhibits rich behavior. Let $`L_n(t)=𝑑vv^nP(v,t)`$ be the moments of the fluid velocity distribution. Multiplying the equations (5)–(7) by $`v^n`$ and integrating, the moments obey the recursive equations

$`{\displaystyle \frac{d}{dt}}L_n+a_nL_n`$ $`=`$ $`{\displaystyle \underset{j=2}{\overset{n2}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{j}}\right)p^j(1p)^{nj}L_jL_{nj},`$ (29)
$`{\displaystyle \frac{d}{dt}}M_n+b_nM_n`$ $`=`$ $`{\displaystyle \underset{j=0}{\overset{n2}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{j}}\right)q^j(1q)^{nj}M_jL_{nj},`$ (30)
with

$$a_n(p)=1p^n(1p)^n,b_n(q)=1q^n.$$
(31)
Asymptotically, the fluid moments decay exponentially according to

$$L_n(t)e^{a_n(p)t}.$$
(32)
Using this asymptotics we analyze the behavior of the impurity moments. The second moment, i.e. the impurity temperature, was already shown to behave similar to the fluid temperature when $`a_2(p)<b_2(q)`$. The fourth moment behaves similarly to the fourth moment of the fluid when $`a_4(p)<b_4(q)`$, and generally the first $`n`$ moments are proportional to each other, $`M_2L_2,\mathrm{},M_{2n}L_{2n}`$, when $`a_{2k}(p)<b_{2k}(q)`$ for $`k=1,\mathrm{},n`$. In particular, all respective impurity and fluid moments are proportional to each other when $`a_{2k}(p)<b_{2k}(q)`$ is valid for every $`k`$. It is easy to see that above inequalities hold when $`p+q<1`$, which is always obeyed when $`m<1`$.

However, in the complementary parameter range, $`a_n(p)>b_n(q)`$ for sufficiently large $`n`$. Such moments are no longer governed by the fluid and

$$M_n(t)e^{(1q^n)t}.$$
(33)
The corresponding moment ratio diverges asymptotically: $`M_n/L_n\mathrm{}`$ as $`t\mathrm{}`$. Interestingly, the same behavior (33) is found in the heavy impurity phase, as will be shown below. Therefore, the two phases are not entirely distinct. A series of transitions affecting moments of decreasing order occurring at increasing masses,

$$m_1>m_2>\mathrm{}>m_{\mathrm{}},$$
(34)
signals the transition to the heavy impurity phase. When $`mm_n`$, the ratio $`M_{2k}/L_{2k}`$ diverges asymptotically for all $`kn`$. This generalizes the second moment transition occurring at $`m_1m_{}`$. The transition masses

$$m_n=\frac{r_q+\frac{1}{2}\left[(1r_p)^{2n}+(1+r_p)^{2n}\right]^{\frac{1}{2n}}}{1\frac{1}{2}\left[(1r_p)^{2n}+(1+r_p)^{2n}\right]^{\frac{1}{2n}}}$$
(35)
are found from $`q^{2n}=p^{2n}+(1p)^{2n}`$ and Eq. (4). In particular, $`m_{\mathrm{}}=lim_n\mathrm{}m_n=(1+r_p+2r_q)/(1r_p)`$, and sufficiently light impurities ($`m<m_{\mathrm{}}`$) mimic the fluid completely. Additionally, since $`m_n>m_{\mathrm{}}1`$, the impurity must be heavier than the fluid for any transition to occur.

The above analysis of the light impurity phase suggests that in higher dimensions, the impurity velocity distribution might be similar to that of the fluid. In higher dimensions we again assume that the Lorentz-Boltzmann equation admits a scaling solution. The corresponding equation in Fourier space (12) then considerably simplifies. Following the earlier treatments of the fluid case we extract the high-energy tail from the small-$`k`$ behavior of the Fourier transform of the scaled velocity distribution. The outcome is that both velocity distributions have the same algebraic high-energy tail

$$𝒬(w)𝒫(w)w^\sigma ,$$
(36)
with the exponent $`\sigma `$ calculated in . Such analysis also yields the ratio of the prefactors governing this algebraic decay.

### C The heavy impurity phase

When the mass of the impurity is equal to or larger than the critical mass, $`mm_{}`$, the velocities of the fluid particles are asymptotically negligible compared with the velocity of the impurity. Hence, in the $`t\mathrm{}`$ limit, fluid particles become stationary as viewed by the impurity. Therefore, one can set $`𝐮_20`$ in the collision rule (2):

$$𝐯=𝐮(1q)(𝐮𝐧)𝐧.$$
(37)
This process is somewhat analogous to a Lorentz gas. However, in the granular impurity system, a heavy particle scatters of a static background of lighter particles, while in the Lorentz gas the scatterers are infinitely massive. Despite this difference, the mathematical descriptions of the two problems are similar. Specifically, the collision rule for the inelastic Lorentz gas  is obtained from (37) by a mere replacement of the factor $`(1q)`$ with $`(1+r_q)`$.

Let us first consider the one-dimensional case where an explicit solution of the velocity distribution is possible. Setting $`u_20`$ in the delta function in the Lorentz-Boltzmann equation (7), integration over the fluid velocity $`u_2`$ is trivial, $`𝑑u_2P(u_2,t)=1`$, and integration over the impurity velocity $`u_1`$ gives

$$\frac{}{t}Q(v,t)+Q(v,t)=\frac{1}{q}Q(\frac{v}{q},t).$$
(38)
This equation can be solved directly by considering the stochastic process the impurity particle experiences. In a sequence of collisions, the impurity velocity changes according to $`v_0qv_0q^2v_0\mathrm{}`$ with $`v_0`$ the initial velocity. After $`n`$ collisions the impurity velocity decreases exponentially, $`v_n=q^nv_0`$. The collision rate is unity, and hence, the average number of collisions experienced till time $`t`$ equals $`t`$. Furthermore, the collision process is random, and therefore, the probability that the impurity undergoes exactly $`n`$ collisions up to time $`t`$ is Poissonian $`t^ne^t/n!`$. Thus, the velocity distribution function reads

$$Q(v,t)=e^t\underset{n=0}{\overset{\mathrm{}}{}}\frac{t^n}{n!}\frac{1}{q^n}Q_0\left(\frac{v}{q^n}\right),$$
(39)
where $`Q_0(v)`$ is the initial velocity distribution of the impurity. Indeed, one can check that this properly normalized solution satisfies Eq. (38).

Interestingly, the impurity velocity distribution function is a time-dependent combination of “replicas” of the initial velocity distribution. Since the corresponding argument is stretched, compact velocity distributions display an infinite set of singularities, a generic feature of the Maxwell model.

The impurity velocity distribution exhibits interesting asymptotic behaviors. Consider for simplicity the uniform initial velocity distribution: $`Q_0(v)=1`$ for $`|v|<1/2`$ and $`Q_0(v)=0`$ otherwise. The solution (39) reduces to a finite sum, with $`nN=\frac{\mathrm{ln}(2v)}{\mathrm{ln}q}`$. In the physically interesting limits $`t\mathrm{}`$ and $`v0`$, the sum on the right-hand side of Eq. (39) simplifies to respectively $`e^{t/q}`$ or $`(t/q)^N/N!`$ when the number of terms is above or below the threshold value $`N=t/q`$. The magnitude of $`Q(v,t)`$ above the threshold greatly exceeds the magnitude below the threshold, so $`Q(v,t)`$ appears to approach a step function. A refined analysis shows that the width of the front widens diffusively so the front remains smooth although its relative width vanishes. Specifically, one finds the following traveling-wave like scaling solution

$$Q(v,t)e^{t+t/q}\mathrm{\Psi }(\eta ),$$
(40)
with the following wave form and coordinate

$$\mathrm{\Psi }(\eta )=\frac{1}{\sqrt{\pi }}_{\mathrm{}}^\eta 𝑑xe^{x^2},\eta =\frac{\frac{q}{\mathrm{ln}q}\mathrm{ln}(2v)t}{\sqrt{2qt}}.$$
(41)
Note, however, that the large velocity tail ($`\eta \mathrm{}`$), ignored in Eq. (41), provides actually the dominant contribution to the moments. This is an unusual traveling wave form in the sense that the argument is the logarithm of the velocity rather then the velocity itself, a reflection of the exponentially decaying velocity.

In contrast to the velocity distribution, the moments $`M_n(t)=𝑑vv^nQ(v,t)`$ exhibit a much simpler behavior. Indeed, from Eq. (38) one finds that every moment is coupled only to itself, $`\frac{d}{dt}M_n=(1q^n)M_n`$. Solving this equation we recover Eq. (33); in the heavy impurity phase, however, it holds for all $`n`$. Therefore the moments exhibit multiscaling asymptotic behavior. The decay coefficients, characterizing the $`n`$-th moment, depend on $`n`$ in a nonlinear fashion. This multiscaling behavior excludes scaling solutions with sharp tails (stretched exponentials decays and faster). While collisions with the fluid are sub-dominant, they still lead to corrections to the leading asymptotic behavior.

We turn now to arbitrary spatial dimensions $`d`$. The impurity velocity changes according to Eq. (37) with the impact direction $`𝐧`$ chosen randomly. The corresponding Lorentz-Boltzmann equation reads

$`{\displaystyle \frac{}{t}}Q(𝐯,t)+Q(𝐯,t)`$ $`=`$ $`{\displaystyle 𝑑𝐧𝑑𝐮Q(𝐮,t)}`$ (43)
$`\times \delta \left[𝐯𝐮+(1q)(𝐮𝐧)𝐧\right].`$
Moments of the velocity distribution

$$M_n(t)=𝑑𝐯v^nQ(𝐯,t),$$
(44)
can be obtained directly. We focus on the even moments of the distribution. Indeed, they satisfy the following evolution equation

$$\frac{d}{dt}M_{2n}=\left(1\xi ^n\right)M_{2n}$$
(45)
where $`\mu =\mathrm{cos}^2\theta =(\widehat{𝐮}𝐧)^2`$, $`\xi \xi (q,\mu )=1(1q^2)\mu `$, and $``$ is the shorthand notation for the angular integration: $`f_0^1𝒟\mu f(\mu )`$. Since $`d𝐧\mathrm{sin}^{n2}\theta d\theta `$, the (normalized) integration measure $`𝒟\mu `$ is

$$B(\frac{1}{2},\frac{d1}{2})𝒟\mu =\mu ^{\frac{1}{2}}(1\mu )^{\frac{d3}{2}}d\mu $$
(46)
where $`B(a,b)`$ is the beta function. For example, $`1=1`$, and $`\mu =1/d`$.

From the evolution equations (45), the moments are found to decay exponentially with time

$$M_{2n}(t)=M_{2n}(0)e^{\left(1\xi ^n\right)t}.$$
(47)
In particular, $`1\xi =(1q^2)/d`$, and thus, the temperature decay $`\mathrm{\Theta }(t)=\mathrm{\Theta }(0)e^{(1q^2)t/d}`$ is recovered. There is a little discrepancy with the exact temperature of Eq. (16) which depends on the initial fluid temperature, $`T_0`$. This is a remnant of a transient regime where the two velocity scales are comparable. Generally, the time dependence is correct. Note also that in one dimension $`\xi =q^2`$ and hence Eq. (47) reduces to our earlier result. In the infinite dimension limit, $`\mu 0`$, and all moments decay according to $`e^t`$. However, in general, the moments exhibit multiscaling asymptotic behavior, and knowledge of the typical velocity is insufficient to fully characterize the entire velocity distribution. Indeed, writing $`M_{2n}M_2^{\alpha _n}`$, the exponents $`\alpha _n=(1\xi ^n)/(1\xi )`$ have a nontrivial spectrum.

The moments directly give a formal exact solution of the Fourier transform of the impurity velocity distribution (10). We consider isotropic situations where $`G(𝐤,t)G(k^2,t)`$ with $`k|𝐤|`$. Expanding the transform in powers of $`k^2`$ and substituting the moment result (47) yields

$$G(𝐤,t)=e^t\underset{n=0}{\overset{\mathrm{}}{}}\frac{\left(k^2\right)^n\mu ^n}{(2n)!}M_{2n}(0)e^{\xi ^nt}.$$
(48)
While this is an explicit solution, it is not too illuminating. First, it is in Fourier space, and second, it involves the complicated angular averages $`\xi ^n`$.

Nevertheless, it can be shown that the solution remains a time dependent combination of properly modified replicas of the initial distribution. Indeed, either from Eq. (48) or directly from the Fourier transform equation $`\frac{}{t}G(𝐤,t)+G(𝐤,t)=G(𝐤\xi ,t)`$, the solution can be rewritten in the form

$$G(k^2,t)=e^t\underset{n=0}{\overset{\mathrm{}}{}}\frac{t^n}{n!}G_n(k^2).$$
(49)
Here, $`G_0(k^2)`$ is the initial Fourier transform, and the “building blocks” $`G_n`$ are obtained from a recursive procedure of angular integration $`G_{n+1}(k^2)=G_n(k^2\xi )`$.

## III The Improved Inelastic Maxwell Model

In the IMM, the rates for fluid-fluid and fluid-impurity collisions were identical (and therefore, set to unity for convenience). For granular fluids, however, the collision rate is proportional to the relative velocity . Therefore, one can improve the Maxwell model by replacing the actual collision rate with an average collision rate proportional to the average relative velocity. The simplest choice of the average relative velocity is $`\sqrt{(𝐯_1𝐯_2)^2}\sqrt{(T_1+T_2)/2}`$. For fluid-fluid and impurity-fluid collisions we thus obtain $`\sqrt{T}`$ and $`\sqrt{(T+\mathrm{\Theta })/2}`$, respectively. Thus, different collision rates multiply the collision integrals in the Boltzmann and Lorentz-Boltzmann equations (5)–(7). In the fluid case, this overall prefactor merely affects the time dependence of the temperature. As will be shown below, in the impurity case, the above phase transition is suppressed in the IIMM, although secondary transitions corresponding to higher order moments remain.

Let us again start with the behavior of the temperature. The fluid temperature satisfies

$$\frac{d}{dt}T=\sqrt{T}\left[\frac{2p(1p)}{d}T\right].$$
(50)
Solving this equation, we recover Haff’s cooling law $`T(t)=T_0[1+t/t_0]^2`$, with $`T_0`$ the initial temperature and $`t_0=d/[p(1p)T_0^{1/2}]`$ .

The corresponding rate equation for the impurity temperature $`\mathrm{\Theta }`$ is

$$\frac{d}{dt}\mathrm{\Theta }=\sqrt{\frac{T+\mathrm{\Theta }}{2}}\left[\frac{1q^2}{d}\mathrm{\Theta }+\frac{(1q)^2}{d}T\right].$$
(51)
Since we are primarily interested in the temperature ratio, $`S=\mathrm{\Theta }/T`$, we study this quantity directly. It evolves according to

$`{\displaystyle \frac{1}{\sqrt{T}}}{\displaystyle \frac{d}{dt}}S`$ $`=`$ $`\sqrt{{\displaystyle \frac{1+S}{2}}}\left[{\displaystyle \frac{1q^2}{d}}S+{\displaystyle \frac{(1q)^2}{d}}\right]`$ (52)
$`+`$ $`{\displaystyle \frac{2p(1p)}{d}}S.`$ (53)
In the inelastic Maxwell model, the gain and the loss terms were comparable, both increasing linearly with $`S`$. Here, in contrast, the loss term, which grows as $`S^{3/2}`$, eventually overtakes the gain term that grows only linearly with $`S`$. Therefore, $`Sc`$ where $`c`$ is the root of the cubic equation

$$\sqrt{\frac{1+c}{2}}\left(c\frac{1q}{1+q}\right)=\frac{2p(1p)}{1q^2}c.$$
(54)
Consequently, there is only one phase, the light impurity phase. We note that the ratio $`c`$ is independent of the spatial dimension $`d`$. Intuitively, since the impurity collision rate relatively increases with the impurity temperature, the impurity energy dissipation rate increases, thereby limiting the (relative) growth of the impurity temperature.

Generally, there is no equipartition of energy except for a particular surface in the space $`(m,r_p,r_q)`$. Energy equipartition occurs when $`mc=1`$. Using the relations (3)–(4), we find the equipartition mass

$$m_{\mathrm{eq}}=\frac{1}{2}+\frac{1}{2}\sqrt{1+8\left(\frac{1r_q^2}{1r_p^2}\right)^2}.$$
(55)
Figure 3 plots $`m_{\mathrm{eq}}=m_{\mathrm{eq}}(r_p,r_q)`$.

Qualitatively, our findings in the light impurity phase of the IMM extend to the IIMM. For example, both velocity distributions follow scaling forms and the large-velocity tails of both distributions are the same. In the one dimensional case, explicit expressions for the impurity scaling function are possible, and as the treatment follows closely that outlined in the light impurity phase, we briefly outline the results. In 1D, the impurity velocity distribution approaches a scaling solution $`Q(v,t)T^{1/2}𝒬\left(vT^{1/2}\right)`$. The corresponding Fourier transform reads $`G(k,\tau )=g\left(|k|T^{1/2}\right)`$. The only difference with the above inelastic Maxwell model is that the collision terms are proportional to $`\beta ^1=\sqrt{(1+c)/2}`$. Consequently, Eq. (22) generalizes as follows

$$\beta p(1p)zg^{}(z)+g(z)=(1+az)e^{az}g(qz).$$
(56)
Seeking a series solution of the form (23), leads to the following recursion relations for the coefficients

$$A_n=\frac{(1q)q^{n1}\beta p(1p)}{1q^nn\beta p(1p)}A_{n1},$$
(57)
with $`A_0=A_1=1`$. Again, the velocity distribution is a combination of powers of Lorentzians as in Eq. (26):

$$𝒬(w)=\frac{2}{\pi }\underset{n=2}{\overset{\mathrm{}}{}}B_n\left(1+w^2\right)^n,$$
(58)
where $`B_n`$ are linear combinations of the coefficients $`A_n`$’s given by the same expressions as in Sec. II. In particular, the large-velocity tail is generic $`𝒬(w)w^4`$.

Given the algebraic form of the velocity distributions, we examine the asymptotic behavior of moments of the velocity distribution. Moments of the fluid and the impurity, $`L_n`$ and $`M_n`$, respectively, evolve according to a straightforward generalization of Eqs.(29)–(30),

$`{\displaystyle \frac{d}{d\tau _p}}L_n+a_nL_n`$ $`=`$ $`{\displaystyle \underset{j=2}{\overset{n2}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{j}}\right)p^j(1p)^{nj}L_jL_{nj},`$ (59)
$`{\displaystyle \frac{d}{d\tau _q}}M_n+b_nM_n`$ $`=`$ $`{\displaystyle \underset{j=0}{\overset{n2}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{j}}\right)q^j(1q)^{nj}M_jL_{nj}.`$ (60)
Here $`a_n(p)`$ and $`b_n(q)`$ are given by Eq. (31) and the collision counters,

$`\tau _p={\displaystyle _0^t}𝑑t^{}\sqrt{T},\tau _q={\displaystyle _0^t}𝑑t^{}\sqrt{(T+\mathrm{\Theta })/2},`$
play the role of time in Eq. (59) and (60), respectively. Since the two temperatures are asymptotically proportional to each other, $`\tau _q\tau _p/\beta `$. The fluid moments decay according to $`L_ne^{a_n(p)\tau _p}t^{2a_n/a_2}`$ . Inserting the asymptotics $`L_ne^{a_n(p)\beta \tau _q}`$ into Eq. (60) and performing the same analysis as in the IMM we find that when $`\beta a_n(p)<b_n(q)`$, the impurity moments are enslaved to the fluid moments, i.e., $`M_nL_n`$ asymptotically. Otherwise, sufficiently large impurity moments behave differently than the fluid moments, viz. $`L_ne^{b_n\tau _q}t^{2b_n/\beta a_n}`$. Although the primary transition affecting the second moment does not occur, secondary transitions affecting larger moments do occur at a series of masses, as in Eq. (34). In the IIMM $`m_1m_{}`$ diverges, but other masses remain finite. The transition masses $`m_n`$ are found by solving $`\beta a_{2n}(p)=b_{2n}(q)`$ simultaneously with Eq. (54), and then applying Eq. (4). For example, for completely inelastic collisions ($`r_p=r_q=0`$) one finds $`m_2=1.65`$. Qualitatively, these transitions imply that some velocity statistics of the impurity, specifically large moments, are no longer governed by the fluid.

## IV Summary

We have studied dynamics of impurities in granular fluids using the inelastic Maxwell model. In general, there is breakdown of energy equipartition as two different temperatures characterize the impurity and the fluid. We analyzed two different models. First, we considered identical fluid-fluid and fluid-impurity overall collision rates. In this case a phase transition, marked by a dissipation dependent critical impurity mass, occurs. Breakdown of equipartition is moderate in one phase and extreme in the other with the asymptotic temperature ratio diverging. Sufficiently light impurities are governed by the fluid, and their scaled velocity distribution has similar extremal statistics as has the fluid. In one dimension, the scaled impurity velocity distribution is given by a series containing powers of Lorentzians of all orders (the infinite sum truncates for a few special values of the impurity restitution coefficient). Interestingly there is a series of increasing transition masses that corresponds to divergence of the ratio of velocity moments of decreasing order. These masses are always larger than unity and the largest such mass corresponds to the temperature. When the impurity mass is larger than this critical mass, the impurity and the fluid effectively decouple. Next, we considered collision rates that account for the temperature difference between the fluid and the impurity. In this case, the primary transition corresponding to the temperature is suppressed but further transitions corresponding to divergence of higher moments may still occur.

Many more quantities should be tractable in the framework of the Maxwell model. For example, one can study velocity correlations as well as the velocity autocorrelation functions. Furthermore, corrections to the leading asymptotic behavior, especially in the heavy impurity phase, can be systematically evaluated. Additionally, this model can be applied to particles of different sizes by modifying the collision cross sections.

We reiterate that the inelastic Maxwell model exactly describes only a mean-field collision process where the collision partners are random and the impact directions are chosen according to a uniform distribution. On the other hand, it is an uncontrolled approximation when applied to real granular fluids. At best, this model may be used to approximate sufficiently low velocity moments. Still, these results show that the relative collision rates play an important role. Comparison with the hard sphere Boltzmann equation, Molecular Dynamics simulations, and experiment is needed to gauge the utility of either of the two approximate models we treated. In practice, the impurity problem involves an additional level of difficulty. Gathering meaningful statistics requires many replica systems with one impurity or alternatively, the vanishing volume fraction limit of a mixture.

This research was supported by DOE (W-7405-ENG-36), NSF(DMR9978902), and ARO (DAAD19-99-1-0173).


format: text

text:
# Large-scale QRPA calculation of 𝐸1-strength and its impact on the neutron capture cross section

## 1 Introduction

The radiative neutron captures by exotic nuclei are known to be of fundamental importance in the rapid neutron-capture process (or r-process) invoked to explain the origin of approximately one half of the nuclides heavier than iron observed in nature. The r-process is believed to take place in environments characterized by high neutron densities ($`N_n>\mathrm{\hspace{0.17em}\hspace{0.17em}10}^{20}\mathrm{cm}^3`$), so that successive neutron captures can proceed into neutron-rich regions well off the $`\beta `$-stability valley. If the temperatures or the neutron densities characterizing the r-process are low enough to break the $`(n,\gamma )(\gamma ,n)`$ equilibrium, the r-abundance distribution depends directly on the neutron capture rates by the so-produced exotic neutron-rich nuclei . The neutron capture rates are commonly evaluated within the framework of the statistical model of Hauser-Feshbach. This model makes the fundamental assumption that the capture process takes place with the intermediary formation of a compound nucleus in thermodynamic equilibrium. In this approach, the Maxwellian-averaged $`(n,\gamma )`$ rate at temperatures of relevance in r-process environments strongly depends on the electromagnetic interaction, i.e the photon de-excitation probability.

The photon transmission coefficient from a compound-nucleus excited state is dominated by the $`E1`$ transition which is classically estimated within the Lorentzian representation of the giant dipole resonance (GDR), at least for medium- and heavy-mass nuclei. Reaction theory relates the $`\gamma `$-transmission coefficient for excited states to the ground state photoabsorption assuming the giant resonance to be built on each excited state. Experimental photoabsorption data confirms the simple semi-classical prediction of a Lorentzian shape at energies around the resonance energy $`E_{GDR}`$. For this reason, for practical applications involving neutron-capture cross section calculations, the $`E1`$-strength function is exclusively estimated so far on the basis of a form or another of the Lorentzian model. To explain the measured strength functions and $`\gamma `$-ray intensities , it is however required to bring correction to the low-energy behavior of the strength. This is achieved by introducing an energy-dependent width  or an energy- and temperature-dependent width  guided by more detailed finite Fermi system descriptions  of the $`E1`$-strength function at energies below the neutron separation energy $`S_n`$. These corrections were shown to improve significantly the predictions of the experimental radiation widths and gamma-ray spectra .

In addition to the generalized Lorentzian model , improved description of the $`E1`$-strength function can be derived from the semi-classical thermodynamic pole approach , the microscopic random phase approximation (RPA) calculations (based on Hartree-Fock  or relativistic mean field  ground state descriptions) or the theory of finite Fermi systems . None of the so-called microscopic approaches provides at the moment a complete set of $`E1`$-strength function for practical applications, so that global predictions for unknown nuclei can only be taken from phenomenological approximation of Lorentzian type. However, extrapolating these approximations (including the location of the resonance energy) to exotic neutron-rich nuclei might not be reliable. In particular, Catara et al.  showed in the framework of RPA calculations that the neutron excess affects the spreading of the isovector dipole strength, as well as the centroid of the strength function. The energy shift is found larger than the one accounted by the usual $`A^{1/6}`$ or $`A^{1/3}`$ dependence given by the phenomenological liquid drop approximations. Various calculations  on the exotic neutron-rich nuclei predict that the existence of a neutron mantle could introduce a new type of collective mode corresponding to an out-of-phase vibration of the neutron-proton core against the neutron mantle. The restoring force for this soft dipole vibration is predicted to be smaller than that of the GDR. The total $`E1`$ strength in this so-called pygmy resonance (PR) is small (around a few percent of the total GDR strength), but, if located well below the neutron separation energy, can significantly increase the radiative neutron capture cross section . Van Isacker et al.  estimated the PR energy and strength relative to that of the GDR within the incompressible fluid model of Goldhaber-Teller, assuming the dynamics of the dipole oscillations to be governed by the neutron-proton interaction. Such calculations were then systematically applied to all the neutron-rich nuclei and shown to affect significantly (up to factors of about 100) the neutron capture rates of the $`2<S_n[\mathrm{MeV}]<\mathrm{\hspace{0.17em}\hspace{0.17em}4}`$ closed shell nuclei . Large uncertainties in the description of the PR (in particular its energy and strength) remain and only more fundamental descriptions could shed light on its existence, as well as its relative importance and impact on neutron capture rates. In addition, different measurements suggest that some particular enhancement of the $`E1`$-strength could be located at low energies even on stable nuclei, a feature that cannot be predicted by the phenomenological approach. In particular, dipole transitions to bound states investigated by means of the nuclear resonance fluorescence confirmed the systematic existence of a pygmy $`E1`$-resonance at energies below the neutron separation energy . PR are observed in $`fp`$-shell nuclei as well as in heavy spherical nuclei near closed shells (Zr, Mo, Ba, Ce, Sn and Pb).

For various type of applications, and more particularly for nucleosynthesis applications, but also accelerator driven systems, a great effort must be made to improve the reliability of the $`E1`$-strength predictions, in particular for experimentally unreachable nuclei and in the low-energy region (typically below the neutron separation energy). Generally speaking, the more microscopic the underlying theory, the greater will be one’s confidence in the extrapolations out towards the neutron-drip line, provided, of course, the predictions are capable of reproducing experimentally known data with a sufficient degree of accuracy. For this reason, use should be made preferentially of microscopic or semi-microscopic global predictions based on sound and reliable nuclear models which, in turn, can compete with more phenomenological highly-parameterized models in the reproduction of experimental data. Microscopic approaches (like the RPA) are almost never used for practical applications, because of the numerical difficulty associated with large-scale predictions and the fine tuning of the model to reproduce accurately experimental data, especially when considering a large data set. We present here the first attempt to construct a complete set of $`E1`$-strength function derived from Quasi-particle RPA (QRPA) calculations based on an effective nucleon interaction of Skyrme-type. In Sect. 2, the QRPA calculations are presented. Different Skyrme interactions are compared in their predictions of the $`E1`$-strength and GDR energies. In Sect. 3, the QRPA shortcoming in broadening the $`E1`$-strength is cured by folding the QRPA strength on a Lorentzian-type function. A folding prescription is provided to reproduce experimental photoabsorption and available experimental $`E1`$ strengths at low energies, in particular average resonance capture (ARC) data. Section 4 is devoted to the analysis of exotic neutron-rich nuclei and the impact of the newly-determined strength on the neutron-capture cross sections and rates of interest for astrophysical applications. Conclusions are given in Sect. 5

## 2 QRPA calculations of the $`E1`$-strength

Microscopic models based on the RPA are well known for their accuracy to describe collective modes such as giant resonances . Due to the mean field approximation, they are able to calculate the strength function on the whole nuclear mass range in a short amount of time. An improved RPA model is the QRPA which takes the pairing effect into account. This effect is expected to be important for open-shell nuclei. In our case the ground state is calculated with a Hartree-Fock+BCS (HF+BCS) model. It should be noted that it is the first time that large scale calculations on more than 6000 nuclei are performed in order to calculate not only their ground state (with HF+BCS) but also their excited states (with the QRPA). This allows to characterize the structure model on the whole nuclear chart, with the effective nucleon-nucleon interaction as only input, and opens the exciting perspective to study the properties of the effective interaction with respect to the nuclear excitations. However the BCS approximation is known to not properly describe the pairing correlation for many drip-line nuclei . Ideally, for very neutron-rich nuclei along the drip-line, the model should take into account the coupling between the nuclear mean-field, the pairing and the particle continuum. Different approaches, such as the coordinate space HFB theory  or the continuum QRPA equations derived from the Time-Dependent-HFB theory  have been developed. However, large-scale calculations with such sophisticate continuum-QRPA models are far from being feasible at the present time. For this reason, we will restrict ourselves to the QRPA based on HF+BCS to get a qualitative insight of this approach on the neutron capture cross section by exotic neutron-rich nuclei. This study also provides interesting information on how calculated cross-sections obtained with two different approaches compare, namely the phenomenological Lorentzian (exclusively used nowadays in cross section calculation) and the microscopically rooted QRPA.

The QRPA model here employed  has already been successful in describing low-lying collective states in the oxygen , sulfur and argon  isotopic chains. The main feature of this model is its self-consistency in the sense that the mean field and the residual interaction are calculated from the same Skyrme effective interaction. The ground state is calculated within the HF+BCS approximation. The HF and the BCS equations are both solved in a self-consistent procedure. The only inputs to the calculations are the Skyrme interaction and the value of the constant pairing gap $`\mathrm{\Delta }`$. For the sake of simplicity, we choose a gap given by :

$$\mathrm{\Delta }=12A^{\frac{1}{2}}\mathrm{MeV}.$$
(1)
Because of the constant gap approximation we use a cutoff beyond which subshells do not participate to the pairing effect. This cutoff was chosen to take into account all the subshells of the major shell where the Fermi level belongs to. The pairing gap (1) is also used for magic nuclei. In this case the effect is not important because of the large shell spacing, so that the giant dipole resonance strength is not affected. Note that the present calculation with constant pairing gap is self-consistent only in the particle-hole channel, but not in the particle-particle channel. On top of the HF+BCS ground state, the excited $`J^\pi `$=1- states are calculated within the QRPA model. The QRPA equations are solved in the configuration space so as to exhaust the energy weighted sum rule (EWSR). All QRPA calculations are performed in the spherical approximation. The deformation effects are introduced in a phenomenological way as explained in Sect. 3. The results on strongly deformed nuclei should therefore be taken with circumspection. In E1-strength calculations a well-known spurious state appears, which corresponds to a translation of the whole nucleus. In our case, since we do not use a projection method or a modified operator, the spurious state is unambiguously identified in each E1-strength calculation : it lies close to 0 MeV (in a range between 0 and 2 MeV) and takes almost all the isoscalar strength (around 90%). Once identified this state is eliminated from the E1-strength.

The most significant observable to check the validity of the model and the accuracy of the interaction is the GDR energy. Several Skyrme forces have been used, namely SIII , SGII , SLy4  and MSk7 . In all cases, the pairing interaction is treated within the simple above-described prescription (Eq. (1)). In order to select the most adequate interaction to the $`E1`$ calculations, we compare the calculated centroid and measured GDR energies for 48 spherical nuclei . The root mean square deviations between the calculated centroid and measured energies are, respectively, 2267 keV for SIII, 573 keV for SGII, 457 keV for SLy4 and 564 keV for MSk7. The results are accurate for SGII, SLy4 and MSk7 interactions, as shown in Fig. 1 for the specific case of SLy4. Except for the lightest ($`A<\mathrm{\hspace{0.17em}\hspace{0.17em}80}`$) nuclei, all GDR energies are predicted within a few hundred keVs. On this basis we decide to use SLy4 for the large-scale calculation of the $`E1`$-strength and the corresponding radiative neutron capture. However, it should be noted that the agreement obtained with MSk7 is also good, although for some given nuclei (in particular 130Te and 133Cs) the GDR energy is larger than the experimental value by more than 1 MeV. This force has an effective mass of M/M=1.05 and was fitted exclusively to the experimental masses of 1888 nuclei, with a root mean square error of 738 keV . This result shows that an effective interaction like MSk7 characterized by a relatively large isoscalar and isovector effective mass (both equal to 1.05 of the nucleon mass) can reproduce the nuclear mass as well as GDR data. A low isovector effective mass of M$`{}_{v}{}^{}{}_{}{}^{}`$/M $`0.70`$ might therefore not be a necessary condition to reproduce photoabsorption cross sections in contrast to .

## 3 Folding procedure of the $`E1`$-strength

As explained in Sect. 2, the QRPA provides an accurate description of the GDR centroid and the fraction of the EWSR exhausted by the $`E1`$ mode. However, it is necessary to go beyond this approximation scheme in order to describe the damping properties of the collective motion. The GDR is known experimentally to have a large energy width and therefore a finite lifetime. Different microscopic theories exist (see e.g ). However, we will restrict ourselves, for the sake of simplicity and applicability to large-scale calculations, to an empirical broadening of the GDR. Such a broadening is obtained by folding the QRPA strength by a normalized Lorentzian function. Different types of Lorentzian functions as well as prescription for the Lorentzian width are considered. The folding prescription significantly affects the extrapolation of the $`E1`$-strength at low energies. In these conditions, ARC data at low excitation energies are of particular interest to constrain the folding procedure. The $`E1`$-strength is finally obtained by broadening each QRPA resonance strength at an energy $`E_i`$ by the normalized Lorentzian function

$$f_L(E,E_i,\gamma _i)=\frac{2}{\pi }\frac{\gamma _iE^2}{(E^2E_i^2)^2+\gamma _i^2E^2}.$$
(2)
A width $`\gamma _i=\mathrm{\Gamma }_{GDR}/2`$ (where the GDR width $`\mathrm{\Gamma }_{GDR}`$ is taken from experimental data  when available or, otherwise, from empirical systematics ) is found to lead to results in good agreement with experimental GDR widths and ARC data. The resulting folded GDR peak energies obtained with the SLy4 Skyrme force are compared with experimental values in Fig. 1 for spherical and deformed nuclei (see below for the treatment of deformation). As shown in Fig. 2 in the specific case of the spherical 144Nd nucleus, the ARC data for the 143Nd(n,$`\gamma `$)144Nd reaction  are also correctly reproduced at low energies. In addition, we compare in Fig. 3 the QRPA predictions with the compilation of experimental $`E1`$ strength functions at low energies ranging from 4 to 8 MeV  for nuclei from $`{}_{}{}^{25}\mathrm{Mg}`$ up to $`{}_{}{}^{239}\mathrm{U}`$. The data set includes resolved-resonance measurements, thermal-captures measurements and photonuclear data. In a certain number of cases the original experimental values need to be corrected, typically for non-statistical effects, so that only values recommended by  are considered in Fig. 3. QRPA predictions are found to be globally in good agreement with experimental data at low energies in the whole nuclear chart and of the same degree of accuracy as the recent calculation within the thermodynamic pole approximation .

Note that the GDR width is expected to be temperature-dependent . A temperature-dependent width $`\gamma _i`$ affects the folded $`E1`$-strength function only at high temperatures ($`T>\mathrm{\hspace{0.17em}\hspace{0.17em}2}`$ MeV) or if some significant strength is concentrated at very low energies ($`E_i<\gamma _i`$), but does not influence the predictions of the low-energy data shown in Figs. 2-3. Note also that introducing a temperature-dependent spreading width $`\gamma _i`$ in Eq. 2 does not give rise to a zero $`E0`$ limit of the $`E1`$-strength function due to the specific characteristics of the Lorentzian function considered in the folding procedure.

In the case of deformed spheroidal nuclei, the GDR is known to split into two major resonances as a result of the different resonance conditions characterizing the oscillations of protons against neutrons along the axis of rotational symmetry and an arbitrary axis perpendicular to it. In the phenomenological approach, the Lorentzian-type formula is generalized to a sum of two Lorentzian-type functions of energies $`E_{GDR}^l`$ and width $`\mathrm{\Gamma }_{GDR}^l`$ , such that

$`E_{GDR}^1+2E_{GDR}^2=3E_{GDR}`$ (3)
$`E_{GDR}^2/E_{GDR}^1=0.911\eta +0.089`$
where $`\eta `$ is the ratio of the diameter along the axis of symmetry to the diameter along an axis perpendicular to it. In turn, when not available experimentally, the width $`\mathrm{\Gamma }_{GDR}^l`$ of each peak can be expressed as a function of the respective energy $`E_{GDR}^l`$ . A similar splitting of the resonance strength is used in the folding procedure given by Eq. 2, each resonance energy $`E_i`$ being divided with an equal strength according to Eq. 3 and characterized by a width $`\gamma _i`$ proportional to $`\mathrm{\Gamma }_{GDR}^l`$ following the above-mentioned relation. Note that distributing the strength equally between the two resonance peaks (and not twice more on the high energy peak as done in the phenomelogical approach) is found to give optimal location and relative strength of both centroid energies as observed experimentally. It is shown in Fig. 1 that the splitting between both energy peaks can be rather well reproduced by the present correction. The strength at low energies is supposed to be affected in a similar way, although there is no experimental indication as such. This prescription affects the low-lying strength in a way similar to the folding procedure, i.e to spread the strength around the corresponding energy.

We illustrate in Fig. 4 how the QRPA calculation reproduces the photoabsorption cross sections of four (spherical and deformed) nuclei in the whole GDR region . The deformation parameter $`\eta `$ is extracted from the HF+BCS mass calculation of . In particular, the GDR energy and width, as well as the double peak structure observed experimentally are rather well reproduced by the QRPA calculation including the above-described folding and deformation prescriptions. In the case of deformed nuclei, the relative strength between both splitted peaks is also in good agreement with photonuclear data, in particular for 181Ta. Note that the double peak structure in deformed nuclei, is not solely due to deformation effects, but also to a E1-strength unequally distributed between 10 and 15 MeV. Although the strongest spherical QRPA peak is located around E=13 MeV, extra strength is found in the whole $`10<E[\mathrm{MeV}]<\mathrm{\hspace{0.17em}\hspace{0.17em}15}`$ energy range and is responsible for the asymmetries of the peaks. Experimentally, the $`E1`$ distribution is generally (at least in the rare-earth region) characterized by a larger strength at high energies than at low energies (see e.g. 181Ta). However, an inversion of this GDR strength asymmetry with a low-lying peak stronger than the high-energy one can be found in other region and is even recommended by the Obninsk evaluated photonuclear library , in particular for 234U, 235U (see Fig. 4), 239Pu and 241Pu. This recommended inversion, though not clearly observed in the 1976 experimental data of  illustrated in Fig. 4, is predicted by our QRPA calculation with the above-mentionned correction for deformed nuclei. Note that distributing unequally the strength between the two splitted peaks (for example with a strength twice larger in the high-energy state) would not reproduce as well experimental data as shown in Fig. 4.

Finally, little experimental systematics exists about the exact location and strength of the PR. Information obtained by the nuclear resonance fluorescence technique  suggests that the PR energy decreases with increasing masses, being around 8 MeV in the $`fp`$-shell region, 6.5 MeV around $`N=50,Z=50`$ and 5.5 MeV for Pb. In contrast, the PR strength increases with increasing masses, namely $`B(E1)=2060`$ (in $`10^3`$ e2fm2) for $`fp`$-shell nuclei, $`7090`$ in the $`N=50`$ region, $`80180`$ in the $`Z=50`$ region, and reaching 250 for 208Pb. Such a systematic behavior is also found by our QRPA calculation, but the low-lying strength is located systematically some 3 MeV higher than observed. Improvements of the QRPA model based on the quasiparticle-phonon coupling could shift the PR peak energy to lower values. In particular, Colò et al. (2001)  found that the observed low-lying dipole strength in the O isotopes can be explained by the QRPA plus phonon coupling model. Large-scale calculations in this framework are extremely computer-time consuming and remain to be performed.

## 4 Extrapolation to neutron-rich nuclei and application to the radiative neutron capture

Large-scale QRPA calculations based on the SLy4 Skyrme force (See Sect. 2) have been performed for all $`8Z110`$ nuclei lying between the proton and the neutron driplines, i.e some 6000 nuclei. All QRPA $`E1`$-distributions are folded by a Lorentzian-type function as explained in Sect. 3, although the validity of the phenomenological folding procedure adopted in the present work cannot be confirmed at this stage. Only further improved calculations including a detailed description of the deformation and damping effects can shed light on this assumption. In the neutron-deficient region, as well as along the valley of $`\beta `$-stability, the resulting QRPA strength functions are very similar to the empirical Lorentzian-like approximation and therefore in good agreement with experimental photoabsorption data (Fig. 4). When dealing with neutron-rich nuclei, the QRPA predictions start deviating from a simple Lorentzian shape. In particular, some extra strength is found to be located at an energy lower than the GDR energy. The more exotic the nucleus, the stronger this low-energy component. This is illustrated in Fig. 5 for the $`E1`$-strength function in the Sn isotopic chain. We should stress here that several sound nuclear models predict different behavior for the two-neutron separation energies of the Sn isotopes with A $`>`$ 132 . We present here the results for our QRPA model in order to illustrate qualitatively the role of a microscopic description of the E1-strength for neutron-rich nuclei. Among the 8 distributions shown in Fig. 5, only the $`A=150`$ one corresponds to a deformed configuration responsible for the double peak structure. For the $`A140`$ neutron-rich isotopes, an important part of the strength is concentrated at low energies ($`E<\mathrm{\hspace{0.17em}\hspace{0.17em}5}7`$ MeV). However, it should be stressed that for nuclei close to the drip-line, some spurious effect may be expected due to the lack of correct continuum treatment by our model. Microscopic models such as RPA and QRPA are well suited to interpret excitations in terms of particle-hole (p-h) configurations, each configuration having a specific weight in the excited state while the pairing effects allow additional configurations with holes above the Fermi level and particles below. In this framework the low energy component arises from two distinct origins, mainly depending on the mass of the nucleus. For rather light or medium nuclei, one or two particular (p-h) configurations are found responsible for the spectacular increase of the low energy part of the $`E1`$-strength in very neutron-rich nuclei such as 82Ni. For heavier elements, this increase is due to more collective modes. For instance in the 258Pb case, around 10 (p-h) configurations are participating to this low energy component with weights varying from 5% to 37%.

In 150Sn, mainly one (p-h) configuration is participating to the low energy component with a 87% weight. Phenomenological models are unable to predict such low energy components, whatever their collectivity is. In particular for 150Sn, all phenomenological systematics (as used for cross section calculation) predict a $`\gamma `$-ray strength peaked around 15 MeV with a full width at half maximum of about 4.5 MeV  which is obviously very different from the microscopic estimate (Fig. 5). The above-described feature of the QRPA $`E1`$-strength function for nuclei with a large neutron excess is found to be qualitatively independent of the adopted effective interaction. In particular, a similar enhanced strength is predicted by the MSk7 force which is characterized by a relatively large nucleon effective mass. Quantitatively, the total strength at low energies with the MSk7 interaction is smaller than in the SLy4 case.

The radiative neutron capture cross section is estimated within the statistical model of Hauser-Feshbach making use of the MOST code . This version benefits in particular from an improved description of the nuclear ground state properties derived from the microscopic Hartree-Fock method , as well as from a reliable nuclear level density prescription based on the microscopic statistical model . The direct capture contribution as well as the possible overestimate of the statistical predictions for resonance-deficient nuclei are effects that could have an important impact on the radiative neutron captures by exotic nuclei , but will not be included in the present analysis. We show in Fig. 6 for the Sn isotopes, the ratio of the radiative neutron capture cross section (at an arbitrary incident energy of 10 MeV exemplifying a high energy regime for which the statistical approximation of the Hauser-Feshbach model remains valid) predicted with the $`E1`$-strength taken from the QRPA to the one taken from the Hybrid Lorentzian-type formula . The temperature-dependent Hybrid formula corresponds to a generalization of the energy- and temperature-dependent Lorentzian formula including an improved description of the $`E1`$-strength function at energies below $`S_n`$ as derived from . The strength of the hybrid formula differs from the QRPA estimate not only in the location of the centroid energy, but also in the low-energy tail determined by the temperature and energy-dependent width. No extra low-lying strength is included in the phenomenological Hybrid formula, but its temperature dependence increases the E1 strength at low energies and is responsible for its non-zero $`E0`$ limit. The newly-derived $`E1`$-strength gives an increase of the cross section by a factor up to 20 close to the neutron drip line. This is due to the shift of the GDR to lower energies compared with the usually adopted $`A^{1/3}`$ rule (see Sect. 1 and Fig. 5), as well as to the appearance of the GDR strength at low energies as explained above. Both effects tend to increase the $`E1`$ strength at energies below the GDR, i.e in the energy window of relevance in the neutron capture process. For less exotic nuclei, the QRPA impact is relatively small, differences being mainly due to the exact position of the GDR energy and the resulting low-energy tail, i.e the energy-dependent width.

The Maxwellian-averaged $`(n,\gamma )`$ rates at a temperature of $`1.510^9`$ K characteristic of the r-process nucleosynthesis are also displayed in Fig. 6. These rates are sensitive to the neutron capture cross section at incident energies around 130 keV, and therefore depend on the $`E1`$-strength in an narrow range of a few hundred keV around $`S_n`$. For nuclei with low neutron separation energies, the resonant contribution around $`S_n`$ can be negligible, so that the stellar rate is principally dominated by captures of a few MeV neutrons. Deviations up to a factor of 5 are found with respect to the Hybrid formula (Fig. 6). The GDR strength predicted at low energies (see Fig. 5) is mainly responsible for the QRPA increase of the stellar reaction rate. In particular, below the neutron shell closure at $`N=82`$, the neutron separation energy is large enough (5–7 MeV), so that the GDR tail resulting from the low-energy QRPA strength enhances the reaction rate by a factor of 2. Above $`N=82`$, the neutron separation energy drops to low values ($`S_n<\mathrm{\hspace{0.17em}\hspace{0.17em}3}`$ MeV) and the reaction rates increases proportionnally to the $`E1`$ strength shifted to lower energies for the most exotic nuclei. For example in the $`{}_{}{}^{144}\mathrm{Sn}(n,\gamma )^{145}\mathrm{Sn}`$ case, the stellar rate and the cross section at 10 MeV are reduced by a factor of 2 and 3, respectively, if the $`E1`$ strength below 7 MeV is neglected (Fig. 5). These results show qualitatively the important effect the extreme neutron excess can have on the neutron capture cross section.

Similarly, the QRPA calculations have been used to estimate the neutron capture cross sections and Maxwellian-averaged rates for the 4000 nuclei with $`8Z84`$ lying between the neutron and the proton drip lines. They are compared in Figs. 7-8 to the predictions based on the Hybrid formula. Increases up to a factor of 30 are obtained in the cross section at 10 MeV and a factor of 8 in the rate at $`1.510^9`$ K, in both cases for the most exotic neutron-rich nuclei. A small decrease of the neutron capture cross sections and rates can also be observed, especially in the deformed neutron-deficient region ($`S_n10`$ MeV, $`Z70`$ and $`N105`$) where the neutron separation energy is large and the QRPA predicts, in the respective energy range of relevance, an $`E1`$-strength smaller than does the Hybrid formula. In this case, the neutron capture cross section with the Hybrid formula can be up to twice larger than the one obtained with the QRPA strength. This effect is closely related to the folding procedure adopted here, leading to a zero limit of the strength function at zero energy in contrast to the temperature-dependence of the Hybrid prescription leading to a non-zero limit . Qualitatively, similar results are obtained with the MSk7 predictions of the $`E1`$-strength. Quantitatively, the increase of the neutron capture cross section by nuclei close to the neutron drip line is slightly smaller with MSk7 than with SLy4. These qualitative results points out the necessity to perform more predictive calculations, such as continuum-QRPA, to have a reliable description of nuclei along the drip-line. The impact of the newly-derived rates on the r-process nucleosynthesis will be studied in a forthcoming paper.

## 5 Conclusions

The total photon transmission coefficient from a compound nucleus excited state is one of the key ingredients for statistical cross section evaluation. The cross section for the radiative capture of low-energy neutrons strongly depends on the low-energy tail of the giant dipole resonance. Recent experimental and theoretical work predict a particular enhancement of the $`E1`$-strength at low energies, in particular for neutron-rich nuclei. This low-energy component of the $`E1`$-strength can only be derived from microscopic calculations. Microscopic descriptions of the $`E1`$-strength function have become available, mainly in the QRPA approach. None of the previous QRPA calculations, however, provides a complete set of $`E1`$-strength functions to be used for practical applications.

For these reasons, we present here the first attempt to construct a complete set of $`E1`$-strength function derived from QRPA calculations based on an effective nucleon-nucleon interaction of Skyrme-type. The QRPA predictions using the SLy4 Skyrme force are in close agreement with photoabsorption data. A folding procedure is applied to the QRPA strength to take the damping of the collective motion into account. A prescription for the folding procedure and for the deformation effects is proposed and shown to lead to $`E1`$-strength functions in fair agreement with the available experimental data, in particular at low energies.

The folded QRPA strength is used for the evaluation of the neutron capture cross section within the framework of the statistical model of Hauser-Feshbach. The cross sections at energies around 10 MeV are found to be enhanced by a factor up to 30 when approaching the neutron drip line. The low-energy contribution to the $`E1`$ strength predicted by the QRPA calculations is held responsible for such an increase. No empirical prescription of the $`E1`$-strength as traditionally used in reaction calculation can predict this increase of the radiative capture. This emphasizes the crucial role of nuclei along the drip-lines, and more predictive calculations in this area are called for. The neutron capture rate at temperatures of relevance to nucleosynthesis applications is however less sensitive to the GDR strength calculated here. Increase up to a factor of 8 are obtained. The QRPA strength is tabulated and made available to the scientific community via the nuclear astrophysics library at http://www-astro.ulb.ac.be.

Further improvements might be necessary for a better description of the low-energy tail of the GDR. In particular, the particle-vibration coupling is also known to affect the low-energy strength and could contribute to an extra increase of the radiative neutron capture rate by exotic nuclei. A more fundamental determination of the GDR width, i.e of the damping of the collective motion, as well as a more consistent treatment of deformation effects also remain to be studied in the future. A proper treatment of the pairing and the coupling to the continuum for nuclei along the drip-line is also crucial to have a quantitative idea of the increase of the neutron capture cross-sections.

Acknowledgments S.G. is FNRS Research Associate. This work has been performed within the scientific collaboration (Tournesol) between the Wallonie–Bruxelles Community and France.