text:
# Quasi-Mössbauer effect in two dimensions

## I Introduction

In nuclear reactions, the atoms usually are treated as free, i.e., unbound. The chemical bonds between each other are neglected. This is justified with the argument that the typical energy in a nuclear process is generally much higher than the chemical binding energies ($`1`$ to $`10`$ eV) and moreso than that of lattice vibrations ($`10^2`$ to $`10^1`$ eV) wertheim . In return, most of the detailed qualities of the nucleus are of secondary importance in solid state physics.

These assumptions are not correct without limitations. Consider for example the absorption or emission of $`\gamma `$-photons by nuclei. During absorption processes (emission analogously) the nucleus first is in a ground state and makes a transistion into an excited state. The energy difference between the ground and the excited state is not constant. It has a certain distribution which is linked to the life-time of the excited state by Heisenberg’s uncertainty relation. For these reactions the difference has values between $`10^{11}`$ to $`10^6`$ eV wertheim . Thus, it is considerably smaller than the energy of chemical bonds.

During the absortion event, the atom would recoil due to the conservation of the sums of the momenta of the atom before the reaction and the absorbed photon. The thermal movement of an atom in a gaseous absorber leads to a Doppler broadening of the absorption line and the recoil to a shift of its maximum. Typical recoil-energies lie in the range from $`10^5`$ to $`10^{+2}`$ eV wertheim . This interval also comprises the energy of chemical bonds, phonons, and even energetically lower areas.

For atoms bound in a system and photon momenta below the threshold of the destruction of the chemical bonds, there are different aspects coming into play. For molecular systems additional rotational and vibrational degrees of freedom can be important. Simultaneously, the Doppler broadening and the recoil of the system decrease due to the enhanced mass of the system.

Making the transition to a crystal, the Doppler broadening and the recoil of the entire system become negligible. The entire recoil-energy goes into lattice vibrations. Further, there is the possibility that the momentum of the photon is taken up by the entire crystal without producing any phonons. Such an event is called recoil-free. In a rigid crystal, the absorption spectrum would just show the pure absorption spectrum of the nucleus, because no phonons can be produced. For photons with very low energy, i.e., less than $`1.5\times 10^{+5}`$ eV the recoil-energy is below the energy of the phonons. Classically, no vibrational modes of the crystal could be excited, thus the energy and momentum conservation could not be satisfied and the related processes could not take place. However, observations show that they do take place and that in those situations a part of the events are recoil-free. Energy and momentum conservation are guaranteed when averaging over many events lipkin . The effect and the fraction of recoil-free events have been named after its discoverer R. L. Mössbauer moessbauer .

The above discussion also plays a rôle in the scattering of x-rays or neutrons from atoms. There, the distinction between elastic and inelastic events is also known. Mössbauer used also results by Lamb about neutron scattering lamb . This is the reason for the name Mössbauer-Lamb fraction.

Of interest is not only the total intensity not contained in the Mössbauer peak but also its distribution onto the remaining spectrum. Here, Visscher visscher  derived an analytical result for the observable spectrum at zero temperature in a three dimensional crystal. He also performed numerical calculations for finite temperatures. Systems with almost two-dimensional magnetic interactions are known (interaction strength $`10^4`$-times stronger in the plane than between the planes). Now, it is interesting to identify systems with low-dimensional elastic interactions. Here, one could think about ultrathin films sputtered onto the surface of a solid, a liquid, or the boundary layer between two liquids. Thus a mechanical decoupling to the substrate should be achieveable at least in directions parallel to the surface.

In the following section we discuss the theoretical description of these phenomena by deriving and solving a differential equation for the spectrum. Subsequently, in section III, we evaluate the solution in harmonic approximation for the isotropic Debye model in three and two dimensions at zero and finite temperature. In two dimensions at finite temperature we find that there are no recoil-free events. The Mössbauer effect ceases to exist in the strict sense. However there is still a central peak made up of almost recoil-free events, a Quasi-Mössbauer effect. Ways how to identify two dimensional elastic behavior are discussed. In the final section we summarise our paper.

## II Theory

Starting from the absorption spectrum of a nucleus inside a rigid crystal lattice $`S_{rig}(\omega _\gamma )`$, we write the observed spectrum in an elastic lattice $`S(\omega _\gamma )`$ as a convolution of the first with a shift spectrum: $`S=IS_{rig}`$. Here, $`\omega _\gamma `$ is the circular frequency of the converted photon. The shift spectrum $`I(\omega )`$ comprises all the information on the modifications of the observed spectrum $`S(\omega _\gamma )`$ due to the dynamics of the lattice. Up to an overall factor $`I(\omega )`$ is the spectrum one would observe, if the nucleus’ line width was exactly zero. A normalisation to one $`𝑑\omega I(\omega )=1`$ guarantees the conservation of the integrated intensity: $`𝑑\omega S(\omega )=𝑑\omega S_{rig}(\omega )`$.

In a crystal with $`N`$ Bravais lattice sites, $`A`$ atoms per elementary cell, and where the atoms are elongated in $`d^{}`$ dimensions there are $`M=NAd^{}`$ vibrational modes with circular frequencies $`\mathrm{\Omega }_J`$ with $`J\{1,2,\mathrm{},M\}`$. In general, the dimensionality of the crystal $`d`$ is the same as that of the elongations of the atoms $`d^{}`$ which does not exclude the possibility of doing calculations with $`dd^{}`$.

The state of the crystal before (i) and after (f) the emission or absorption process can be described by a set of $`M`$ numbers each $`\{n\}^{i,f}=\{n_J^{i,f}:J\{1,2,\mathrm{},M\}\}`$. The set of all possible changes of the state is constructed by combining every possible initial state with every possible final state (direct product). The different changes contribute to the shift spectrum $`I(\omega )`$ with different weights $`A(\{n\}^i\{n\}^f)`$. The shift spectrum can be expressed as the weighted mean over all changes conserving the total energy:

$`I(\omega )={\displaystyle \underset{\{n\}^i}{}}{\displaystyle \underset{\{n\}^f}{}}\delta [\mathrm{}\omega (E^fE^i)]A(\{n\}^i\{n\}^f).`$ (1)
In a harmonic crystal the total energy is given by:

$`E^{i,f}=\mathrm{}{\displaystyle \underset{J}{}}\left(n_J^{i,f}+{\displaystyle \frac{1}{2}}\right)\mathrm{\Omega }_J.`$ (2)
Due to the linearity of the above expression, the energy difference in equation (1) is determined exclusively by the change

$`m_J=n_J^fn_J^i`$ (3)
of the occupational numbers in the different modes but not from the state $`\{n\}^i`$ of the crystal before the absorption or emission process. The weight $`A(\{n\}^i\{n\}^f)`$ is equal to the absolute probability for the transition from the state $`\{n\}^i`$ to the state $`\{n\}^f`$. It is given by the product of the probability $`A(\{n\}^i)`$ of the crystal to be in a state $`\{n\}^i`$ before the process and the conditional probability $`A_{\{n\}^i}(\{n\}^f)`$ to reach a final state $`\{n\}^f`$ coming from an initial state $`\{n\}^i`$:

$`A(\{n\}^i\{n\}^f)=A(\{n\}^i)A_{\{n\}^i}(\{n\}^f).`$ (4)
Substitution of equations (2) and (4) into equation (1) leads to:

$`I(\omega )`$ $`=`$ $`{\displaystyle \underset{\{m\}}{}}\delta \left(\mathrm{}\omega \mathrm{}{\displaystyle \underset{J}{}}m_J\mathrm{\Omega }_J\right){\displaystyle \underset{\{n\}^i}{}}A(\{n\}^i)A_{\{n\}^i}(\{n\}^f),`$ (5)
where the elements of $`\{n\}^f`$ are to be replaced in accordance with equation (3). For a harmonic crystal in thermal equilibrium before the elementary process, the terms belonging to a certain mode are independent from the configuration of the rest of the system. After defining the thermally averaged – i.e., averaged with the thermal weight $`P_{n_J^i}^J`$ for having $`n_J^i`$ phonons in mode $`J`$ before the process 222In thermal equilibrium: $`A(\{n\}_i)=_JP_{n_J^i}^J`$. – probability for the production of $`m_J`$ phonons in mode $`J`$:

$`A_{m_J}^J={\displaystyle \underset{{}_{i}{}^{}n_{J}^{}}{}}P_{{}_{i}{}^{}n_{J}^{}}^JA_{{}_{i}{}^{}n_{J}^{}_in_J+m_J}^J,`$ (6)
where $`A_{{}_{i}{}^{}n_{J}^{}_in_J+m_J}^J`$ is the probability for a change of the phonon number in mode $`J`$ from $`{}_{i}{}^{}n_{J}^{}`$ to $`{}_{i}{}^{}n_{J}^{}+m_J`$ 333The probability for the change of two configurations of the crystal into another is given by the product of the probabilities for the respective changes in all the modes: $`A_{\{n\}^i}(\{n\}^f)=_JA_{n_J^in_J^i+m_J}^J`$., the shift spectrum can be reexpressed as:

$`I(\omega )`$ $`=`$ $`{\displaystyle \underset{\{m_J\}}{}}\delta \left(\mathrm{}\omega \mathrm{}{\displaystyle \underset{J}{}}m_J\mathrm{\Omega }_J\right){\displaystyle \underset{J^{}}{}}A_{m_J^{}}^J^{}.`$ (7)
The probabilities are normalised to one: $`_{m_J}A_{m_J}^J=1`$. Consider the shift spectrum $`I^{(J)}(\omega )`$ for only $`J`$ active modes. The shift spectrum $`I^{(0)}(\omega )`$ for a nucleus in an entirely rigid crystal lattice is given by a Dirac-$`\delta `$-distribution: $`\delta (\omega )`$. The shift spectrum $`I^{(M)}(\omega )`$ with all modes in action is identical to the actual shift spectrum $`I(\omega )`$. From equation (7) we get a recursion equation for the partial shift spectra:

$`I^{(J)}(\omega )={\displaystyle \underset{m_J}{}}A_{m_J}^JI^{(J1)}(\omega m_J\mathrm{\Omega }_J).`$ (8)
An analogous equation is valid for the changes $`\mathrm{\Delta }I^{(J)}=I^{(J)}(\omega )I^{(J1)}(\omega )`$ in connection with the shift amplitudes $`a_{m_J}^J=A_{m_J}^J\delta _{m_J,0}`$:

$`\mathrm{\Delta }I^{(J)}(\omega )={\displaystyle \underset{m_J}{}}a_{m_J}^JI^{(J1)}(\omega m_J\mathrm{\Omega }_J).`$ (9)
The shift amplitudes $`a_{m_J}^J`$ differ from the probabilities $`A_{m_J}^J`$ only for $`m_J=0`$. There the respective shift amplitude is smaller by one. For example $`A_{m_J}^J=0`$ for all $`m_J0`$ and $`A_0^J=1`$ signifies that no phonons are created or annihilated in mode $`J`$; so does $`a_{m_J}^J=0`$ for all $`m_J`$. Hence $`A_0^J`$ is a meassure for remaining in the previous state and $`a_0^J`$ indicates the relative change. Because of the additional Kronecker tensor $`\delta _{i,j}`$ and the normalisation of the probabilities $`A_{m_J}^J`$, the sum over all shift amplitudes amounts to zero: $`_{m_J}a_{m_J}^J=0`$. Writing the shift spectrum $`I^{(J)}(\omega )`$ as a sum over all its changes and replacing the changes $`\mathrm{\Delta }I^{(J)}(\omega )`$ of the shift spectrum by the recursion (9) yields:

$`I^{(J)}(\omega )=I^{(0)}(\omega )+{\displaystyle \underset{J^{}=1}{\overset{J}{}}}{\displaystyle \underset{m_J^{}}{}}a_{m_J^{}}^J^{}I^{(J^{}1)}(\omega m_J^{}\mathrm{\Omega }_J^{}).`$ (10)
In case, several active branches are present, the label of the mode $`J`$ is composed of a branch index $`j`$ and a Bloch vector $`\stackrel{}{k}`$: $`J=(j,\stackrel{}{k})`$. In a transition from a discrete to a continuous description the summation over all modes $`J`$ is replaced by a summation over the different modal branches and an integration over $`k`$-space. The $`d`$-dimensional integral is substituted by a one-dimensional integration over the frequencies $`\mathrm{\Omega }`$ of the modes. As a correction, the density of modes $`\eta _j(\mathrm{\Omega })`$ for the corresponding branch has to be taken into account. In the following, all investigations will be carried out for a single active branch. So, with the volume $`V`$ of the crystal:

$`{\displaystyle \underset{J}{}}V{\displaystyle \frac{d^dk}{(2\pi )^d}}{\displaystyle 𝑑\mathrm{\Omega }\eta (\mathrm{\Omega })}.`$ (11)
With the dispersion relation $`\mathrm{\Omega }(\stackrel{}{k})`$ the density could be given by ashcroftmermin :

$`\eta (\mathrm{\Omega })=V{\displaystyle \delta [\mathrm{\Omega }\mathrm{\Omega }(\stackrel{}{k})]\frac{d^dk}{(2\pi )^d}}.`$ (12)
It can be interpreted as the total number of modes in the infinitessimally small interval $`[\mathrm{\Omega },\mathrm{\Omega }+d\mathrm{\Omega }]`$. The modes are counted in the entire crystal and not only in a unit volume. So, the density of modes grows proportionally to the number $`N`$ of Bravais lattice sites. It is to vanish for frequencies outside the intervall $`[0,\mathrm{\Omega }_M]`$. When connected to the discrete description this must also be the interval in which lie all the modal frequencies $`\mathrm{\Omega }_J`$ with $`J\{1,2,\mathrm{},M\}`$. In the Debye model, $`\mathrm{\Omega }_M`$ is equal to the Debye frequency $`\mathrm{\Omega }_D`$. For simple models, the alternative to marking the shift spectrum $`I^{(J)}(\omega )`$ and the shift amplitudes $`a_{m_J}^J`$ with the mode index $`J`$ is to express them as functions of the modal frequency $`\mathrm{\Omega }_J`$: $`I^{(J)}(\omega )=I(\mathrm{\Omega }_J,\omega )`$ and $`a_{m_J}^J=a_m(\mathrm{\Omega }_J)`$. During the transition from the discrete to the continuous description the discrete parameter $`\mathrm{\Omega }_J`$ is replaced by the continuous $`\mathrm{\Omega }`$. Values at $`\mathrm{\Omega }=\mathrm{\Omega }_J`$ for all $`J\{1,2,\mathrm{},M\}`$ are preserved: $`\mathrm{\Omega }_J,J\{1,2,\mathrm{},M\}\mathrm{\Omega }[0,\mathrm{\Omega }]`$, $`I(\mathrm{\Omega }_J,\omega )I(\mathrm{\Omega },\omega )`$, and $`a_m(\mathrm{\Omega }_J)a_m(\mathrm{\Omega })`$. Carrying out the continuum limit in equation (10) leads to:

$`I(\mathrm{\Omega },\omega )=\delta (\omega )+{\displaystyle _0^\mathrm{\Omega }}𝑑\mathrm{\Omega }^{}\eta (\mathrm{\Omega }^{}){\displaystyle \underset{m}{}}a_m(\mathrm{\Omega }^{})I(\mathrm{\Omega }^{},\omega m\mathrm{\Omega }^{}).`$ (13)
Here – for the sake of simplicity – we anticipate a result that will be derived later on: The shift amplitudes $`a_m(\mathrm{\Omega })`$ for fixed argument are proportional to integer negative powers of the number $`N`$ of Bravais lattice sites: $`a_m(\mathrm{\Omega })N^{|m|}`$ for all $`m0`$ and $`a_0(\mathrm{\Omega })N^1`$. The shift amplitudes occur only combined in a product with the density of modes. As already stated above, the density of modes is proportional to the number $`N`$ of Bravais-lattice sites. So, only the three terms with $`|m|1`$ survive taking the continuum limit. Two of the three shift amplitudes that are left can be expressed as functions of the third: They can be eliminated with the help of the normalisation of the shift amplitudes and because we are assuming thermal equilibrium: $`a_m(\mathrm{\Omega })=e^{m\beta \mathrm{}\mathrm{\Omega }}a_{+m}(\mathrm{\Omega })`$, respectively, where $`\beta =k_BT`$. Partial differentiation of equation (13) with respect to the maximum frequency $`\mathrm{\Omega }`$ together with the above facts leads to a linear, first-order, non-local, partial differential equation:

$`{\displaystyle \frac{}{\mathrm{\Omega }}}I(\mathrm{\Omega },\omega )=\eta (\mathrm{\Omega })a_{+1}(\mathrm{\Omega })\{e^{\beta \mathrm{}\mathrm{\Omega }}I(\mathrm{\Omega },\omega +\mathrm{\Omega })(e^{\beta \mathrm{}\mathrm{\Omega }}+1)I(\mathrm{\Omega },\omega )+I(\mathrm{\Omega },\omega \mathrm{\Omega })\},`$ (14)
with the starting condition $`I(0,\omega )=\delta (\omega )`$. A Fourier transformation from the variable $`\omega `$ to the variable $`\varphi `$ and a subsequent division through the Fourier transformed shift spectrum $`\stackrel{~}{I}(\mathrm{\Omega },\varphi )`$ leads to an ordinary first-order differential equation:

$`{\displaystyle \frac{d}{d\mathrm{\Omega }}}\mathrm{ln}\{\stackrel{~}{I}(\mathrm{\Omega },\varphi )\}=\eta (\mathrm{\Omega })a_{+1}(\mathrm{\Omega })\{e^{\beta \mathrm{}\mathrm{\Omega }}(e^{+i\varphi \mathrm{\Omega }}1)+(e^{i\varphi \mathrm{\Omega }}1)\}.`$ (15)
The actual shift spectrum $`I(\omega )`$ is equal to $`I(\mathrm{\Omega }=\mathrm{\Omega }_M,\omega )`$. Taking into account the Fourier transformed boundary condition $`\stackrel{~}{I}(\varphi =0)=1`$, we find by direct integration:

$`\stackrel{~}{I}(\varphi )=\mathrm{exp}\left\{{\displaystyle _0^{\mathrm{\Omega }_M}}𝑑\mathrm{\Omega }\eta (\mathrm{\Omega })a_{+1}(\mathrm{\Omega })[e^{\beta \mathrm{}\mathrm{\Omega }}(e^{+i\varphi \mathrm{\Omega }}1)+(e^{i\varphi \mathrm{\Omega }}1)]\right\}.`$ (16)
Given the existence of the resulting integrals, the integral in the exponent can be seperated into a $`\varphi `$-dependent and a $`\varphi `$-independent part. A Fourier retransformation yields the shift spectrum $`I(\omega )`$ in momentum representation:

$`I(\omega )=e^F{\displaystyle \frac{d\varphi }{2\pi }e^{i\varphi \omega }\mathrm{exp}\left\{_0^{\mathrm{\Omega }_M}𝑑\mathrm{\Omega }\eta (\mathrm{\Omega })a_{+1}(\mathrm{\Omega })[e^{(+i\varphi \beta \mathrm{})\mathrm{\Omega }}+e^{i\varphi \mathrm{\Omega }}]\right\}}.`$ (17)
The factor $`e^F`$ with $`F=_0^{\mathrm{\Omega }_M}𝑑\mathrm{\Omega }\eta (\mathrm{\Omega })a_{+1}(\mathrm{\Omega })\{e^{\beta \mathrm{}\mathrm{\Omega }}+1\}`$ is equal to the Mössbauer-Lamb fraction. The exponential series is uniformly convergent on every bounded area of the complex plane. So, as long as the magnitude of the exponent is finite, the exponential series can be integrated term by term:

$`I(\omega )=e^F{\displaystyle \underset{\nu =0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{\nu !}}{\displaystyle \frac{d\varphi }{2\pi }e^{i\varphi \omega }\left\{_0^{\mathrm{\Omega }_M}𝑑\mathrm{\Omega }\eta (\mathrm{\Omega })a_{+1}(\mathrm{\Omega })[e^{(+i\varphi \beta \mathrm{})\mathrm{\Omega }}+e^{i\varphi \mathrm{\Omega }}]\right\}^\nu }.`$ (18)
Reinterpretation of the inner integration as an additional Fourier transformation and decomposition of the sum into single terms $`I_\nu (\omega )`$ results in:

$`I_\nu (\omega )={\displaystyle \frac{e^F}{\nu !}}{\displaystyle \frac{_{\varphi \omega }}{2\pi }}\left\{_{\mathrm{\Omega }\varphi }[\eta (\mathrm{\Omega })a_{+1}(\mathrm{\Omega })e^{\beta \mathrm{}\mathrm{\Omega }}]+_{\mathrm{\Omega }+\varphi }[\eta (\mathrm{\Omega })a_{+1}(\mathrm{\Omega })]\right\}^\nu ,`$ (19)
where $`_{\omega +\varphi }`$ denotes a Fourier transformation and where $`(2\pi )^1_{\varphi \omega }`$ is its inverse. The evaluation of the first two terms is independent of the form of the density of modes or the shift amplitude:

$`I_0(\omega )`$ $`=`$ $`e^F\delta (\omega )`$
$`I_1(\omega )`$ $`=`$ $`e^F\left\{\eta (\omega )a_{+1}(\omega )e^{+\beta \mathrm{}\omega }+\eta (+\omega )a_{+1}(+\omega )\right\}.`$ (20)
Terms with $`\nu 2`$ could be reexpressed with the help of the binomial theorem and by noting that products become convolutions after a Fourier transformation. Then one sees that the term with index $`\nu `$ contains the effect of events on the shift spectrum during which $`\nu `$ phonons are converted, i.e., created or annihilated.

## III Application

Now, we would like to evaluate explicitely the shift spectrum $`I(\omega )`$ in harmonic approximation for the isotropic Debye model. The Debye model is based on the assumption that the crystal is a $`d`$-dimensional hypercubic lattice. So, it is a Bravais-crystal in which no optical modes exist. Isotropy leads to degenerate modes on all $`d^{}=d`$ branches. The frequencies follow a linear dispersion relation: $`\mathrm{\Omega }(\stackrel{}{k})=sk`$ with the speed of sound $`s`$ and the absolute value $`k`$ of the wave-number vector $`\stackrel{}{k}`$. Additionally, the integration over the first Brillouin zone is replaced by an integration over the volume of a sphere in $`k`$-space. The size of the sphere is chosen in such a way that it contains $`N`$ allowed wave-number vectors where $`N`$ is the number of ions in the crystal. Evaluation of equation (12) under these assumptions leads to:

$`\eta (\mathrm{\Omega })=d^2{\displaystyle \frac{N}{\mathrm{\Omega }_D}}\left({\displaystyle \frac{\mathrm{\Omega }}{\mathrm{\Omega }_D}}\right)^{d1}\theta (\mathrm{\Omega }_D\mathrm{\Omega })\theta (\mathrm{\Omega }),`$ (21)
with the Debye frequency $`\mathrm{\Omega }_D=sk_D`$. As mentioned above, the Debye frequency replaces the maximal frequency $`\mathrm{\Omega }_M`$. The transition probabilities which are the thermally averaged absolute squares of the transition matrix-elements $`\{n_J^{}+m_J^{}\}|e^{iqx}|\{n_J^{}\}`$ are to be evaluated in a harmonic system. Here, $`q`$ stands for the momentum transfer and $`x`$ for the spatial coordinate along its direction. After some calculations one finds:

$`A_{m_J}^J=e^{+m_J\beta \mathrm{}\mathrm{\Omega }_J/2}\mathrm{exp}\left\{\rho ^J\mathrm{coth}\left({\displaystyle \frac{\beta \mathrm{}\mathrm{\Omega }_J}{2}}\right)\right\}\mathrm{I}_{m_J}\left\{\rho ^J\left[\mathrm{sinh}\left({\displaystyle \frac{\beta \mathrm{}\mathrm{\Omega }_J}{2}}\right)\right]^1\right\}=a_{m_J}^J+\delta _{J,0},`$
(22)
with the parameter:

$`\rho ^J={\displaystyle \frac{1}{N}}{\displaystyle \frac{(\mathrm{}q)^2}{2M_{ion}\mathrm{}\mathrm{\Omega }_J}},`$ (23)
and where $`\mathrm{I}_m`$ stands for the modified Bessel function of the first kind and $`M_{ion}`$ for the mass of the ion whose nucleus absorbs the $`\gamma `$-photon. For small arguments, we could use the lowest order approximation to the Bessel functions: $`\mathrm{I}_m(w)=(w/2)^m/m!+𝒪(w^{m+2})`$. However, the argument of the Bessel function is divergent for a vanishing circular frequency $`\mathrm{\Omega }`$. Due to the periodic boundary conditions there exists a minimal circular frequency $`\mathrm{\Omega }_1`$ which, for a constant volume density of atoms $`n`$, depends on the total number of atoms according to: $`\mathrm{\Omega }_1=\mathrm{\Xi }N^{1/d}`$. The constant of proportionality $`\mathrm{\Xi }`$ depends on the geometry of the lattice and for a $`d`$-dimensional hypercubic lattice is given by $`\mathrm{\Xi }=2\pi sn^{1/d}`$. For large $`N`$, the argument of the Bessel function becomes: $`w=2q^2N^{(2d)/d}/(M_{ion}\beta \mathrm{\Xi }^2)`$. In three dimensions it goes to zero as $`N`$ goes to infinity. In two dimensions it is independent of the number of atoms $`N`$, so for the validity of the lowest order approximation we have to postulate additionally: $`2q^2/(M_{ion}\beta \mathrm{\Xi }^2)<<1`$. In one dimension the argument $`w`$ is not boundend from above and so the lowest order approximation is not adequate. In every case where the lowest order approximation is justified, we have:

$`a_{m_J}^J={\displaystyle \frac{e^{+m_J\beta \mathrm{}\mathrm{\Omega }_J/2}}{|m_J|!}}\left\{{\displaystyle \frac{\rho ^J}{2}}\left[\mathrm{sinh}\left({\displaystyle \frac{\beta \mathrm{}\mathrm{\Omega }_J}{2}}\right)\right]^1\right\}^{|m_J|}N^{|m_J|}`$ (24)
for $`|m_J|1`$ and

$`a_0^J=\rho ^J\mathrm{coth}\left({\displaystyle \frac{\beta \mathrm{}\mathrm{\Omega }_J}{2}}\right)N^1`$ (25)
for $`m_J=0`$. The density of modes $`\eta (\mathrm{\Omega })`$ only contributes one additional factor of $`N`$. So, if $`|m_J|>1`$ the negative powers of the number of ions in the crystal cannot be compensated and the corresponding terms do not contribute. This result has already been used to derive equation (14). For vanishing absolute temperature, the above expansion is always legitimate and one finds:

$`a_{+1}^J=+\rho ^Ja_0^J=\rho ^Ja_1^J=0.`$ (26)
Finally, the transition to the continuum has to be performed as described above. Now, equation (16) can be reexpressed depending on the dimension $`d`$ of the crystal:

$`\stackrel{~}{I}_d(\varphi )=\mathrm{exp}\left\{{\displaystyle \frac{\sigma _d}{\mathrm{\Omega }_D^{d1}}}{\displaystyle _0^{\mathrm{\Omega }_D}}𝑑\mathrm{\Omega }\mathrm{\Omega }^{d2}\left[\mathrm{coth}\left\{{\displaystyle \frac{\beta \mathrm{}\mathrm{\Omega }}{2}}\right\}(\mathrm{cos}\{\varphi \mathrm{\Omega }\}1)i\mathrm{sin}\{\varphi \mathrm{\Omega }\}\right]\right\}`$ (27)
with the parameter: $`\sigma _d=(dq\mathrm{})^2/(2M_{ion}\mathrm{}\mathrm{\Omega }_D)`$. The aforementioned separation into a $`\varphi `$-dependent and a $`\varphi `$-independent part without leaving a non-integrable pole in the $`\varphi `$-dependent part is possible in three dimensions at any temperature and in two dimensions for zero temperature. The analytic expression for the exponent $`F_d`$ of the Mössbauer-Lamb fraction is given by:

$`F_d={\displaystyle \frac{\sigma _d}{\mathrm{\Omega }_D^{d1}}}{\displaystyle _0^{\mathrm{\Omega }_D}}𝑑\mathrm{\Omega }\mathrm{\Omega }^{d2}\mathrm{coth}\left\{{\displaystyle \frac{\beta \mathrm{}\mathrm{\Omega }}{2}}\right\}`$ (28)
which for zero temperature is equal to $`F_d^0(\mathrm{\Omega }_D)=\sigma _d/(d1)`$ and for non-zero temperature and in three dimensions reads:

$`F_3=\sigma _3\left[{\displaystyle \frac{\mathrm{ln}(1e^{+\beta \mathrm{}\mathrm{\Omega }_D})}{\beta \mathrm{}\mathrm{\Omega }_D}}{\displaystyle \frac{\mathrm{dilog}(e^{+\beta \mathrm{}\mathrm{\Omega }_D})+\mathrm{dilog}(1e^{+\beta \mathrm{}\mathrm{\Omega }_D})}{(\beta \mathrm{}\mathrm{\Omega }_D)^2}}+{\displaystyle \frac{\pi ^2}{6(\beta \mathrm{}\mathrm{\Omega }_D)^2}}\right]`$ (29)
where the dilogarithm is defined as in chapter 27.7 of reference abramowitzstegun . For the subset of cases investigated up to the present, the exponent is bounded and the Fourier retransformation can be carried out term by term for the exponential series. The integrated contributions $`\stackrel{~}{I}_{d,\nu }(\varphi =0)`$ follow a Poissonian distribution with a mean value equal to $`F_d`$. So, this parameter determines how fast the series converges. For zero temperature the retransformations for all the terms of the series can be carried out exactly:

$`I_{3,\nu }^0(\omega )={\displaystyle \frac{e^{\sigma _3/2}}{\mathrm{\Omega }_D}}{\displaystyle \underset{\mu _1=0}{\overset{\nu }{}}}{\displaystyle \underset{\mu _2=0}{\overset{\mu _1}{}}}{\displaystyle \frac{(1)^{\mu _1}\sigma _{3}^{}{}_{}{}^{\nu }}{\nu !(2\nu \mu _21)!}}\left({\displaystyle \genfrac{}{}{0pt}{}{\nu }{\mu _1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{\mu _1}{\mu _2}}\right)\theta \left({\displaystyle \frac{\omega }{\mathrm{\Omega }_D}}\mu _1\right)\left({\displaystyle \frac{\omega }{\mathrm{\Omega }_D}}\mu _1\right)^{2\nu \mu _21}`$ (30)
(see figure 1) and

$`I_{2,\nu }^0(\omega )={\displaystyle \frac{e^{\sigma _2}}{\mathrm{\Omega }_D}}{\displaystyle \underset{\mu =0}{\overset{\nu }{}}}{\displaystyle \frac{(1)^\mu \sigma _{2}^{}{}_{}{}^{\nu }}{\nu !(\nu 1)!}}\left({\displaystyle \genfrac{}{}{0pt}{}{\nu }{\mu }}\right)\theta \left({\displaystyle \frac{\omega }{\mathrm{\Omega }_D}}\mu \right)\left({\displaystyle \frac{\omega }{\mathrm{\Omega }_D}}\mu \right)^{\nu 1}`$ (31)
(see figure 2). The integrations in the three dimensional case and at finite temperatures can be carried out after approximating the hyperbolic cotangent by:

$`\mathrm{coth}\left\{{\displaystyle \frac{\beta \mathrm{}\mathrm{\Omega }}{2}}\right\}1+2{\displaystyle \frac{e^{\beta \mathrm{}\mathrm{\Omega }/2}}{\beta \mathrm{}\mathrm{\Omega }}}.`$ (32)
An expansion into partial fractions leads to (see figure 3):

$`I_{3,\nu }(\omega ){\displaystyle \underset{\mu _1=0}{\overset{\nu }{}}}{\displaystyle \underset{\mu _2=0}{\overset{\mu _1}{}}}{\displaystyle \underset{\lambda _1=1}{\overset{\nu }{}}}`$ $`\mathrm{\Lambda }_{1,\lambda _1}^{}`$ $`(\omega ^{})^{\lambda _11}e^{+\beta \mathrm{}\mathrm{\Omega }_D\omega ^{}}\theta (\omega ^{})+`$
$`+{\displaystyle \underset{\lambda _2=1}{\overset{2(\nu \mu _2)}{}}}`$ $`\mathrm{\Lambda }_{2,\lambda _2}^{}`$ $`(+\omega ^{})^{\lambda _21}\theta (+\omega ^{})+`$
$`+{\displaystyle \underset{\lambda _3=1}{\overset{\nu \mu _2}{}}}`$ $`\mathrm{\Lambda }_{3,\lambda _3}^{}`$ $`(+\omega ^{})^{\lambda _31}e^{\beta \mathrm{}\mathrm{\Omega }_D\omega ^{}}\theta (+\omega ^{})`$ (33)
with the factors independent of $`\omega ^{}=2\mu _2\mu _1+\omega /\mathrm{\Omega }_D`$ defined as:

$`\mathrm{\Lambda }_{(1,2,3),\lambda }^{}=e^F{\displaystyle \frac{1}{\mathrm{\Omega }_D}}{\displaystyle \frac{\sigma _{3}^{}{}_{}{}^{\nu }}{\nu !}}\left({\displaystyle \frac{\beta \mathrm{}\mathrm{\Omega }_D}{2}}\right)^{2\nu \mu _1}\left({\displaystyle \genfrac{}{}{0pt}{}{\nu }{\mu _1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{\mu _1}{\mu _2}}\right)e^{\mu _2\beta \mathrm{}\mathrm{\Omega }/2}{\displaystyle \frac{\mathrm{\Lambda }_{(1,2,3),\lambda }}{\mathrm{\Omega }_{D}^{}{}_{}{}^{4\nu \mu _1\lambda }}}{\displaystyle \frac{(1)^\lambda }{(\lambda 1)!}}`$ (34)
and where the coefficients $`\mathrm{\Lambda }_{(1,2,3),\lambda }`$ are given by:

$`\mathrm{\Lambda }_{1,\lambda _1}`$ $`=`$ $`{\displaystyle \frac{1}{(\lambda _11)!}}{\displaystyle \frac{d^{\nu \lambda _1}}{d\varphi ^{\nu \lambda _1}}}|_{i\varphi =+\beta \mathrm{}}{\displaystyle \frac{p^{\mu _1\mu _2}(\varphi )}{(i\varphi )^{2(\nu \mu _2)}(i\varphi \beta \mathrm{})^{\nu \mu _2}}}`$
$`\mathrm{\Lambda }_{2,\lambda _2}`$ $`=`$ $`{\displaystyle \frac{1}{(\lambda _21)!}}{\displaystyle \frac{d^{2(\nu \mu _2)\lambda _2}}{d\varphi ^{2(\nu \mu _2)\lambda _2}}}|_{\varphi =0}{\displaystyle \frac{p^{\mu _1\mu _2}(\varphi )}{(i\varphi \beta \mathrm{})^{\nu \mu _2}(+i\varphi \beta \mathrm{})^\nu }}`$
$`\mathrm{\Lambda }_{3,\lambda _3}`$ $`=`$ $`{\displaystyle \frac{1}{(\lambda _31)!}}{\displaystyle \frac{d^{\nu \mu _2\lambda _3}}{d\varphi ^{\nu \mu _2\lambda _3}}}|_{i\varphi =\beta \mathrm{}}{\displaystyle \frac{p^{\mu _1\mu _2}(\varphi )}{(i\varphi )^{2(\nu \mu _2)}(+i\varphi \beta \mathrm{})^\nu }}`$ (35)
with the polynomial in $`\varphi `$:

$`p(\varphi )=i\varphi ^3(\beta \mathrm{}\mathrm{\Omega }_D+e^{\beta \mathrm{}\mathrm{\Omega }_D/2})\varphi ^2\beta \mathrm{}(1e^{\beta \mathrm{}\mathrm{\Omega }_D/2}/2)i\varphi (\beta \mathrm{})^2(\beta \mathrm{}\mathrm{\Omega }_D/4)(\beta \mathrm{})^3/4`$
For the two-dimensional case at finite temperatures, the separation into a $`\varphi `$-dependent and a $`\varphi `$-independent part is no longer possible without leaving behind a non-integrable pole. Actually, the Mössbauer-Lamb fraction $`e^F`$ is exactly equal to zero because its exponent diverges. Hence, there is no recoil-free emission in this case. For further evaluation, we start by separating off the shift spectrum for zero temperature from that for finite temperature:

$`\stackrel{~}{I}_2(\varphi )=\stackrel{~}{I}_2^0(\varphi )\mathrm{exp}\left\{2{\displaystyle \frac{\sigma _2}{\mathrm{\Omega }_D}}{\displaystyle _0^{\mathrm{\Omega }_D}}𝑑\mathrm{\Omega }{\displaystyle \frac{e^{\beta \mathrm{}\mathrm{\Omega }}}{1e^{\beta \mathrm{}\mathrm{\Omega }}}}(\mathrm{cos}\{\varphi \mathrm{\Omega }\}1)\right\}`$ (37)
With an approximation analogous to that in equation (32) one finds:

$`\stackrel{~}{I}_2(\varphi )\stackrel{~}{I}_2^0(\varphi ){\displaystyle \frac{(\beta \mathrm{}\mathrm{\Omega }_D/2)^2}{(\beta \mathrm{}\mathrm{\Omega }_D/2)^2+(\varphi \mathrm{\Omega }_D)^2}}\mathrm{exp}\left\{2{\displaystyle \frac{\sigma _2}{\beta \mathrm{}\mathrm{\Omega }_D}}[E_1(\beta \mathrm{}\mathrm{\Omega }_D/2)\mathrm{Re}\{E_1[(\beta \mathrm{}/2\pm i\varphi )\mathrm{\Omega }_D]\}]\right\},`$
where the exponential integral function $`E_1`$ is defined as in chapter 5 of reference abramowitzstegun :

$`E_1(z)={\displaystyle _z^{\mathrm{}}}𝑑t{\displaystyle \frac{e^t}{t}}.`$ (39)
The following integrations cannot be performed analytically, due to the exponential integral function $`E_1`$. For the two dimensional case we had already postulated that $`\sigma _2`$ is to be small compared to one. In the physically interesting case of absolute temperatures $`T`$ small compared to the Debye temperature $`\mathrm{\Theta }_D=\mathrm{}\mathrm{\Omega }_D/k_B`$, the exponent as such becomes extremely small and the exponential function can be approximated by one. In the low-temperature approximation it remains to perform the Fourier retransformation of the power of the Lorentz distribution. However, the asymptotic form of this function for large $`\varphi `$ is given by $`[\beta \mathrm{}/(2\varphi )]^{2\sigma _2/(\beta \mathrm{}\mathrm{\Omega }_D)}`$. Due to the small absolute values of the exponent this function decays slower than $`\varphi ^1`$ which renders the Fourier retransformation difficult. Let us introduce an additional factor which ensures the convergence of the integral under investigation. If this manoeuvre is performed in a suitable way, the observable spectrum $`S(\omega _\gamma )`$ can be obtained directly. When choosing a Lorentz distribution as a model for the natural line shape of the nucleus, the observabel spectrum is just the convolution of the former with the shift spectrum $`I(\omega )`$ in momentum representation. This corresponds to a multiplication of the Fourier transformed shift spectrum $`\stackrel{~}{I}(\varphi )`$ with the Fourier transformed Lorentz distribution, i.e., a Laplace distribution also known as a doubly-exponential distribution. What remains to be executed is:

$`S_2(\omega )S_0I_2^0(\omega ){\displaystyle \frac{d\varphi }{2\pi }e^{i\varphi (\omega \omega _0)}e^{\zeta |\varphi |}\left(\frac{(\beta \mathrm{}\mathrm{\Omega }_D/2)^2}{(\beta \mathrm{}\mathrm{\Omega }_D/2)^2+(\varphi \mathrm{\Omega }_D)^2}\right)^{\frac{\sigma _2}{\beta \mathrm{}\mathrm{\Omega }_D}}},`$ (40)
with the resonance frequency of the nucleus $`\omega _0`$, the natural line-width $`\zeta `$ (half-width half-maximum) and the integrated intensity $`S_0`$. The integral can be interpreted as a Laplace transformation and yields (see erdelyi  and figure 5):

$`S_2(\omega ){\displaystyle \frac{\beta \mathrm{}S_0}{2\pi }}I_2^0(\omega )\mathrm{Re}\left\{\mathrm{L}_{\frac{1}{2}\frac{\sigma _2}{\beta \mathrm{}\mathrm{\Omega }_D}}\left({\displaystyle \frac{\beta \mathrm{}[\zeta i(\omega \omega _0)]}{2}}\right)\right\},`$ (41)
with:

$`\mathrm{L}_l(z)=2^{l1}\sqrt{\pi }\mathrm{\Gamma }\left(l+{\displaystyle \frac{1}{2}}\right)z^l[\mathrm{H}_l(z)\mathrm{Y}_l(z)],`$ (42)
where $`\mathrm{H}_l(z)`$ denotes the Struve function (see chapter 12 of reference abramowitzstegun ), $`\mathrm{Y}_l(z)`$ stands for a Bessel function of the second kind also known as Weber function (see chapter 9 of reference abramowitzstegun ), and $`\mathrm{\Gamma }(z)`$ represents the $`\mathrm{\Gamma }`$-function. In the limit of vanishing line width $`\zeta `$, the real part of the argument of $`\mathrm{L}_l(z)`$ is equal to zero. Now, the approximation to what shall be called the Quasi-Mössbauer line of the shift spectrum is the convolution of the zeroth order term of the shift spectrum for zero temperature $`I_{2,\nu =0}^0(\omega )`$ with the line shape given by equation (41) with $`S_0`$, $`\omega _0`$, and $`\zeta `$ set to zero (see figure 4):

$`I_{2,\nu =0}^T(\omega ){\displaystyle \frac{e^{\sigma _2}}{\mathrm{\Omega }_D}}{\displaystyle \frac{\beta \mathrm{}\mathrm{\Omega }_D}{2\pi }}\mathrm{Re}\left\{\mathrm{L}_{\frac{1}{2}\frac{\sigma _2}{\beta \mathrm{}\mathrm{\Omega }_D}}\left({\displaystyle \frac{\beta \mathrm{}[\zeta i(\omega \omega _0)]}{2}}\right)\right\}.`$ (43)
## IV Results

In three dimensions there is largely agreement between the present work and visscher . Here, a deeper insight has been permitted due to the presentation of the results in an analytic form. The investigations have been extended to incorporate the two dimensional case at zero and at finite temperature respectively. The one-dimensional case at zero temperature is also covered but not investigated any further here.

In the two dimensional case at finite absolute temperatures, the Mössbauer-Lamb fraction according to its definition (17) vanishes exactly. No recoil-free absorption or emission takes place and strictly speaking, there is no Mössbauer line in the spectra. The mathematical investigation carried out here in this case is only applicable for small values of the model parameter $`\sigma _2`$. But even then the execution of the calculations was not possible in the same way as for the other cases (three dimensions and/or zero absolute temperature). An alternative way had to be pursued. Although a Mössbauer line does not exist, a central peak is found. Already in the shift-spectrum representation this Quasi-Mössbauer peak shows a finite width. The width depends on the model parameter $`\sigma _2`$ and the temperature $`T`$.

The disappearance of the Mössbauer-Lamb fraction in two dimensions and at finite absolute temperature signifies that no entirely recoil-free events are possible and that there cannot be a Mössbauer peak. Its replacement, the Quasi-Mössbauer line must be based on quasi-recoil-free processes. These are events during which only very little energy is transformed into lattice vibrations. This could be due to the predominant creation and annihilation of low-energy phonons. Another more general possibility is based on a dominance of multi-phonon events with small net energy. In the two-dimensional case, as opposed to the three-dimensional, the density of modes does not drop rapidly enough for small absolute values of the circular frequency in order to compensate the pole of the shift amplitudes. This indicates that the first of the two explanations given above is the decisive one but does not rule out the second. The different mechanism of formation in two dimensions leads to a central peak of finite width depending on the temperature. In the representation as observable spectrum, i.e., after the convolution of the shift spectrum with the natural line shape, the Mössbauer and the Quasi-Mössbauer line appear qualitatively similar (see figure 4). Quantitatively they are not. They differ in peak height, line width, and integrated intensity. Given the qualitative similarity, the calculations would have to be reperformed with a more realistic than the Debye model in order to obtain more reliable quantitative results.

While the phonon wings might render the determination of the central line’s integrated intensity cumbersome, they provide the best means to identify two-dimensional elastic behavior (see figure 5). For narrow spectra, in three dimensions, the phonon wing passes through a relative minimum almost reaching down to zero close to the central peak. In two dimensions a nearly horizontal passage can be seen. This means also that the phonon wings for positive and negative circular frequencies are continuously connected at zero frequency in three dimensions (see also figure 3). In two dimensions there is a step (see figure 5).

Room for extension of the present work lies in the investigation of slowly decaying spectra in two dimensions.

## Acknowledgements

Help by R. J. Jelitto is gratefully acknowledged.


format: text

text:
# THE 𝑟-PROCESS IN SUPERNOVAE: IMPACT OF NEW MICROSCOPIC MASS FORMULAS

## 1 INTRODUCTION

The origin of the rapid neutron-capture ($`r`$-process) nuclei is still a mystery. One of the underlying difficulties is that the astrophysical site (and consequently the astrophysical conditions) in which the r-process takes place has not been identified. Although some scenarios such as “neutrino-driven winds” from nascent neutron stars (Woosley et al., 1994; Takahashi et al., 1994; Qian & Woosley, 1996; Cardall & Fuller, 1997; Otsuki et al., 2000; Wanajo et al., 2001; Thompson, Burrows, & Meyer, 2001), “prompt supernova explosions” from small iron cores (Sumiyoshi et al., 2001) or from O-Ne-Mg cores (Wanajo et al., 2003), and “neutron star mergers” (Freiburghaus, Rosswog, & Thielemann, 1999) show some promise, each of them faces severe problems and cannot at the present time be called for to explain the production and galactic enrichment of the $`r`$-process nuclei observed in nature.

Another underlying difficulty is due to the uncertainties in the theoretical predictions of nuclear data far from the $`\beta `$-stability, for which essentially no experimental data exist. In particular, mass predictions for neutron-rich nuclei play a key role since they affect all the nuclear quantities of relevance in the $`r`$-process, namely the neutron capture, photodisintegration and $`\beta `$-decay rates, as well as the fission probabilities . Although most of the recent mass formulas show fits to experimental masses of similar quality (characterized by an rms error of about $`0.7`$ MeV), the mass extrapolations far from the valley of $`\beta `$-stability can differ from each other quite significantly (for a recent review, see Lunney, Pearson, & Thibault, 2003). Recently, Hartree-Fock mass formulas with fully microscopic approaches have been constructed (Goriely, Tondeur, & Pearson, 2001; Samyn et al., 2002; Goriely et al., 2002; Samyn et al., 2003; Goriely et al., 2003). The latest Hartree-Fock-Bogoliubov formula, labeled HFB-2 up to HFB-7(Goriely et al., 2002; Samyn et al., 2003; Goriely et al., 2003), are among the most accurate mass formulas, predicting the 2135 measured masses with a root-mean-square error around 0.670 MeV for nuclei with $`N,Z8`$, i.e with the same accuracy or even better than the one obtained with droplet-like mass formulas (Hilf et al., 1976; Möller et al., 1995). A brief comparison of these mass models is given in § 2.

The purpose of this study is to examine the effects of the newly-derived microscopic mass formula on the $`r`$-process nucleosynthesis and analyse to what extent a solar-like distribution can be obtained. Most of previous works devoted to the study of nuclear mass formulas and their impact on the $`r`$-process nucleosynthesis were based on site-independent approaches, the so-called canonical model assuming a constant temperature and neutron number density during the neutron-capture phase (Goriely & Arnould, 1996; Cowan et al., 1999; Goriely & Clerbaux, 1999; Goriely & Arnould, 2001; Schatz et al., 2002). This approach is understandable since the astrophysical $`r`$-process site has not been unambiguously identified. Another approach consists in considering promising sites and in modifying some of the relevant characteristics to force a successful $`r`$-process. In this case, the free parameter space is usually reduced and an easier analysis can be performed. In addition, some works based on such site-specific approaches have demonstrated that the abundance distributions are sensitive to the astrophysical conditions adopted, in particular during the “freezeout” phase (Surman et al., 1997; Freiburghaus et al., 1999; Surman & Engel, 2001; Wanajo et al., 2002).

We adopt, here, for the physical conditions the semi-realistic astrophysical model of the “prompt supernova explosion” from the collapsing O-Ne-Mg core by Wanajo et al. (2003) (§ 3). The $`r`$-process nucleosynthesis with the HFB-2 mass formula in each outgoing mass trajectory is then calculated with a nuclear reaction network code. The mass-averaged yields over the mass shells relevant for the $`r`$-process is compared with the $`r`$-process abundance patterns in the solar system and in an extremely metal-poor star (CS 22892-052), as well as with those obtained with other mass formulas, more specifically the HFB-7 mass prediction and the extensively used droplet formulas of Hilf et al. (1976) and Möller et al. (1995) (§ 4). Uncertainties stemming from the astrophysical conditions and extra nuclear ingredients are discussed in § 5. A summary follows in § 6.

## 2 MICROSCOPIC MASS MODELS

Among the ground state properties, the atomic mass is obviously the most fundamental quantity and influences the $`r`$-process abundance predictions mainly through the $`(n,\gamma )(\gamma ,n)`$ competition taking place in the neutron-rich region. The calculation of the reaction rates also requires the knowledge of other ground state properties, such as the deformation, density distribution, single-particle level scheme, pairing force and shell correction energies. Recently, impressive progress has been made experimentally. This situation results largely from recent measurements with Penning-trap or Schottky spectrometers which have enlarged the region of known masses, in particular towards the neutron-deficient side of the valley of nuclear stability (for a recent review, see Lunney, Pearson, & Thibault, 2003). The new Atomic Mass Evaluation (Audi & Wapstra, 2001) contains 2214 measured masses, i.e 250 more than the 1995 one (Audi & Wapstra, 1995). A more accurate mass determination is also available now for about 132 nuclides originally included in the 1995 compilation. Out of the 382 new experimental masses, 337 are located in the proton-rich region of the nuclear chart and only 45 in the neutron rich region. As far as nuclei directly involved in the $`r`$-process are concerned, almost no experimental mass data exist, and theory must fill the gap.

Attempts to estimate nuclear masses go back to the liquid-drop Weizsäcker mass formula. Improvements to this original model have been brought little by little, leading to the development of macroscopic-microscopic mass formulas, such as the droplet model (DM) (e.g Hilf et al., 1976) and the ‘finite-range droplet model’ (FRDM) of Möller et al. (1995). In this framework, the macroscopic contribution to the masses and the microscopic corrections of phenomenological nature are treated independently, both parts being connected solely through a parameter fit to experimental masses. Despite the great empirical success of these formulas (e.g. FRDM fits the 2135 $`Z8`$ experimental masses (Audi & Wapstra, 2001) with an rms error of 0.676 MeV), it suffers from major shortcomings, such as the incoherency of the link between the macroscopic part and the microscopic correction, the instability of the mass prediction to different parameter sets, or the instability of the shell corrections. As a consequence, its reliability when extrapolating far from experimentally known masses is severely limited. Although these models have been extensively used for astrophysics applications (Woosley et al., 1994; Takahashi et al., 1994; Sumiyoshi et al., 2001; Wanajo et al., 2001, 2002, 2003), there is an obvious need to develop a mass model that is more closely connected to the basic nuclear interaction properties.

A new major progress has been achieved recently within the Hartree-Fock method (Goriely, Tondeur, & Pearson, 2001; Samyn et al., 2002, 2003; Goriely et al., 2002, 2003). It is now demonstrated that this microscopic approach, making use of a Skyrme force fitted to essentially all the mass data, is not only feasible, but can successfully compete with the most accurate droplet-like formulas available nowadays in the reproduction of measured masses. This holds true not only when the pairing force is described in the BCS approximation, but also when the Bogoliubov method is adopted (HFB model), which has the advantage of ensuring the self-consistency of the treatment of the nuclear single-particle and pairing properties. These large-scale HFB calculations are based on the conventional 10-parameter Skyrme force along with a 4-parameter $`\delta `$-function pairing force. The Skyrme and pairing parameters are derived from a fit to the full data set of 2135 measured masses of nuclei with $`Z,N8`$, leading to an rms error of the order of 0.674 MeV for the HFB-2 mass table.

Despite the success of the HFB-2 mass formula, a series of studies of possible modifications to the basic force model and to the method of calculation was initiated all within the HFB framework in order to test the reliability of the mass extrapolations (Samyn et al., 2003; Goriely et al., 2003). For this reason, a set of additional 5 new mass tables, referred to as HFB-3 to HFB-7 were designed to analyse the sensitivity of the mass fit and extrapolations to the prescription used for the center-of-mass correction, the density-dependence of the pairing interaction and the amplitude of the effective nucleon mass. The new mass fits are of the same quality as the HFB-2 mass predictions. In addition, it is found that globally the extrapolations out to the neutron-drip line of all these different HFB mass formulas are essentially equivalent. Figure 1 (lower panel) compares the HFB-2 and HFB-7 masses for all nuclei with $`8Z110`$ lying between the proton and neutron driplines. Although HFB-2 and HFB-7 are obtained with significantly different Skyrme forces (in particular, HFB-2 is characterized by an density-independent pairing force and an effective isoscalar mass $`M_s^{}=1.04`$, while HFB-7 has an density-dependent pairing force and $`M_s^{}=0.8`$), deviations smaller than about 2 MeV are obtained for nuclei with $`Z82`$.

In contrast, more discrepancies are seen between HFB and droplet-type masses (Figure 1), especially for superheavy nuclei. For lighter species, the mass differences remain below some 5-10 MeV, although significantly different shell and deformation effects are predicted. Most particularly, the HFB mass formulas show a weaker neutron-shell closure close to the neutron drip line with respect to droplet-like models as FRDM (for a more detailed discussion, see Goriely et al., 2002). Future improved microscopic predictions, in particular in mean field, but also shell model approaches, as well as future experiments will hopefully shed light on the behavior of the $`N=82`$ and 126 neutron shell closures at large isospins, in a similar way as for the $`N=20`$, 28, and 50 magic numbers (Lunney, Pearson, & Thibault, 2003). These major differences in the nuclear structure properties may affect the nucleosynthesis predictions, as shown in § 4.

## 3 SUPERNOVA MODEL AND THE $`r`$-PROCESS

Recent comprehensive spectroscopic studies of extremely metal-poor stars (Hill et al., 2002; Cowan et al., 2002; Sneden et al., 2003), as well as chemical evolution studies (Ishimaru & Wanajo, 1999; Ishimaru et al., 2004; Argast et al., 2003), in the Galactic halo suggest that the $`r`$-process is a primary nucleosynthesis process and that the astrophysics site could be associated with core-collapse supernovae. In fact, most of the recent studies on the nucleosynthesis $`r`$-process have been based on the “neutrino-driven wind” scenario, in which the free nucleons accelerated by the intense neutrino flux near the neutrino sphere of a core-collapse supernova assemble to heavier nuclei (Woosley et al., 1994; Qian & Woosley, 1996; Cardall & Fuller, 1997; Otsuki et al., 2000; Wanajo et al., 2001; Thompson, Burrows, & Meyer, 2001).

In the present study, we use the parametrized model of the “prompt supernova explosion” from an $`810M_{}`$ progenitor star (with a $`1.38M_{}`$ O-Ne-Mg core) by Wanajo et al. (2003). The reason is that this model leads to $`r`$-abundance distributions that have been shown to be relatively similar to the solar distribution, at least if an artificial enhancement of the shock-heating energy is assumed. In addition, this scenario does not suffer from the problematic overproduction of $`A90`$ nuclei seen in the neutrino-driven wind model (Woosley et al., 1994; Wanajo et al., 2001). Finally, the physical conditions in hydrodynamical (i.e., prompt) explosions may not be significantly affected by highly uncertain neutrino heating (or cooling) processes (Hillebrandt, Nomoto, & Wolff, 1984) as they are in the neutrino wind model. The site modeling is consequently easier and the final abundance predictions less sensitive to the neutrino physics. All these reasons justify the present study and the comparison of the nucleosynthesis results with the solar $`r`$-abundance distribution. It is obvious, however, that a similar study based on the neutrino-driven wind scenario is needed since the prompt explosion scenario also faces severe problems that need to be resolved (Wanajo et al., 2003). In particular, it should be recalled here that in consistent realistic prompt explosion models, only a weak explosion is obtained and no $`r`$-processing. An artificial enhancement of the shock energy is needed to drive an energetic explosion and provide the necessary conditions for a successful $`r`$-process.

Wanajo et al. (2003) performed core-collapse simulations with a one-dimension, Newtonian hydrodynamic code. The highly neutronized ejecta ($`Y_{e0}`$, initial electron mole fraction, $`0.140.20`$) in the (artificial) energetic explosion (model Q6) is subject to a strong production of $`r`$-process nuclei. In this study, we use the mass trajectories of the ejected material from $`M_{\mathrm{ej}}=0.08M_{}`$ (zone number 60, $`Y_{e0}=0.45`$) to $`0.31M_{}`$ (zone number 105, $`Y_{e0}=0.14`$) in mass coordinate (the surface of the O-Ne-Mg core is at mass coordinate zero), as shown in Figure 2. This mass range was shown by Wanajo et al. (2003) to give a solar-like $`r`$-abundance distribution assuming that part of the ejected material with low-$`Y_e`$ (layers 106 to 132 of Wanajo et al., 2003) was artificially assumed to be reaccreted onto the compact object, and consequently not to contribute to the galactic enrichment. The temperature and density histories of some mass elements are presented in Figure 3.

Adopting model Q6 in Wanajo et al. (2003) for the physical conditions, the $`r`$-process abundances are obtained by solving an extensive nuclear reaction network code. The network consists of $`4400`$ species, all the way from single neutrons and protons up to the fermium ($`Z=100`$) isotopes (Figure 4). We include all relevant reactions, i.e. $`(n,\gamma )`$, $`(p,\gamma )`$, $`(\alpha ,\gamma )`$, $`(p,n)`$, $`(\alpha ,n)`$, $`(\alpha ,p)`$, and their inverse. All reaction rates are calculated within the statistical model of Hauser-Feshbach making use of experimental masses (Audi & Wapstra, 2001) whenever available or the HFB-2 mass predictions (Goriely et al., 2002) otherwise. The photodisintegration rates are deduced from the reverse $`(n,\gamma )`$ rates applying the reciprocity theorem with the nuclear masses considered. The $`\beta `$-decay and $`\beta `$-delayed neutron emission rates are taken from the gross theory (GT2) of Tachibana, Yamada, & Yoshida (1990), obtained with the ETFSI (extended Thomas-Fermi plus Strutinsky integral, Aboussir et al., 1995) $`Q_\beta `$ predictions. Other nuclear inputs are the same as in Wanajo et al. (2003).

A word of caution is required concerning some inconsistencies in the nuclear input considered here, and most particularly concerning the use of ETFSI masses to estimate $`\beta `$-decay rates within the GT2 model. As mentioned in § 2, nuclear masses influence the $`r`$-abundance distribution mainly through the $`(n,\gamma )(\gamma ,n)`$ competition taking place in the neutron-rich region. Thanks to the balance theorem, such a competition remains relatively independent of the model adopted to estimate the neutron capture and photodisintegration rates. In contrast, the sensitivity of the absolute value of the reaction and $`\beta `$-decay rates on nuclear masses can depend significantly on the model considered. For example, within the statistical GT2 model, the $`\beta `$-decay rate scale as the fifth power of the $`Q_\beta `$ value, while in the continuum QRPA model (e.g Borzov & Goriely, 2003), the rates depend on the possible spin-isospin excitation within the $`Q_\beta `$ window but are not directly affected by the nuclear mass as such. Ultimately, for each mass model, a consistent estimate of the reaction and $`\beta `$-decay rates based on the corresponding mass model predictions for masses, but also deformation, single-particle properties, …should be performed. In order to simplify here the discussion on the impact of masses, but also to avoid the complicate model dependence in the estimate of the reaction and $`\beta `$-decay rates, we will restrict ourselves to consider in all cases only one set of $`\beta `$-decay rates, namely the GT2 rates with ETFSI masses.

The calculation in each mass trajectory is started at $`T_9=9`$ (where $`T_9T/10^9`$ K). The initial composition is taken to be that of the nuclear statistical equilibrium with the matter density $`\rho `$ at $`T_9=9`$, and consists mostly of free nucleons and alpha particles. The calculation is terminated after 10 seconds from the start in each mass trajectory. Snapshots of the $`r`$-process calculation in the mass-trajectory 99 ($`1.11M_{}`$ in mass coordinate and $`Y_{e0}=0.16`$) are shown in Figure 4. The abundances (mole fractions) are color coded, while the isotopes included in the reaction network are shown by dots with the stable and meta-stable isotopes by larger dots. Contours of constant two-neutron-separation energies divided by two ($`S_{2n}/2`$) from 1 to 8 MeV are superimposed. The abundance curve as a function of the mass number is also shown in each upper left panel.

The upper panel of Figure 4 shows the time slice when the neutron-to-seed ratio, $`Y_n/Y_h`$, decreases to $`1`$, where $`Y_n`$ is the abundance of free neutrons and

$$Y_h\underset{Z>2,A}{}Y(Z,A)$$
(1)
that of nuclei heavier than helium. At this time, the $`r`$-process flow closely follows an iso-$`S_{2n}`$ curve given by $`S_{2n}/2=S_a^01.6`$ MeV (red line), where

$$S_a^0(\mathrm{MeV})\left(34.075\mathrm{log}N_n+\frac{3}{2}\mathrm{log}T_9\right)\frac{T_9}{5.04},$$
(2)
$`N_n`$ is the number density of free neutrons in units of cm<sup>-3</sup>. The condition $`Y_n/Y_h=1`$ is referred to as the neutron exhaustion.

For the adopted $`(n,\gamma )`$ rates, a quasi-equilibrium between neutron capture and photodisintegration continues until $`\tau _\beta /\tau _n`$ decreases to $`1`$ as can be seen in the middle panel of Figure 4, where

$$\tau _\beta \left[\frac{1}{Y_h}\underset{Z>2,A}{}\lambda _\beta (Z,A)Y(Z,A)\right]^1$$
(3)
and

$$\tau _n\left[\frac{\rho Y_n}{Y_h}\underset{Z>2,A}{}\lambda _n(Z,A)Y(Z,A)\right]^1$$
(4)
are the averaged $`\beta `$-decay and neutron-capture lifetimes, and $`\lambda _\beta (Z,A)`$ and $`\lambda _n(Z,A)`$ the $`\beta `$-decay and neutron-capture rates of the nucleus $`(Z,A)`$, respectively. At this time, the abundance distribution is still near the $`S_a^02.6`$ MeV line. We refer hereafter the condition $`\tau _\beta /\tau _n=1`$ as freezeout, and the epoch between neutron exhaustion and freezeout as the freezeout phase. This freezeout phase is of particular relevance to determine the final abundance pattern as discussed in § 4. After the freezeout, the global abundance pattern is not significantly affected anymore, although it can still be somewhat smoothed out by photodisintegration and $`\beta `$-delayed neutron emission, as can be seen in the lower panel of Figure 4 ($`\tau _\beta /\tau _n0.2`$).

## 4 IMPACT OF MASS PREDICTIONS ON THE $`r`$-PROCESS

The final $`r`$-process yields are mass-averaged over all mass trajectories between 1.08 and $`1.30M_{}`$ in mass coordinate. In Figure 5a, the mass-averaged abundances (line) are compared with the solar $`r`$-process abundance pattern (dots, Goriely, 1999) that is scaled to match the height of the third $`r`$-process peak. For comparison, identical calculations were performed by replacing our standard HFB-2 masses by the HFB-7 (Goriely et al., 2003, Figure 5b), FRDM (Möller et al., 1995, Figure 5c), and DM (Hilf et al., 1976, Figure 5d) predictions.

Roughly speaking, an overall agreement between calculated and solar abundance curves can be seen in Figures 5a-d for nuclei with $`A=100200`$. Discrepancies common to all cases are, however, observed, in particular the deficiencies for nuclei in the vicinity of the first $`r`$-process peak ($`A=80`$), the shifts with respect to the solar curve of the second and third $`r`$-process peaks ($`A=130`$ and 195) to lower mass numbers ($`3`$ units in mass number) and the second to third peak height which is significantly smaller than observed in the solar distribution. We do not discuss further on the deficiencies near $`A=80`$, since theoretically the results are highly dependent on the adopted mass cut of the ejecta (see § 3 and Wanajo et al., 2003, for more details) and observationally there has been an increasing number of evidences that the light ($`A<130`$) $`r`$-process nuclei have a different origin from the heavy ($`A>130`$) ones (e.g., Hill et al., 2002; Sneden et al., 2003). For the shifts of the peak locations, as well as the second to third peak ratio, it is conceivable that this deficiency does not originate from the mass model, but rather from inadequate astrophysical or nuclear description, since it is observed for all mass formulas. This point will be further discussed in § 5.

A few significant differences in the abundance patterns can, however, be observed near the second and third peaks when use is made of the Hartree-Fock models (HFB-2 and HFB-7) on one side and the droplet models on the other side. First, the underproduction of nuclei at $`A115`$ and 140 is more pronounced with the FRDM (Figure 5c) and DM (Figure 5d) masses than with the HFB-2 (Figure 5a) or HFB-7 (Figure 5b) masses. Second, the abundances near $`A=130`$ in the HFB cases are spread out in contrast to what is observed in the solar $`r`$-abundances. Third, the abundance curves near the third peak with the HFB masses are widened and the valley at $`A=183`$ as observed in the solar $`r`$-distribution is significantly shifted to lower mass numbers.

These differences reflect the model properties of iso-$`S_{2n}/2`$ curves shown in Figure 6, along which the $`r`$-process proceeds (see also Figure 4). Major local differences between the HFB (Figures 6a, b) and the droplet masses (Figures 6c, d) are found near the neutron magic numbers $`N=82`$ and 126. The Hartree-Fock masses show weaker shell-closures, i.e. smoother iso-$`S_{2n}/2`$ curves, at $`N=82`$ and 126. This reduced shell effect (which actually is not totally quenched at the neutron dripline in contrast to other predictions, e.g Pearson, Nayak, & Goriely, 1996) is responsible for spreading the second and third abundance peaks. In contrast, the large droplet shell effect is clearly seen at $`N=82`$ and 126 by the steep $`S_{2n}/2`$ character even close to the neutron dripline and gives rise to sharp abundance peaks.

Figure 7 illustrates the formation process of the second and third peaks with the HFB-2 masses in more details. The left panels show the snapshots of the abundance distributions for the trajectory 90 ($`1.15M_{}`$ in mass coordinate and $`Y_{e0}=0.23`$), and the right panels for the trajectory 99 ($`1.11M_{}`$ in mass coordinate and $`Y_{e0}=0.16`$) at the neutron exhaustion (top panels), freezeout (middle panels), and later time (bottom panels). The trajectory 90 (left panels in Figure 7) gives a major contribution to the second-peak formation. The abundances are widely distributed already at neutron exhaustion (left-top panel), owing to the weaker shell gaps at $`N=82`$, in particular for low $`S_{2n}/2`$. This results in a split of the second abundance peak as seen in Figure 5a. The shift of the abundance distribution to the heavier mass number during the freezeout phase is small, although the abundance curve is significantly smoothed by photodisintegration (and slightly by $`\beta `$-delayed emission). The reason is that the freezeout takes place at a low $`S_a^0`$-value ($`2.9`$ MeV, left-middle panel in Figure 7).

The formation process of the third peak in the trajectory 99 can be seen in the right panels of Figure 7. The $`N=126`$ shell effect in the neutron-rich region is not reduced as the $`N=82`$ one, though it still remains weaker than in the droplet approach. The splitting of the peak abundances is not observed, but the abundance curve is significantly broaden by photodisintegrations during the freezeout phase. The freezeout takes place when $`S_a^0`$ is still low (2.6 MeV), and for this reason the peak position at $`A=190`$ is not shifted during the freezeout phase.

Some of the deficiencies observed in the $`r`$-abundance distribution obtained with HFB masses are seen to be cured when considering droplet masses (Figure 5). In particular, the DM masses characterized by isospin-independent $`N=82`$ and $`N=126`$ shell effects give rise to an abundance curve which is seen to be globally in better agreement with the solar pattern than do the HFB curves. It would be tempting to judge the quality of the nuclear mass predictions on the basis of such a comparison. However, it must be recalled here that such $`r`$-process calculations in prompt explosion scenarios (like in any other site studied so far) remain strongly affected by other input uncertainties in the astrophysical modeling as well as extra nuclear ingredients. More specifically, the total mass finally ejected (i.e not reaccreted onto the proto-neutron star) remains unknown and strongly dependent on the astrophysical modeling of the explosion. The same holds for the dynamical timescale of the outgoing material. Finally, major difficulties still exist in predicting reliably the $`\beta `$-decay half-lives of exotic neutron-rich nuclei. These uncertainties are studied in § 5, where it is shown that a modification of either the physical conditions of the explosion or the $`\beta `$-decay rates can influence drastically the agreement with the solar $`r`$-distribution.

Finally, it could also be argued that the outcome of the prompt explosion nucleosynthesis in one specific star might not need to match the solar $`r`$-process pattern. Recent spectroscopic studies show that the elemental abundance patterns between the second and third peaks of $`r`$-process-enhanced ultra-metal-poor stars in the Galactic halo are in excellent agreement with the solar $`r`$-process abundances. In particular, the agreement extends to the third peak elements Os, Ir, and Pt for HD 115444, HD 126238 (Sneden et al., 1998), BD +173248 (Cowan et al., 2002), and CS 22892-052 (Sneden et al., 2003, Figure 8), and CS 31082-001 (Hill et al., 2002). These observations led these authors to the conclusion that there is only one single $`r`$-process site that contributes to the production of heavy ($`A>130`$) $`r`$-process nuclei and that the corresponding distribution is similar to the solar $`r`$-process pattern (Sneden et al., 1996). Instead of comparing the predicted abundances with the solar content, we can therefore use such observations and estimate to what extent the disagreements seen in Figure 5 and discussed above are also reflected in the comparison with the observed elemental distribution. The comparison is shown in Figure 8 for the most-studied, extremely metal-poor star CS 22892-052 (Sneden et al., 2003). Even in this case, it is difficult to discriminate between the predictions obtained with the four different mass formulas considered here. This is quite normal since the observed distribution is reduced to some 30 elements (in comparison with some 130 nuclei in the solar material). And above all, the observed pattern could well be explained by the nuclear properties invariance in the $`56Z76`$ region rather than from specific astrophysics conditions as discussed in Goriely & Arnould (1997). It is therefore more informative to compare the $`r`$-abundance prediction with the solar pattern, keeping in mind that there is no stringent observational constraint that the final distribution should be strictly solar. Nevertheless, future abundance determinations of as many extremely metal-poor stars as possible will be of particular importance to shed light on nucleosynthesis yields from only one or a few supernova events.

## 5 UNCERTAINTIES IN THE $`r`$-PROCESS CALCULATIONS

In this section, we examine possible astrophysical and nuclear modifications leading to a better agreement of the HFB-2 abundance curve with the solar $`r`$-process pattern. Regarding the prompt supernova explosions considered here, we do not expect the entropy per baryon of the shocked material to be very different from the one deduced from the simulation ($`10`$ in units of the Boltzmann constant), since many previous simulations resulted in similar values (Hillebrandt, Nomoto, & Wolff, 1984; Sumiyoshi et al., 2001). However, the dynamical timescale of the outgoing matter (without changing the entropy) can significantly differ from model to model. For example, the dynamical timescales of the mass trajectories responsible for $`r`$-process in Sumiyoshi et al. (2001) are typically a few times shorter than those applied in this study (Wanajo et al., 2003), because of their different treatments of electron capture. It is also conceivable that other effects, such as convection, asymmetric explosion, or reverse shock from the outer envelope, cause an acceleration or a deceleration of the outgoing material.

In order to test the impact of a change in the dynamical timescales, we simply modify the density and temperature profiles of each trajectory, so that $`\rho ^{}(t)=\rho (t/f_t)`$ and $`T^{}(t)=T(t/f_t)`$, i.e the dynamical timescale is multiplied by a factor of $`f_t`$. Figures 9 and 10 show the abundance curves for $`f_t=0.2,0.5`$ (fast), 1 (standard), 2, 3, and 5 (slow) obtained with the mass-trajectories 90 and 99, which are responsible for the formation of the second and third peaks, respectively. The abundance curves near the second (Figure 9) and third (Figure 10) peaks are seen to agree fairly well with to the solar $`r`$-process pattern when the timescale is multiplied by a factor of $`23`$ (Figures 9d-e and Figures 10d-e), and to disagree for the fast trajectories (Figures 9a-b and Figures 10a-b).

This temporal effect is due to the significant abundance evolution during the feezeout phase for the slowly outgoing material, as shown in Figure 11 (same as Figure 7, but for $`f_t=3`$). At the neutron exhaustion, the abundances are dispersed near the second (left-top panel) and third (right-top panel) peaks, respectively, owing to the weak shell gaps at $`N=82`$ and 126 for low $`S_{2n}/2`$ ($`23`$ MeV), as in the standard case (Figure 7). The freezeout takes place (middle panels), however, at much higher $`S_a^0`$ values (4.2 and 3.7 MeV for the trajectories 90 and 99, respectively) than those in the standard case (2.9 and 2.6 MeV, Figure 7). The reason is that the temperatures at the freezeout are higher ($`T_9=1.4`$ and 1.2 for trajectories 90 and 99, respectively) owing to the slow expansion with $`f_t=3`$ than those in the standard case ($`T_9=1.0`$ and 0.93). As a result, the abundance peaks shift to $`A=130`$ (trajectory 90) and 194 (trajectory 99), i.e mostly at the same position as those in the solar $`r`$-distribution. Concomitantly, the width of the abundance curve near the second and third peaks is reduced significantly due to the larger shell gaps at $`S_a^04`$ MeV than that at $`S_a^03`$ MeV, and becomes similar to the one observed in the solar pattern (Figures 9e and 10e).

Assuming all the nuclear inputs are correct and the $`r`$-process indeed originates from the collapse of O-Ne-Mg core, this would imply that the expansion of the outgoing material must be significantly slower than obtained by Wanajo et al. (2003). A possible explanation could be found in the deceleration of the ejecta by the reverse shock from the outer envelope if the hydrogen and helium layers survive the mass loss till the onset of the core collapse. This effect is absent from the simulation of Wanajo et al. (2003) with a bare O-Ne-Mg white dwarf. This reverse shock effect would also reduce the total overproduction of $`r`$-processed material per event in the prompt explosion scenario (see Wanajo et al., 2003). It is interesting to note that such a significant deceleration of the outgoing material takes place in the neutrino-powered explosion due to the stalled shock wall once launched by a core bounce (Woosley et al., 1994). More detailed numerical simulations of the collapsing O-Ne-Mg core including the outer envelope are needed to quantify this effect.

The final mass-averaged abundance curve corresponding to mass trajectories slowed by a factor of three ($`f_t=3`$) is shown with the scaled solar $`r`$-process abundances in Figure 12a. We find a good agreement between the calculated and solar $`r`$-process patterns, in particular near the second and third peaks. The steep valley at $`A183`$ observed in the solar $`r`$-process abundances is also well reproduced. However, an underproduction at $`A115`$ and 140 appear.

As far as nuclear uncertainties are concerned, $`\beta `$-decays rates, neutron-captures rates and fission probabilities also play a key role (neutrino reactions are not important in the prompt explosion, see Wanajo et al., 2003). A change in their uncertain prediction can affect the nucleosynthesis predictions. In the present scenario, fission does not play an important role however, since it affects only a small fraction of the total $`r`$-processed matter (Wanajo et al., 2003). The prediction of neutron capture rates, especially for exotic neutron-rich nuclei, remains difficult, in particular in the prediction of the direct capture contribution and of the $`\gamma `$-ray strength and neutron-nucleus optical model potential for exotic neutron-rich nuclei as required in the statistical reaction model of Hauser-Feshbach (for more details, see Goriely, 1997, 2003). Much work remains to be done to estimate the neutron capture rates reliably. For this reason, their impact on the $`r`$-process nucleosynthesis is postponed to a future study.

Recent microscopic calculations of the $`\beta `$-decay rate within the energy density functional plus continuum QRPA formalism shows that the first forbidden transitions might have been underestimated in the past and can increase the $`\beta `$-decay rate by a factor of about 2 along the neutron-rich $`N=82`$ isotone and of about 3-10 along the neutron-rich $`N=126`$ isotone compared with the calculation based exclusively on the allowed transitions (Borzov, 2003). Detailed microscopic estimates of the $`\beta `$-decay rates (including forbidden transitions) for all nuclei of relevance in $`r`$-process calculations, as well as future experiments with radioactive ion beams are deeply awaited to improve this fundamental input in the description of the $`r`$-process nucleosynthesis.

To estimate the influence of $`\beta `$-decays, we show in Figure 12b the mass-averaged $`r`$-process yields obtained by multiplying all the $`\beta `$-decay rates by a factor of three (i.e., reducing $`\tau _\beta `$ by a factor of three). Interestingly, no significant differences are seen between Figures 12a and 12b. For $`\beta `$-decay rates faster by a factor of three, the freezeout (corresponding to $`\tau _\beta =\tau _n`$) takes place at higher temperatures and thus at higher $`S_a^0`$ value, which has globally the same effect as slowing down the outgoing material by the same factor. This test calculation shows the crucial role played by $`\beta `$-decay rates when considering non-parametrized density and temperature profiles, i.e thermodynamic conditions fixed by (semi-)realistic models.

## 6 SUMMARY

We have examined the $`r`$-process nucleosynthesis with the state-of-the-art microscopic mass formula, HFB-2 (Goriely et al., 2002), to see its effects on the theoretical $`r`$-process abundance prediction. The yields of $`r`$-process species were calculated with a nuclear reaction network code ($`4400`$ isotopes), adopting favorable physical conditions in the prompt explosion model of a collapsing O-Ne-Mg core by Wanajo et al. (2003). This model led to production of the solar-like $`r`$-abundance pattern in the previous study (Wanajo et al., 2003). However, an artificial enhancement of the shock-heating energy was needed to obtain requisite physical conditions for successful $`r`$-processing.

Due to its weak shell effect at the neutron magic numbers in the neutron-rich region, the microscopic mass formulas (HFB-2 and HFB-7) give rise to a spread of the abundance distribution in the vicinity of the $`r`$-process peaks ($`A=130`$ and 195). While this effect resolves the large underproduction at $`A115`$ and 140 obtained with droplet-type mass formulas, large deviations compared to the solar pattern are found near the third $`r`$-process peak. When use is made of droplet mass predictions, sharp $`r`$-process peaks are systematically found, owing to their strong shell effect for neutron magic numbers even in the neutron-rich region. However, due to the numerous uncertainties still affecting the astrophysics models as well as the prediction of extra nuclear ingredients, it would be highly premature to judge the quality of the mass formula on the basis of such a comparison.

We found that abundance peaks similar to the one observed in the solar system could be recovered if the dynamical timescales of the mass trajectories are increased by a factor of $`23`$ (without any change in the entropy) or by decreasing systematically the $`\beta `$-decay half-lives by the same factor. These changes might be conceivable when considering the current uncertainties in the astrophysics as well as in the nuclear $`\beta `$-decay model. Much effort in the astrophysics and nuclear modeling remain to be devoted to improve the difficult description of the $`r`$-process nucleosynthesis.

This work was supported by a Grant-in-Aid for Scientific Research (13640245, 13740129) from the Ministry of Education, Culture, Sports, Science, and Technology of Japan. M.S. and S.G. are FNRS Research Fellow and Associate, respectively.


text:
# Multifractality and nonextensivity at the edge of chaos of unimodal maps

## 1 Introduction

The thermodynamic framework for the characterization of multifractals  has proved to be a useful and insightful tool for the analysis of many complex system problems in which this geometrical feature plays a prominent role in establishing physical behavior. This has clearly been the case of dynamical systems, turbulence, growth models, sand pile models, etc. . As it has recently been pointed out, multifractality is also of relevance to critical  and glassy dynamics . In the former case the geometry of critical fluctuations is found to be multifractal and this has a hold on their dynamical properties , whereas in the latter case slow dynamics develops as phase-space accessibility is reduced and is finally confined to a multifractal subspace . Concurrently, in the study of the same types of systems the nonextensive generalization of the Boltzmann Gibbs (BG) statistical mechanics ,  has recently raised much interest and provoked considerable debate  as to whether there is firm evidence for its applicability in circumstances where a system is out of the range of validity of the canonical BG theory. Slow dynamics is plausibly related to hindrance of movement in phase space that leads to ergodicity breakdown, and this in turn has been linked to the failure of BG statistics . Here we briefly examine the relationship between multifractal properties and nonextensive dynamical evolution at the edge of chaos, a nonergodic state, in unimodal maps, as part of a wider analysis , that might contribute to reveal the physical source of the nonextensive statistics.

## 2 Nonextensivity of the Feigenbaum attractor

An important class of multifractals are the strange attractors that are generated by iteration of nonlinear maps. Amongst these, a prototypical example is the Feigenbaum attractor, also known as the edge of chaos, generated by the logistic map , $`f_{\mu ,2}(x)=1\mu x^2`$,$`1x1`$, when the control parameter $`\mu `$ takes the value $`\mu _{\mathrm{}}=1.40115\mathrm{}`$ . Recently, the predictions of the nonextensive theory have been rigorously proved   for this state that marks the transition from periodic to chaotic motion and that is characterized by the vanishing of the ordinary Lyapunov exponent $`\lambda _1=0`$. At the edge of chaos (and also at infinitely many other critical states of the same map such as the pitchfork and tangent bifurcations where $`\lambda _1=0`$ ) the sensitivity to initial conditions $`\xi _t\left|dx_t/dx_0\right|`$ (where $`x_t`$ is the orbit position at time $`t`$ given the initial position $`x_0`$ at time $`t=0`$) generally acquires power law instead of exponential time dependence. The nonextensive formalism suggests  that $`\xi _t`$ when $`\lambda _1=0`$ obeys the $`q`$-exponential expression,

$$\xi _t=\mathrm{exp}_q(\lambda _qt)[1(q1)\lambda _qt]^{1/(q1)},$$
(1)
containing the entropic index $`q`$ and the $`q`$-generalized Lyapunov coefficient $`\lambda _q`$. The BG exponential form $`\xi _t=\mathrm{exp}(\lambda _1t)`$ is recovered when $`q1`$. Further, it was conjectured  that the ordinary Pesin identity $`\lambda _1=K_1`$, where $`K_1`$ is the Kolmogorov-Sinai (KS) rate of entropy increase  is replaced by $`\lambda _q=K_q`$, where $`K_qt=S_q(t)S_q(0)`$, $`t`$ large, and where

$$S_q\underset{i}{}p_i\mathrm{ln}_q\left(\frac{1}{p_i}\right)=\frac{1_i^Wp_i^q}{q1}$$
(2)
is the Tsallis entropy ($`\mathrm{ln}_qy(y^{1q}1)/(1q)`$ is the inverse of $`\mathrm{exp}_q(y)`$). In the limit $`q1`$ $`K_q`$ becomes $`K_1t^1[S_1(t)S_1(0)]`$ where $`S_1(t)=_{i=1}^Wp_i(t)\mathrm{ln}p_i(t)`$. In Eq. (2) $`p_i(t)`$ is the probability distribution obtained from the relative frequencies with which the positions of an ensemble of trajectories occur within cells $`i=1,\mathrm{},W`$ at iteration time $`t`$. Based on the renormalization group (RG) transformation $`𝐑f(x)\alpha f(f(x/\alpha ))`$ where $`\alpha 2.5029`$ is the Feigenbaum’s universal constant that measures the power-law period-doubling spreading of iterate positions, a rigorous analytical confirmation of Eq. (1) and also of $`\lambda _q=K_q`$ have been obtained  . Specifically, $`\lambda _q`$ and $`q`$ are simply given by $`\lambda _q=\mathrm{ln}\alpha /\mathrm{ln}21.3236`$ and $`q=1\mathrm{ln}2/\mathrm{ln}\alpha 0.2445`$.

On the other hand, a remarkable link between the well-known  geometrical description of multifractals and the nonextensive statistics has been derived and corroborated numerically not only for the Feigenbaum attractor but for the edge of chaos of the generalization of the logistic map to nonlinearity $`z>1`$, $`f_{\mu ,z}(x)=1\mu \left|x\right|^z`$,$`1x1`$ . This consists of the expression

$$\frac{1}{1q}=\frac{1}{D_{\mathrm{}}}\frac{1}{D_{\mathrm{}}},$$
(3)
where $`D_{\mathrm{}}`$ and $`D_{\mathrm{}}`$ are the extreme values of the generalized dimension spectrum $`D_\beta `$, and correspond to the smallest and largest length scales of the multifractal set . For the $`z`$-logistic map one has $`D_{\mathrm{}}=\mathrm{ln}2/z\mathrm{ln}\alpha `$ and $`D_{\mathrm{}}=\mathrm{ln}2/\mathrm{ln}\alpha `$ (notice that $`\alpha `$ depends on $`z`$) . In view of the clear association of the nonextensive statistics with the properties at the edge of chaos of the logistic map we have pursued the study of this critical state by means of the standard thermodynamic tools that characterize multifractals. Here we advance some of the results we have found .

## 3 Sensitivity to initial conditions

We determine the index $`q`$ and the $`q`$-Lyapunov coefficient $`\lambda _q`$ at the edge of chaos of the $`z`$-logistic maps through direct evaluation of the sensitivity to initial conditions $`\xi _t`$. By making use of the properties  $`g(0)=(\alpha )^ng^{(2^n)}(0)=1`$, $`g(1)=(\alpha )^ng^{(2^n)}(\alpha ^n)=\alpha ^1`$ and

$$g^{}(1)=(1)^ng^{}(\alpha ^n)g^{}(g(\alpha ^n))\mathrm{}g^{}(g^{(2^n1)}(\alpha ^n))=\alpha ^{z1},$$
(4)
of the fixed-point map, $`g(x)=lim_n\mathrm{}(\alpha )^nf_{\mu ,z}^{(2^n)}(x/\alpha ^n)`$, into the definition of $`\xi _t`$,

$$\xi _t\left|\frac{dx_t}{dx_0}\right|=\left|\frac{dg^{(2^n1)}(x)}{dx}\right|_{x=x_0},$$
(5)
with the choice $`x_0=0`$, and for iteration times of the form $`t=2^n1`$, we obtain $`\xi _t=\alpha ^{(z1)n}`$. This result can be re-expressed as

$$\xi _t=\mathrm{exp}_q(\lambda _qt)[1(q1)\lambda _qt]^{1/(q1)},t=2^n1,$$
(6)
with

$$\lambda _q=\frac{1}{1q}=(z1)\frac{\mathrm{ln}\alpha }{\mathrm{ln}2}.$$
(7)
For details see Ref. . Thus, without reference to the multifractal properties of the attractor, we have corroborated the result in Eq. (3)  and obtained additionally the $`q`$-Lyapunov coefficient $`\lambda _q=(1q)^1`$. Alternatively, $`\xi _t`$ can be obtained from $`\left|\mathrm{\Delta }x_t\right|/\left|\mathrm{\Delta }x_0\right|`$ in the limit $`\mathrm{\Delta }x_00`$ where $`\mathrm{\Delta }x_t`$ is the distance between two orbits at time $`t`$ when they were a distance $`\mathrm{\Delta }x_0`$ at $`t=0`$. Following the procedure described in Ref.  with the choices $`\left|\mathrm{\Delta }x_0\right|=\alpha ^j\alpha ^i`$, $`i,j>0`$ and $`t=2^n1`$, $`n=0,1,\mathrm{}`$, we obtain $`\left|\mathrm{\Delta }x_t\right|=`$ $`\left|\mathrm{\Delta }x_0\right|\alpha ^{(z1)n}`$ . Next we study a multifractal property of the attractor and make contact again with Eq. (3).

## 4 Scaling properties at the edge of chaos

In Fig. 1a we have plotted (in logarithmic scales) the first few (absolute values) of iterated positions $`\left|x_t\right|`$ of the orbit at $`\mu _{\mathrm{}}`$ starting at $`x_0=0`$ where the labels indicate iteration time $`t`$. Notice the structure of horizontal bands, and that, in the top band lie half of the attractor positions (odd times), in the second band a quarter of the attractor positions, and so on. The top band is eliminated by functional composition of the original map, that is by considering the orbit generated by the map $`g^{(2)}(0)`$ instead of $`g(0)`$. Successive bands are eliminated by considering the orbits of $`g^{(2^k)}(0),k=1,2\mathrm{}`$. The top band positions (odd times) can be reproduced approximately by the positions of the band below it (times of the form $`t=2+4n`$, $`n=0,1,2\mathrm{}`$) by multiplication of a factor equal to $`\alpha `$, e.g. $`\left|x_1\right|\alpha \left|x_2\right|`$. Likewise, the positions of the second band are reproduced by the positions of the third band under multiplication by $`\alpha `$, e.g. $`\left|x_2\right|\alpha \left|x_4\right|`$, and so on. In Fig. 1b we show a logarithmic scale plot of $`1`$ $`\left|x_t\right|`$ that displays a band structure similar to that in Fig. 1a, in the top band lie again half of the attractor positions (even times) and below the other half (odd times) is distributed in the subsequent bands. This time the positions in one band are reproduced approximately from the positions of the band lying below it by multiplication of a factor $`\alpha ^z`$, e.g. $`1\left|x_3\right|\alpha ^z(1\left|x_5\right|)`$. The scaling properties amongst iterate positions merely follow from repeated composition of the map (for more details see Ref. ) and suggest the construction of an RG transformation described below. For this purpose it is convenient to write the two scale factors $`\alpha `$ and $`\alpha ^z`$ as $`2^\eta _{\mathrm{}}`$ and $`2^\eta _{\mathrm{}}`$, respectively, where $`\eta _{\mathrm{}}=\mathrm{ln}\alpha /\mathrm{ln}2`$ and $`\eta _{\mathrm{}}=z\mathrm{ln}\alpha /\mathrm{ln}2`$.

## 5 ’Take away half’ RG transformation

Consider the positions of the orbit at $`\mu _{\mathrm{}}`$ starting at $`x_0=0`$ up to time $`t=2^N1`$ and from them extract the set of lengths $`l_j^{(0)},j=1\mathrm{},2^N`$, that join closest adjacent positions. Each position is used only once so that the set $`l_j^{(0)}`$ covers the attractor but not all the interval $`1x1`$. To each length $`l_j^{(0)}`$ assign the equiprobability $`p_j^{(0)}=2^N`$. Following the standard procedure  we define the partition function

$$Z_N^{(0)}(\beta ,\tau )=\underset{j}{}[p_j^{(0)}]^\beta [l_j^{(0)}]^\tau ,$$
(8)
and determine the (finite size) generalized dimensions $`D_{\beta ,N}\tau /(\beta 1)`$ by requiring $`Z_N^{(0)}(\beta ,\tau )=1`$, ($`D_\beta =lim_N\mathrm{}D_{\beta ,N}`$). In Fig. 2 we show results for $`D_{\beta ,N}`$ for several values of $`z`$.

We introduce now repeated length scale and restoring probability transformations

$$l_j^{(k)}=2^{k\eta }l_j^{(0)}\text{and}p_j^{(k)}=2^kp_j^{(0)},k=1\mathrm{},2^N,$$
(9)
on the partition function $`Z_N^{(0)}(\beta ,\tau )`$. That is, to transform $`Z_N^{(0)}(\beta ,\tau )`$ into $`Z_N^{(1)}(\beta ,\tau )`$ replace $`l_j^{(0)}`$ and $`p_j^{(0)}`$ in Eq. (8) by $`l_j^{(1)}=2^\eta l_j^{(0)}`$ and $`p_j^{(1)}=2^1p_j^{(0)}`$, respectively, and proceed similarly any number $`k`$ of times. In the limit $`N\mathrm{}`$ the transformation eliminates one half of the multifractal points each time it is applied. By requiring $`Z_N^{(k)}(\beta ,\tau )=1`$ we obtain, again when $`N\mathrm{}`$, the following result for the rescaling exponent $`\eta `$ of the interval lengths

$$\eta =\frac{\beta }{\tau }=\frac{\beta }{\beta 1}\frac{1}{D_\beta }.$$
(10)
As shown in Fig. 3, $`\eta `$ approaches two limiting values, $`\eta _{\mathrm{}}=D_{\mathrm{}}^1`$ and $`\eta _{\mathrm{}}=D_{\mathrm{}}^1`$ as $`\beta \mathrm{}`$ and $`\beta \mathrm{}`$, respectively, and exhibits a singularity at $`\beta =1`$. Except for values of $`\beta `$ in the vicinity of $`\beta =1`$ the exponent $`\eta `$ is always close to either $`\eta _{\mathrm{}}`$ or $`\eta _{\mathrm{}}`$ which correspond to the scaling factors $`l_j^{(k)}=\alpha ^kl_j^{(0)}`$ and $`l_j^{(k)}=\alpha ^{zk}l_j^{(0)}`$ described in the previous section. As also shown in Fig. 3, there appears a gap in the permissible values of $`\eta `$,

$$\mathrm{\Delta }\eta \eta _{\mathrm{}}\eta _{\mathrm{}}=\frac{1}{D_{\mathrm{}}}\frac{1}{D_{\mathrm{}}},$$
(11)
that is noticeably identical to the dynamical result Eq. (7), i.e. $`\lambda _q=(1q)^1=\mathrm{\Delta }\eta `$.

Repeated application of the RG transformation when the limit $`N\mathrm{}`$ is taken before $`k\mathrm{}`$ leaves the generalized dimensions invariant, i.e. $`D_\beta =D_\beta ^{(k)}`$ for all $`k`$. However when both $`k,N\mathrm{}`$ with $`k=N`$ the RG transformation leads to a fixed point $`D_\beta ^{(k)}D_\beta ^{(\mathrm{})}`$ that has the form of a step function, $`D_\beta ^{(\mathrm{})}=\mathrm{ln}2/2\mathrm{ln}\alpha `$ for $`\beta <1`$ and $`D_\beta ^{(\mathrm{})}=\mathrm{ln}2/2z\mathrm{ln}\alpha `$ for $`\beta >1`$, with a divergence at $`\beta =1`$. See Fig. 4. .

## 6 Summary and discussion

Assisted by the known properties of the fixed-point map $`g(x)`$ we derived the expression for the sensitivity to initial conditions $`\xi _t`$ for the $`z`$-logistic map at the edge of chaos. As a consequence we corroborated the known expression for the entropic index $`q`$ in terms of $`\alpha (z)`$  and obtained in addition the closely related expression for the $`q`$-Lyapunov coefficient $`\lambda _q`$ (known previously only for the particular case $`z=2`$  ). Our results are rigorous. The expression obtained is precisely that first derived, heuristically, by Lyra and Tsallis , Eq. (3), that relates $`q`$ to the smallest $`D_{\mathrm{}}`$ and largest $`D_{\mathrm{}}`$ generalized dimensions of a multifractal set. In order to probe this connection between dynamic and geometric properties of the multifractal attractor we introduced an RG approach that consists of sequentially removing remaining halves $`2^k,k=1,2,\mathrm{}`$ of points of the multifractal set while keeping the value of its partition function $`Z_N^{(k)}(\beta ,\tau )`$ fixed. We determined the length rescaling index $`\eta `$ as a function of $`\beta `$ and observed that there is a gap $`\mathrm{\Delta }\eta `$ of permissible values in this quantity equal to $`\lambda _q=1/(1q)`$, as the edges of the gap correspond to $`\eta _{\mathrm{}}=D_{\mathrm{}}^1`$ and $`\eta _{\mathrm{}}=D_{\mathrm{}}^1`$. This result can perhaps be better appreciated if we note that the orbit positions $`\left|x_t\right|`$ evolve according to two basic alternating movements observable at time intervals of the form $`2^n`$, one drives the iterate towards $`\left|x\right|=1`$ as $`\alpha ^{zn}=2^{\eta _{\mathrm{}}n}`$ and the other one towards $`x=0`$ as $`\alpha ^n=2^{\eta _{\mathrm{}}n}`$. See the main diagonal subsequence positions in Figs. 1a and 1b. The net growth in the distance between two orbits is the ratio of these two power laws, that is $`\left|\mathrm{\Delta }x_t\right|`$ $`=`$ $`\left|\mathrm{\Delta }x_0\right|2^{\mathrm{\Delta }\eta n}`$ (see Section 3). Thus, Eq. (3) is obtained as the result of orbits that successively separate and converge, respectively, in relation to $`D_{\mathrm{}}`$ and $`D_{\mathrm{}}`$. The power law for $`\xi _t`$ can only persist if $`D_{\mathrm{}}D_{\mathrm{}}`$.

Acknowledgments. We were partially supported by CONACyT grant P40530-F (Mexican agency).